0 引言
地下流体注采诱导变形理论体系的构建始于Terzaghi[6]提出的土体固结数学模型，随后Biot[7]建立的经典三维固结理论奠定了现代流固耦合研究的理论基础，但其未明确界定渗流与变形边界条件。1957年Geertsma[8]提出油藏生产过程中地层上覆压力不变与单轴变形的假设，完善了Biot理论在处理油藏开发问题时的边界条件。1973年Geertsma[9]进一步结合Mindlin影响函数[10]，通过对油藏边界进行积分首次建立了均匀压降油藏地表位移解析解，但其仍受限于理想化的几何构型与弹性介质假设。为突破解析方法在复杂地质条件下的应用瓶颈，数值模拟技术逐步成为研究主流。Morita等[11]假设油藏中压降恒定，推导了油藏总应力和有效应力的有限元计算模型。此后大量学者采用Abaqus[12-13]、FLAC3D[14]等有限元力学求解器与Eclipse等油藏数值模拟器迭代耦合实现流固耦合求解，也开发了许多全耦合有限元求解模型[15-16]。21世纪以来，非常规油气的大规模开发推动了研究向复杂裂缝性储层的转变。研究人员围绕等效连续介质[17]、离散裂缝[18⇓-20]、嵌入式离散裂缝[21⇓-23]等模型开展研究，在有限元基础上发展形成了边界元法（BEM）[24-25]、扩展有限元法（XFEM）[26⇓-28]、有限体积法（FVM）[29]、间断有限元法[30]、混合有限元法[31-32]等多种创新算法。尽管这些方法在裂缝局部力学场的精确表征方面取得显著进步，但在油藏变形边界处理上仍沿袭传统理念：将油藏范围作为渗流与变形计算的区域，油藏外部区域简化为油藏边界上恒定的应力或位移条件[33-34]，或者从油藏顶部、底部、侧边分别向外进行有限延伸，扩展变形计算区域后在边界赋予应力或位移条件。然而，开发诱导变形与整个半无限大地层密切相关，这种边界简化会导致计算结果出现偏差，吴大卫等[35]研究发现变形计算区域的大小会显著影响计算结果，因此，目前主流商业软件通常采用扩展变形计算区域的方法，但为了保证计算的可行性，在油藏外会使用比油藏网格大得多的网格单元[36]。
目前主流的有限元开发诱导变形求解方法存在以下主要缺陷：无法完整考虑半无限大地层；渗流边界对地层变形的影响未被充分考量；远离油藏区域的网格过大导致远场变形计算精度较低[37]。为此，本文从半无限大地层视角出发，通过详细力学分析揭示油藏内部渗流与渗流边界对地层变形的作用机理，在此基础上，建立了油气藏开发诱导变形的等效力模型，并提出体积边界元数值求解方法，实现全空间域精确计算，为油气藏开发流固耦合问题准确高效数值模拟奠定了基础。
1 基本方程
取地层内任意微元，其体积相对于孔隙结构足够大、相对于整个地层足够小，将其视为均质材料。研究过程中要保证微元的性质与其在地层中一致，例如微元表面的孔隙与岩石骨架并非同时作用在一个封闭面上（除单独讨论的封闭边界面外），而是孔隙与微元外的孔隙对应，岩石骨架与微元外的岩石骨架对应，两个系统分别对应连接。
应力是作用于单位面积上的内力，即内力的分布集度，因此应力的定义必须明确其对应的面积。本文定义岩石骨架应力对应的面积为岩石表观面积，而非骨架面积。真实骨架应力对应的面积为骨架面积，则本文定义岩石骨架应力与真实骨架应力的关系为：
${{\sigma }_{ij}}=\left( 1-\phi \right){{\sigma }_{\text{s},ij}}$（i=x，y，z；j=x，y，z）
基于上述定义，岩石内任意截面的总应力为[38]：
${{\sigma }_{\text{t},ij}}={{\sigma }_{ij}}-\phi p{{\delta }_{ij}}$
（2）式中，孔隙流体压力为压应力，故流体应力为$-p$。
岩石骨架的平衡方程为[39]：
$\frac{\partial {{\sigma }_{ij}}}{\partial {{x}_{j}}}+{{f}_{i}}=0$
岩石骨架的几何方程为[39]：
${{\varepsilon }_{ij}}=\frac{1}{2}\left( \frac{\partial {{u}_{i}}}{\partial {{x}_{j}}}+\frac{\partial {{u}_{j}}}{\partial {{x}_{i}}} \right)$
线弹性各向同性多孔介质的本构方程为[39]：
${{\sigma }_{\text{t},ij}}=2G{{\varepsilon }_{ij}}+\frac{2G\nu }{1-2\nu }{{\delta }_{ij}}\varepsilon -\left( 1-\frac{{{K}_{p}}}{{{K}_{s}}} \right){{\delta }_{ij}}p$
其中 $\varepsilon ={{\varepsilon }_{xx}}+{{\varepsilon }_{yy}}+{{\varepsilon }_{zz}}$
将（2）式代入（5）式可得：
${{\sigma }_{ij}}=2G{{\varepsilon }_{ij}}+\frac{2G\nu }{1-2\nu }{{\delta }_{ij}}\varepsilon -\left( 1-\phi -\frac{{{K}_{p}}}{{{K}_{s}}} \right){{\delta }_{ij}}p$
2 开发诱导变形等效力模型
2.1 渗流等效力
当地下流体注采引起油藏孔隙压力变化时，岩石骨架应力也会同步变化。由于流体在岩石孔喉中的渗流速度极低，流体可以被认为是准静态的，即处于受力平衡状态。同样，岩石骨架也处于受力平衡的准静态，流体与骨架的受力状态如图1 所示。
根据图1b ，可以建立流体的受力平衡方程：
$A\left[ {{\left. (p\phi ) \right|}_{{{x}_{i}}}}-{{\left. (p\phi ) \right|}_{{{x}_{i}}+\Delta {{x}_{i}}}} \right]+{{\rho }_{f}}gA\Delta {{x}_{i}}\phi \frac{{{\left. D \right|}_{{{x}_{i}}+\Delta {{x}_{i}}}}-{{\left. D \right|}_{{{x}_{i}}}}}{\Delta {{x}_{i}}}-$
${{F}_{i}}A\Delta {{x}_{i}}=0$
（7）式两端同时除以AΔxi，取极限Δxi→0，可得流体作用在骨架上的体积力为：
${{F}_{i}}=-\frac{\partial p\phi }{\partial {{x}_{i}}}+{{\rho }_{f}}g\phi \frac{\partial D}{\partial {{x}_{i}}}$
同理，根据图1c ，可以得到单位表观体积岩石骨架所受的总外力为：
${{f}_{i}}={{F}_{i}}+{{\rho }_{\text{s}}}g\left( 1-\phi \right)\frac{\partial D}{\partial {{x}_{i}}}$
将（8）式代入（9）式，得到：
${{f}_{i}}=-\frac{\partial p\phi }{\partial {{x}_{i}}}+\bar{\rho }g\frac{\partial D}{\partial {{x}_{i}}}$
其中 $\bar{\rho }={{\rho }_{f}}\phi +{{\rho }_{\text{s}}}\left( 1-\phi \right)$
将（10）式代入（3）式，得到骨架应力平衡方程：
$\frac{\partial {{\sigma }_{ij}}}{\partial {{x}_{j}}}-\frac{\partial p\phi }{\partial {{x}_{i}}}+\frac{\partial D}{\partial {{x}_{i}}}\bar{\rho }g=0$
将（4）式、（6）式代入（11）式，可以得到用位移表示的骨架变形微分方程组：
$G{{u}_{i,jj}}+\frac{G}{1-2\nu }{{u}_{j,ij}}-\left( 1-\frac{{{K}_{p}}}{{{K}_{s}}} \right)\frac{\partial p}{\partial {{x}_{i}}}+\frac{\partial D}{\partial {{x}_{i}}}\bar{\rho }g=0$
（12）式与弹性力学中位移法求解的Lame-Navier方程（见（13）式）在形式上完全一致[40]。
$G{{u}_{i,jj}}+\frac{G}{1-2\nu }{{u}_{j,ij}}+{{f}_{\text{e},i}}=0$
其中，fe,i是i方向的渗流等效力，可表示为：
${{f}_{\text{e},i}}=-\left( 1-\frac{{{K}_{p}}}{{{K}_{s}}} \right)\frac{\partial p}{\partial {{x}_{i}}}+\frac{\partial D}{\partial {{x}_{i}}}\bar{\rho }g\text{ }$
由此可见，流体渗流引发多孔介质骨架变形的规律相当于在（14）式渗流等效力作用下的一般弹性变形，该渗流等效力是一种体积力。
在油气田开发领域，更受关注的问题是岩石形态的变化量，相对于初始状态的渗流等效力增量为：
$\Delta f_{\mathrm{e}, i}=-\left(1-\frac{K_{\mathrm{p}}}{K_{\mathrm{s}}}\right) \frac{\partial \Delta p}{\partial x_{i}}+\frac{\partial D}{\partial x_{i}} \Delta \bar{\rho} g \quad\left(x_{i} \in \Omega\right)$
在地下流体注采过程中，流体密度随孔隙压力和流体组成变化较大，而任意微元内岩石骨架质量基本不变。因此，忽略$\Delta \left[ {{\rho }_{s}}(1-\phi ) \right]$，（15）式可简化为：
$\Delta f_{\mathrm{e}, i}=-\left(1-\frac{K_{\mathrm{p}}}{K_{\mathrm{s}}}\right) \frac{\partial \Delta p}{\partial x_{i}}+\frac{\partial D}{\partial x_{i}} \Delta \rho_{\mathrm{f}} g \phi \quad\left(x_{i} \in \Omega\right)$
2.2 边界等效力
油气藏开发诱导地层变形问题的边界条件包括两种，即地层岩石变形边界和流体渗流边界。对于实际地层，完整的变形边界条件应包括地表的自由边界和整个地壳的半无限大边界。通常油藏及其周缘地层弹性参数的变化并不显著，为了简化计算和分析，本文假设地层岩石为均质弹性体，据此可将问题简化为均质半无限自由边界面弹性问题。
流体渗流边界为油藏的封闭边界，封闭边界两侧法向总应力平衡（见图2 ），即：
${{\mathbf{n}}_{\Omega }}$和 分别是封闭边界内部空间和外部空间在封闭边界Γ处的单位内法向向量，并满足 。爱因斯坦求和惯例适用于（17）式—（20）式中循环的拉丁下标。
（17）式可以改写为：
将（4）式和（6）式代入（18）式，可得：
（19）式与弹性力学中以位移表示的界面法向力连续性条件（见（20）式）在形式上完全一致[40]。
其中，pe,i是i方向的边界等效力，可表示为：
由此可见，油藏封闭边界对地层变形的影响相当于在封闭边界虚拟表面上施加边界等效力的一般弹性变形，该边界等效力是一种面积力。
相对于初始状态的边界等效力增量为：
油藏封闭边界外部不发生流体注采，封闭边界外部孔隙压力变化可以忽略（即 ），因此（22）式可简化为：
$\Delta p_{\mathrm{e}, i}=-\left(1-K_{\mathrm{p}} / K_{\mathrm{s}}\right) \Delta p_{\Omega} n_{Q_{\Omega, i}} \quad\left(x_{i} \in \Gamma\right)$
以上研究表明油藏内部渗流与渗流边界对地层变形的影响分别等效于油藏内部和油藏边界处受到外力作用时对地层变形产生的影响。与油藏内部渗流等效的外力为渗流等效力，与渗流边界等效的外力为边界等效力。
2.3 解析求解方法
对于实际地层，完整的变形边界可视为整个地壳的半无限大边界。结合前文的分析，油藏开发引起的地层变形，等价于在油藏内部渗流等效力和油藏封闭边界处边界等效力作用下的半无限大一般弹性变形。因此，对于半无限大变形边界，将Mindlin[10]的半无限大体内部集中力基本解作为格林函数，油藏内部渗流和渗流边界分别转换为半无限大地层受到的渗流等效力和边界等效力两种外力条件。根据叠加原理，地层任意点的变形解可表示为渗流等效力、边界等效力与格林函数的卷积。
$\left\{ \begin{align} & {{{\tilde{u}}}_{z\_x}}=\frac{\alpha -x}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left[ \frac{\gamma -z}{{{R}_{1}}^{3}}+\frac{\left( 3-4\nu \right)\left( \gamma -z \right)}{{{R}_{2}}^{3}}- \right. \\ & \left. \begin{matrix} {} \\ {} \\\end{matrix}\frac{4\left( 1-\nu \right)\left( 1-2\nu \right)}{{{R}_{2}}\left( {{R}_{2}}+\gamma +z \right)}+\frac{6\gamma z\left( \gamma +z \right)}{{{R}_{2}}^{5}} \right] \\ & {{{\tilde{u}}}_{z\_y}}=\frac{\beta -y}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left[ \frac{\gamma -z}{{{R}_{1}}^{3}}+\frac{\left( 3-4\nu \right)\left( \gamma -z \right)}{{{R}_{2}}^{3}}- \right. \\ & \left. \begin{matrix} {} \\ {} \\\end{matrix}\frac{4\left( 1-\nu \right)\left( 1-2\nu \right)}{{{R}_{2}}\left( {{R}_{2}}+\gamma +z \right)}+\frac{6\gamma z\left( \gamma +z \right)}{{{R}_{2}}^{5}} \right] \\ & {{{\tilde{u}}}_{z\_z}}=\frac{1}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left[ \frac{3-4\nu }{{{R}_{1}}}+\frac{8{{\left( 1-\nu \right)}^{2}}-\left( 3-4\nu \right)}{{{R}_{2}}}+ \right. \\ & \left. \begin{matrix} {} \\ {} \\\end{matrix}\frac{{{\left( \gamma -z \right)}^{2}}}{{{R}_{1}}^{3}}+\frac{\left( 3-4\nu \right){{\left( \gamma +z \right)}^{2}}-2\gamma z}{{{R}_{2}}^{3}}+\frac{6\gamma z{{\left( \gamma +z \right)}^{2}}}{{{R}_{2}}^{5}} \right] \end{align} \right.$
$\left\{ \begin{align} & {{{\tilde{u}}}_{x\_x}}=\frac{1}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left\{ \frac{3-4\nu }{{{R}_{1}}}+\frac{1}{{{R}_{2}}}+\frac{{{\left( \alpha -x \right)}^{2}}}{{{R}_{1}}^{3}}+\frac{\left( 3-4\nu \right){{\left( \alpha -x \right)}^{2}}}{{{R}_{2}}^{3}}+ \right. \\ & \begin{matrix} {} \\ {} \\\end{matrix}\frac{2\gamma z}{{{R}_{2}}^{3}}\left[ 1-\frac{3{{\left( \alpha -x \right)}^{2}}}{{{R}_{2}}^{2}} \right]+\left. \frac{4\left( 1-\nu \right)\left( 1-2\nu \right)}{{{R}_{2}}+\gamma +z}\left[ 1-\frac{{{\left( \alpha -x \right)}^{2}}}{{{R}_{2}}\left( {{R}_{2}}+\gamma +z \right)} \right] \right\} \\ & {{{\tilde{u}}}_{x\_y}}=\frac{\left( \alpha -x \right)\left( \beta -y \right)}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left[ \frac{1}{{{R}_{1}}^{3}}+\frac{3-4\nu }{{{R}_{2}}^{3}}-\frac{6\gamma z}{{{R}_{2}}^{5}}-\frac{4\left( 1-\nu \right)\left( 1-2\nu \right)}{{{R}_{2}}{{\left( {{R}_{2}}+\gamma +z \right)}^{2}}} \right] \\ & {{{\tilde{u}}}_{x\_z}}=\frac{\alpha -x}{16\text{ }\!\!\pi\!\!\text{ }G\left( 1-\nu \right)}\left[ \frac{\gamma -z}{{{R}_{1}}^{3}}+\frac{\left( 3-4\nu \right)\left( \gamma -z \right)}{{{R}_{2}}^{3}}-\frac{6\gamma z\left( \gamma +z \right)}{{{R}_{2}}^{5}}+ \right. \\ & \left. \begin{matrix} {} \\ {} \\\end{matrix}\frac{4\left( 1-\nu \right)\left( 1-2\nu \right)}{{{R}_{2}}\left( {{R}_{2}}+\gamma +z \right)} \right] \end{align} \right.$
其中
将以上单位力造成的位移定义为位移格林函数 ，表示在点(x, y, z)处受到i方向单位力作用造成点(α, β, γ)处j方向的位移量。
根据线弹性叠加原理，油藏开发引起的地层相对于初始状态的位移可以表示为等效力增量与位移格林函数的卷积：
根据位移增量，将（27）式代入（4）式即可得到相对于注采初始状态的表观应变，再将应变代入（6）式即可得到相对于注采初始状态的岩石骨架应力。
2.4 体积边界元数值求解方法
实际油藏研究的难点在于油藏形状不规则，孔隙压力变化不均匀，只能通过离散数值模拟得到。因此，有必要对上述模型进行离散化处理。
对于渗流等效力，采用中心差分的方式（油藏边界处采用单侧差分）可将（16）式离散表示为：
对于边界等效力，非边界网格均为零，边界网格根据（23）式可离散表示为：
对于位移变化量，（27）式可离散表示为：
其中 ，为网格体积源， 为网格边界源：
（31）式中对油藏边界网格的格林函数进行面积积分得到了网格边界源，对油藏内部网格的格林函数进行体积积分得到了网格体积源，地层任意点的变形解即可根据（30）式，将网格边界等效力与对应的网格边界源相乘、网格渗流等效力与对应的网格体积源相乘再累加得到。
以上数值求解方法与边界元法思路相似，但不同之处在于本文（30）式所示的积分方程除了要对存在边界等效力（面积力）的油藏渗流边界进行面积积分外，还需对存在渗流等效力（体积力）的油藏内部区域进行体积积分，因此将本文建立的这种积分求解方法称为体积边界元法。
3 模型验证与测试
利用特定时间与初始状态的油藏属性场直接进行变形估计，可以有效反映变形求解器的求解精度。本文首先采用还原Geertsma[9]研究的方法对（27）式所示解析解进行验证，随后采用枯竭油藏变形计算与商业数值求解器VISAGE进行对比测试。
3.1 解析解验证
Geertsma[9]建立了薄层圆形油藏均匀压降时地层沉降的解析解，油藏参数如表1 所示。
表1 用于解析解验证的油藏参数 |
| 参数名称 | 参数值 | 参数名称 | 参数值 |
|---|---|---|---|
| 骨架材料体积模量 | 35 000 MPa | 油藏厚度 | 10 m |
| 岩石表观体积模量 | 2 000 MPa | 油藏顶深 | 1 000 m |
| 岩石表观剪切模量 | 1 200 MPa | 油藏半径 | 1 000 m |
| 岩石表观泊松比 | 0.25 | 原始地层压力 | 30 MPa |
| 岩石孔隙度 | 20% | 油藏压力增量 | -10 MPa |
| 流体密度 （参考压力） | 800 kg/m3 （30 MPa） | 综合压缩系数 | 3×10-3 MPa-1 |
考虑状态方程为：
根据（16）式，由于油藏中的压力是均匀分布的，可得渗流等效力为：
根据（23）式，油藏顶面的边界等效力为：
同理，油藏底面的边界等效力为：
油藏侧面的边界等效力可表示为：
3.2 商业模拟器对比测试
考虑图5 所示的均质方形油藏从初始状态衰竭开发至枯竭状态，模型参数如表2 所示，弹性参数与其他物性参数同表1 。
VISAGE地质力学模型采用网格扩展的方式构建地质力学网格，x、y方向向外延伸油藏尺度的3倍以上；油藏上部延伸至地表；油藏下部考虑至避免薄板效应的厚度，即不少于横向尺度的1/3，具体参数见表3 ；整个地层设置为均匀弹性各向同性材料，Biot系数为1－Kp/Ks，设置油藏压力均匀下降10 MPa，变形边界条件为除地表外固定应变为零。
表3 VISAGE模型参数设置 |
| 算例 | 水平方向 延伸倍数 | 水平方向 延伸网格数 | 油藏上部 延伸网格数 | 油藏下部延伸所至深度/m | 油藏下部 延伸网格数 |
|---|---|---|---|---|---|
| 算例1 | 3 | 5 | 15 | 6 400 | 20 |
| 算例2 | 10 | 20 | 20 | 15 000 | 50 |
| 算例3 | 30 | 30 | 20 | 30 000 | 70 |
| 算例4 | 100 | 50 | 40 | 140 000 | 90 |
根据（27）式可得该问题的精确解，以垂向位移为例，其精确解如图6 所示。对于精确解，图6a 的四周边界存在一定的异常低值，图6c 的四周边界存在一定的异常高值，而图6b 的四周边界则不存在明显异常。这是由于侧向边界力对油藏中部垂向位移影响很小，案例中侧向边界力是指向油藏内的力，会使油藏上部产生向上的位移，油藏下部产生向下的位移，从而形成了图6a 与图6c 中四周边界与中心部分的差异，由于油藏厚度相对于平面尺度很小，因此其作用范围较小，主要影响了最外层求解点。
为评价不同方法的求解精度，定义相对误差为对应网格求解结果和精确解的差值与精确解之比的绝对值。不同方法求解结果如图6 所示，其相对误差如图7 所示。本文建立的体积边界元法求解结果与精确解基本一致，误差分布较为均匀，均在0.025%以内。体积边界元法求解精度远超商业模拟器VISAGE，使用VISAGE时，随着网格扩展范围的增加，油藏内部求解结果会逐步向精确解逼近，但在油藏角部边界处始终存在较大误差，对于油藏顶层和底层网格，其四周的边界网格基本存在10%以上的相对误差，这种程度的误差即使在工程上也难以接受。
传统有限元法在渗流边界处产生高误差的根本原因在于其理论框架无法精确反映渗流边界的存在。具体而言，该方法对渗流边界的处理是认为油藏渗流边界外无渗透性，即渗透率为零。初始条件下，孔隙压力从油藏渗流模拟网格映射到有限元网格（油藏区域），并外推到三维有限元网格的扩展区域。随后三维有限元网格的扩展区域由于不具有渗透性而保持压力恒定，通过油藏区域与不渗透区域间的压力梯度作用来反映渗流边界的作用。这种处理方法就导致边界作用与网格有关，因此在边界位置会存在较大误差。本文所建的体积边界元法则突破了这一理论局限，直接将边界等效力作用于确定的边界面，与网格属性无关，因此在边界处也保持与内部区域一致的数值精度。
不同方法在油藏网格点求解误差的详细对比如表4 所示。可见体积边界元法在油藏网格点上的最大误差仅为0.024%，而商业模拟器VISAGE在油藏网格点上的最大误差均在40%以上，且其最大误差并不会随网格扩展范围的增加而减小，反而是最小误差会随网格扩展范围的增加而减小，但这对于求解结果的可靠性并没有太大帮助，此外商业模拟器求解时间会随网格数量的增加而显著增加。本文建立的体积边界元法实现了远超商业模拟器的求解精度，与商业模拟器精细模型算例4相比，体积边界元法在油藏区域平均误差不到商业模拟器的1/500，而求解时间缩短79%。
表4 不同模拟方法求解油藏网格位移性能对比 |
| 方法 | 网格数 | 求解 时间/min | 最小 误差/% | 最大 误差/% | 平均 误差/% |
|---|---|---|---|---|---|
| 体积边界元法 | 31.0 | 2.5×10-14 | 0.024 | 0.008 2 | |
| VISAGE（算例1） | 40 960 | 0.6 | 0.650 00 | 67.210 | 10.770 0 |
| VISAGE（算例2） | 288 300 | 7.0 | 0.100 00 | 47.380 | 5.150 0 |
| VISAGE（算例3） | 638 780 | 21.5 | 0.027 00 | 45.360 | 4.050 0 |
| VISAGE（算例4） | 2 009 340 | 149.2 | 0.000 33 | 51.530 | 4.150 0 |
本文建立的体积边界元法除了近场求解精度高外，还具有远近场求解精度相同的独特优势。地下油藏开发引起的地表沉降问题也是地质力学研究的重要对象，因此通过地表沉降求解结果来对比不同方法的远场求解精度。由于不同模型油藏外网格大小不同，因此采用油藏在地表的垂直投影区域进行对比分析。不同方法在地表处求解误差的详细对比如表5 所示。可见体积边界元法远场求解精度优势更加显著，几乎可以认为误差为零，而商业模拟器在地表处的平均误差在1%~4%，且远场求解误差并不会随网格扩展范围的增加而明显减小，与商业模拟器精细模型算例4相比，体积边界元法远场求解平均误差不到商业模拟器的1/1×1012。
表5 不同模拟方法求解地表沉降性能对比 |
| 方法 | 顶层网格 厚度/m | 最小 误差/% | 最大 误差/% | 平均 误差/% |
|---|---|---|---|---|
| 体积边界元法 | 2.15×10-13 | 1.95×10-12 | 1.11×10-12 | |
| VISAGE（算例1） | 668.2 | 0.10 | 2.33 | 1.21 |
| VISAGE（算例2） | 464.0 | 1.83 | 3.25 | 2.72 |
| VISAGE（算例3） | 464.0 | 3.00 | 4.08 | 3.68 |
| VISAGE（算例4） | 216.6 | 1.76 | 2.22 | 2.05 |
基于方形油藏均匀压降模型的对比研究表明，本文提出的体积边界元法在求解精度方面显著优于传统有限元数值求解器，无论在油藏网格点还是远场区域，该方法均展现出更优的求解精度和稳定性。值得注意的是，VISAGE等商业模拟器在计算中普遍忽略孔隙流体密度与饱和度动态变化对地层变形的影响，因此上述验证中精确解及体积边界元模型均未考虑流体密度及饱和度的变化。若以考虑流体密度及饱和度变化的精确解为基准，传统商业模拟器在油藏网格点和远场的平均误差将增加约20个百分点（如表6 所示），而体积边界元法可以完整考虑这些影响，始终保持极低的误差水平。此外，表6 中商业模拟器在油藏网格点的求解误差随网格扩展范围的增加反而增大，这是因为随网格扩展范围的增加，商业模拟器的求解结果逐步逼近忽略流体密度及饱和度变化的解，反而导致其与实际情况的偏差进一步增大。
表6 考虑流体密度及饱和度变化时不同方法求解误差对比 |
| 位置 | 方法 | 最小 误差/% | 最大 误差/% | 平均 误差/% |
|---|---|---|---|---|
| 油藏 | 体积边界元法 | 0.25×10-4 | 0.023 | 0.008 |
| VISAGE（算例1） | 0.39 | 83.400 | 25.680 | |
| VISAGE（算例2） | 2.07 | 103.860 | 29.390 | |
| VISAGE（算例3） | 0.43 | 109.230 | 30.650 | |
| VISAGE（算例4） | 1.58 | 112.240 | 31.890 | |
| 地面 | 体积边界元法 | 4.54×10-13 | 2.50×10-12 | 1.67×10-12 |
| VISAGE（算例1） | 18.24 | 18.51 | 18.34 | |
| VISAGE（算例2） | 19.57 | 21.33 | 20.15 | |
| VISAGE（算例3） | 20.54 | 22.73 | 21.27 | |
| VISAGE（算例4） | 18.38 | 21.25 | 19.37 |
4 结论
油藏内部渗流与渗流边界对地层变形的影响分别等效于油藏内部和油藏边界处受到外力作用时对地层变形产生的影响，本文给出了渗流等效力和边界等效力的计算方法，渗流等效力由孔隙压力梯度、流体质量变化及岩石弹性参数决定，边界等效力由孔隙压力变化量和岩石弹性参数决定。地层任意点的变形解可通过渗流等效力、边界等效力与格林函数的卷积得到。离散化后，地层任意点的变形解可通过将网格边界等效力与对应的网格边界源相乘、网格渗流等效力与对应的网格体积源相乘再累加得到，即体积边界元数值求解方法。
体积边界元法具有无需对油藏外地层进行网格划分、求解精度高等优点，完整考虑了油藏渗流边界、油藏内孔隙压力梯度场、孔隙内流体质量变化对地层变形的影响。与商业模拟器精细模型相比，体积边界元法求解精度显著提高，在地表等远场位置，精度优势更加显著。
本文推导过程采用了均质地层假设，尚未考虑地层非均质性及裂缝的影响，但本文模型的可拓展性较强，为后续研究提供了新方向。未来可通过更换不同假设下的修正Mindlin解、建立裂缝边界等效力，将本文模型推广应用于层状非均质地层、裂缝性地层等复杂地层。
符号注释：
a，b，c——x，y，z方向的网格单元序号，网格单元序号增加方向与坐标轴正方向一致；A——微元的横截面积，m2；Ct——综合压缩系数，Pa-1；D——深度，m；Dt——油藏顶深，m；fi——单位表观体积岩石骨架在i方向上所受的总外力，Pa/m；fe,i——i方向的渗流等效力，Pa/m；Fi——i方向流体作用在骨架上的体积力，Pa/m；${{\tilde{F}}_{j\_i}}$——网格边界源，m3/N；grid(a, b, c)——x，y，z方向序号分别为a，b，c的网格；G——岩石表观剪切模量，Pa；g——重力加速度，m/s2；h——油藏厚度，m；Kp——岩石表观体积模量，Pa；Ks——骨架材料体积模量，Pa；L——油藏边长，m；${{\mathbf{n}}_{\Omega }}$， ——封闭边界内部空间和外部空间在封闭边界Γ处的单位内法向向量；p——孔隙流体压力，Pa；${{p}_{0}}$——原始地层压力，Pa；pe,i——i方向的边界等效力，Pa；R——油藏半径，m；ui——岩石在i方向上的位移，m；${{\tilde{u}}_{i\_j}}$——位移格林函数，m/N；${{\tilde{V}}_{j\_i}}$——网格体积源，m4/N；x，y，z——笛卡尔坐标系，m；xi——笛卡尔坐标系中i方向的坐标，m；Δxi——微元长度，m；Δx，Δy，Δz——x，y，z方向网格单元的大小，m；α，β，γ——地层任意点坐标，m；Γ——封闭边界；δij——克罗内克符号，若i≠j，δij=0，若i=j，δij=1；ε——岩石表观体积应变，无因次；${{\varepsilon }_{ij}}$——岩石表观应变，无因次；$\nu$——岩石表观泊松比，无因次；${{\rho }_{f}}$——岩石孔隙内流体的平均密度，kg/m3；${{\rho }_{s}}$——岩石骨架材料的密度，kg/m3；${{\rho }_{ref}}$——参考压力下流体的平均密度，kg/m3；$\bar{\rho }$——饱和流体岩石平均密度，kg/m3；σij——岩石骨架应力（拉为正，压为负），Pa；σs,ij——真实骨架应力，Pa；σt,ij——饱和岩石总应力，Pa；ϕ——岩石孔隙度，%；Ω——封闭边界内部空间； ——封闭边界外部空间。