CN104392131A - 一种水驱沙过程中破碎岩石渗流场计算方法 - Google Patents

Abstract

本发明涉及一种水驱沙过程中破碎岩石渗流场计算方法，根据水、沙渗透过程中的质量守恒方程、动量守恒方程、孔隙压缩方程以及渗透性参量-沙体积分数关系式，建立了渗流过程中物理量关系框图；根据物理量关系框图设计渗流场计算方法，该方法考虑了水驱沙过程中，由于沙质量流失水引起的沙体积分数的变化；考虑了沙体积分数对水和沙的流度、非Darcy流β因子、加速度系数的影响；可计算出各节点渗透性参量；在算法设计中，水和沙的流度、非Darcy流β因子、加速度系数，随孔隙度和沙体积分数变化，时变渗流场利用异步迭代的方法计算，节省储存空间，运算速度快，收敛性和数值稳定性好，计算结果精确。

Description

一种水驱沙过程中破碎岩石渗流场计算方法

技术领域

本发明涉及煤矿地质灾害与防治领域，具体涉及一种水驱沙过程中破碎岩石渗流场计算方法。

背景技术

在毛乌素沙漠，煤层具有埋深小、地表沉沙厚的特点。在采煤过程中，地表的砂层随地表水沿着煤层覆岩向采空区运移。为了研究煤层开采过程中由于覆岩破坏造成突水溃沙的机理、防治突水溃沙的发生或减轻突水溃沙的危害，根据水驱沙试验得到的水相和沙相渗透率、非Darcy流β因子、加速度系数随破碎岩石孔隙中沙子体积分数变化规律，建立一种水驱沙过程中破碎岩石渗流场计算方法。

目前，人们描述和模拟水驱沙过程中破碎岩石孔隙中水、沙两相介质的流动主要参照油藏工程中水驱开发理论和方法。由于沙子本质上属于非连续介质，沙子在破碎岩石孔隙中的迁移与油在油藏中渗流的运动机理、运动方式和阻力特性等方面存在显著的差异。特别是，水和沙之间不存在分界面。因此，参照水驱油方法计算水驱沙过程中破碎岩石的渗流场这种做法欠妥，构建一种有别于水驱油渗流场计算方法的、水沙在破碎岩石中水沙两相渗流计算方法既是必要的。

发明内容

本发明的目的是提供一种水驱沙过程中破碎岩石渗流场计算方法，根据水驱沙试验得到的水相和沙相渗透率、非Darcy流β因子、加速度系数随破碎岩石孔隙中沙子体积分数变化关系，建立水驱沙过程中破碎岩体中水沙渗流场的计算方法

本发明是通过以下技术方案实现的：一种水驱沙过程中破碎岩石渗流场计算方法，在水驱沙过程中破碎岩石孔隙被水和沙两种介质占据，根据水、沙渗透过程中的质量守恒方程、动量守恒方程、孔隙压缩方程以及渗透性参量-沙体积分数关系式，建立渗流过程中物理量关系框图；根据物理量关系框图设计渗流场计算方法。

本发明的基本原理为：

分别用下标1和2表示水相和沙相介质的物理量，破碎岩石孔隙中质量守恒方程为

$m_1{c}_{a1}\frac{\∂V_1}{\∂t}=-\frac{\∂p}{\∂x}-\frac{1}{I_1}{V}_1-m_1{\mathrm{\β}}_1{V_1}^2---\left(1\right)$

和

$m_2{c}_{a2}\frac{\∂V_2}{\∂t}=-\frac{\∂p}{\∂x}-\frac{1}{I_{2e}}{V}_2^n-m_2{\mathrm{\β}}_2{V}_2^2---\left(2\right)$

其中，t为时间，x为空间坐标，p为水压，V渗流速度，m为质量密度，I为流度，I_e为有效流度，β为非Darcy流β因子，c_a为加速度系数，n为粘度指数；

b.假设在破碎岩石孔隙中没有空气，即孔隙体积完全被水和沙占据，记水和沙的体积分数分别Y₁和Y₂，则有如下关系

Y₁+Y₂＝1                 (3)

破碎岩石孔隙中水、沙的质量守恒方程分别为

$\frac{\∂\left(\mathrm{\φ}{Y}_1\right)}{\∂t}+\frac{\∂V_1}{\∂x}=0---\left(4\right)$

和

$\frac{\∂\left(\mathrm{\φ}{Y}_2\right)}{\∂t}+\frac{\∂V_2}{\∂x}=0---\left(5\right)$

将式(4)和式(5)相加，并利用式(3)，可以得到

$\frac{\∂\mathrm{\φ}}{\∂t}=-\frac{\∂}{\∂x}(V_1+V_2)---\left(6\right)$

c.破碎岩石的孔隙具有压缩性，即孔隙度随压力变化，亦即孔隙压缩方程为

$\mathrm{\φ}={\mathrm{\φ}}_{p_0}[1+c_{\mathrm{\φ}}(p-p_0)]---\left(7\right)$

由式(4)、式(5)和式(7)，可以得到

$\frac{\∂p}{\∂t}=-\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}\frac{\∂}{\∂x}(V_1+V_2)---\left(8\right)$

式(1)～(3)、式(5)～(6)和式(8)构成六个变量p，φ，V₁，V₂，Y₁，Y₂的封闭方程组，此方程组中水的质量密度m₁、沙的质量密度m₂、破碎岩石的孔隙压缩系数c_φ和参考压力p₀下的孔隙度可近似当做常量来处理，而水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2则随孔隙度φ和体积分数Y₂变化。

如何利用试验数据建立水相渗透性参量(流度I₁、非Darcy流β因子β₁和加速度系数c_a1)和沙相渗透性参量(有效流度I_2e、非Darcy流β因子β₂和加速度系数c_a2)是本专利的关键所在。

本发明的具体计算方法为：

步骤1)计算得到的若干级孔隙度下的渗透性参量-沙体积分数关系式，利用一元三点不等距Lagrange插住法分别计算出各节点在当前孔隙度下水相流度I₁、非Darcy流β因子β₁、加速度系数c_a1、沙相有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2、渗透性参量-沙体积分数Y₂关系式系数；

步骤1.1)计算渗透性参量-沙体积分数关系

对于某一特定孔隙度的破碎岩石，通过水驱沙试验和线性回归，得到水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2与沙子体积分数Y₂的关系，即

$I_1=A_{10}+A_{11}{Y}_2+A_{12}{Y}_2^2---\left(9a\right)$

$\lg {\mathrm{\β}}_1=A_{20}+A_{21}{Y}_2+A_{22}{Y}_2^2---\left(9b\right)$

$\lg c_{a1}=A_{30}+A_{31}{Y}_2+A_{32}{Y}_2^2---\left(9c\right)$

$I_{2e}=A_{40}+A_{41}{Y}_2+A_{42}{Y}_2^2---\left(9d\right)$

$\lg {\mathrm{\β}}_2=A_{50}+A_{51}{Y}_2+A_{52}{Y}_2^2---\left(9e\right)$

$\lg c_{a2}=A_{60}+A_{61}{Y}_2+A_{62}{Y}_2^2---\left(9f\right)$

将回归系数A_ij(i＝1，2，…，6；j＝0，1，2)组装成矩阵

$\underset{\‾}{A}=\left[\begin{array}{ccc}{A}_{10}& A_{11}& A_{12}\\ A_{20}& A_{21}& A_{22}\\ A_{30}& A_{31}& A_{32}\\ A_{40}& A_{41}& A_{42}\\ A_{50}& A_{51}& A_{52}\\ A_{60}& A_{61}& A_{62}\end{array}\right]---\left(10\right)$

步骤1.2)计算渗透性参量-体积分数回归系数三维数组的构造

对同一粒径破碎岩石设置m_p个不同的孔隙度 通过水驱沙试验得到m_p个不同孔隙度下的回归系数矩阵A ¹，A ²，A ³，…，

将A ¹，A ²，A ³，…，组装成三维矩阵A ^b，其中A ^b的第i层为A ⁱ，为了读数方便，也可将A ¹，A ²，A ³，…，组装成二维分块矩阵A ^p，即

${\underset{\‾}{A}}^p=\left[\begin{array}{c}{\underset{\‾}{A}}^1\\ {\underset{\‾}{A}}^2\\ ...\\ {\underset{\‾}{A}}^{m_p}\end{array}\right]---\left(11\right)$

或

${\underset{\‾}{A}}_{(i+6k-6)j}^p={\underset{\‾}{A}}_{\mathrm{kij}}^b,k=\mathrm{1,2},...,m_p;i=\mathrm{1,2},...6;j=\mathrm{1,2,3}---\left(12\right)$

步骤1.3)渗透性参量的计算

在水驱沙过程中，孔隙度随时间变化，当孔隙度介于φ₁和之间时，应用一元三点Lagrange插值的方法计算系数A_ij(i＝1，2，…，6；j＝0，1，2)，即

$\begin{array}{c}{A}_{\mathrm{ij}}=A_{\mathrm{ij}}^k\frac{(\mathrm{\φ}-{\mathrm{\φ}}_{k+1})(\mathrm{\φ}-{\mathrm{\φ}}_{k+2})}{(\mathrm{\φ}-{\mathrm{\φ}}_{k+1})({\mathrm{\φ}}_k-{\mathrm{\φ}}_{k+2})}+A_{\mathrm{ij}}^{k+1}\frac{(\mathrm{\φ}-{\mathrm{\φ}}_k)(\mathrm{\φ}-{\mathrm{\φ}}_{k+2})}{({\mathrm{\φ}}_{k+1}-{\mathrm{\φ}}_k)({\mathrm{\φ}}_{k+1}-{\mathrm{\φ}}_{k+2})}\\ +A_{\mathrm{ij}}^{k+2}\frac{(\mathrm{\φ}-{\mathrm{\φ}}_k)(\mathrm{\φ}-{\mathrm{\φ}}_{k+1})}{({\mathrm{\φ}}_{k+2}-{\mathrm{\φ}}_k)({\mathrm{\φ}}_{k+2}-{\mathrm{\φ}}_{k+1})},\mathrm{\φ}\∈({\mathrm{\φ}}_k,{\mathrm{\φ}}_{k+2});i=\mathrm{1,2},...,6;j=\mathrm{0,1,2}\end{array}---\left(13\right)$

再根据式(9)计算水相和沙相的渗透性参量。

步骤2)根据破碎岩石的质量、高度、半径和真密度，计算出初始时刻(t＝0)破碎试验的孔隙度，根据沙子的质量、真密度计算出初始时刻沙的体积分数，并计算出水的体积分数；设定渗透条件，写出节点压力和渗流速度的分布；

步骤2.1)设定边界条件和初始条件

根据试验条件，岩样上端x＝0的压力保持恒定，即

p|_x＝0＝p₀    (14)

下端与大气连通，故

p|_x＝h＝p₁＝0    (15)

岩样上端x＝0没有沙子流入，故渗流速度

V₂|_x＝0＝0    (16)

步骤2.2)根据破碎岩石的质量M_rock、高度h、半径a_rock和破碎前的真密度m_rock，可以计算出初始时刻(t＝0)破碎试验的孔隙度

${\mathrm{\φ}|}_{t=0}={\mathrm{\φ}}^{\mathrm{initial}}=\frac{\mathrm{\π}{a}_{\mathrm{rock}}^{2}h-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}{\mathrm{\π}{a}_{\mathrm{rock}}^{2}h}---\left(17\right)$

步骤2.3)根据沙子的质量M_sand和真密度m_sand，可以计算出初始时刻沙的体积分数，即

${Y_2|}_{t=0}=Y_2^{\mathrm{initail}}=\frac{\frac{M_{\mathrm{sand}}}{m_{\mathrm{sand}}}}{\mathrm{\π}{a}_{\mathrm{rock}}^{2}h-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}---\left(18\right)$

步骤2.4)利用式(3)可以计算出初始时刻水的体积分数，即

${Y_1|}_{t=0}=1-{Y_2|}_{t=0}=1-Y_2^{\mathrm{initail}}---\left(19\right)$

步骤2.5)初始压力呈分布

${p|}_{t=0}=p_0+\frac{x}{H}(p_1-p_0)---\left(20\right)$

步骤2.6)初始时刻水渗流速度V₁

${V_1|}_{t=0}=\frac{p_0-p_1}{h}{I}_1---\left(21\right)$

步骤2.7)沙渗流速度为

${V_2|}_{t=0}={\left(\frac{p_0-p_1}{h}{I}_1\right)}^{\frac{1}{n}}---\left(22\right)$

步骤3)根据节点物理量计算单元上孔隙度对坐标的偏导数、压力对坐标的偏导数、体积分数对坐标的偏导数、渗流速度对坐标的偏导数；通过节点附近区域平衡方程的积分，分别计算出压力、孔隙度、体积分数和渗流速度对时间的偏导数，在通过时间步长上的积分，计算出下一时刻节点的压力、孔隙度、体积分数和渗流速度；

步骤3.1)水驱沙过程中水沙两相流动的一维流场算法的核心是将物理量分为两类：单元物理量和节点物理量，将区间[0，h]均分为N个单元，单元长度为节点k坐标为

$x_n\left(k\right)=(k-1)h_e=\frac{k-1}{N}h,k=\mathrm{1,2},...,N+1;---\left(23\right)$

单元k的左、右端点坐标分别为x_k和

节点物理量包括孔隙度φ_n(k)、体积分数Y_1n(k)、Y_2n(k)、渗流速度V_1n(k)、V_2n(k)、压力p_n(k)，k＝1，2，…，N+1，单元物理量包括孔隙度梯度体积分数梯度渗流速度梯度和压力梯度k＝1，2，…，N；

单元k上物理量对坐标偏导数的差分格式分别为

${\left(\frac{\∂\mathrm{\φ}}{\∂x}\right)}_k=\frac{{\mathrm{\φ}}_{k+1}-{\mathrm{\φ}}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(24\right)$

${\left(\frac{\∂Y_1}{\∂x}\right)}_k=\frac{{\left(Y_1\right)}_{k+1}-{\left(Y_1\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(25\right)$

${\left(\frac{\∂Y_2}{\∂x}\right)}_k=\frac{{\left(Y_2\right)}_{k+1}-{\left(Y_2\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(26\right)$

${\left(\frac{\∂V_1}{\∂x}\right)}_k=\frac{{\left(Y_1\right)}_{k+1}-{\left(Y_1\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(27\right)$

${\left(\frac{\∂V_2}{\∂x}\right)}_k=\frac{{\left(Y_2\right)}_{k+1}-{\left(Y_2\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(28\right)$

${\left(\frac{\∂p}{\∂x}\right)}_k=\frac{p_{k+1}-p_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(29\right)$

步骤3.2)基于式(8)构造时刻t^j节点i上压力对时间的偏导数，即在区间上对式(8)积分，得到

${\left(\frac{\∂p}{\∂t}\right)}_i^{t_j}\frac{x_{i+1}-x_{i-1}}{2}=-\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_{i-1}^{t_j}\frac{x_i-x_i}{2}-\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_i^{t_j}\frac{x_{i+1}-x_i}{2}---\left(30\right)$

对于等距节点，式(21)可以简化为

${\left(\frac{\∂p}{\∂t}\right)}_i^{t_j}=-\frac{1}{2}[\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_{i-1}^{t_j}+\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_i^{t_j}]---\left(31\right)$

将式(22)对时间积分一步，得到

$p_i^{t_{j+1}}=p_i^{t_j}-\frac{1}{2}[\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_{i-1}^{t_j}+\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{\∂V_1}{\∂x}+\frac{\∂V_2}{\∂x})}_i^{t_j}]h_t---\left(32\right)$

步骤3.3)基于式(7)构造节点孔隙度的算法，即

${\mathrm{\φ}}_i^{t_{j+1}}={\mathrm{\φ}}_{p_0}[1+c_{\mathrm{\φ}}(p_i^{t_{j+1}}-p_0)].---\left(33\right)$

步骤3.4)基于式(5)和式(3)构造体积分数Y₂和Y₁的算法，即

$\begin{array}{c}{\left(Y_2\right)}_i^{t_{j+1}}={\left(Y_2\right)}_i^{t_j}+\frac{h_t}{2{\left({\mathrm{\φ}}_e\right)}_{i-1}^{t_j}}[{\left(Y_{2e}\right)}_{i-1}^{t_j}{\left(\frac{\∂V_1}{\∂x}\right)}_{i-1}^{t_j}-{\left(Y_{1e}\right)}_{i-1}^{t_j}{\left(\frac{\∂V_2}{\∂x}\right)}_{i-1}^{t_j}]\\ +\frac{h_t}{2{\left({\mathrm{\φ}}_e\right)}_i^{t_j}}[{\left(Y_{2e}\right)}_{i-1}^{t_j}{\left(\frac{\∂V_1}{\∂x}\right)}_i^{t_j}-{\left(Y_{1e}\right)}_{i-1}^{t_j}{\left(\frac{\∂V_2}{\∂x}\right)}_i^{t_j}]\end{array}---\left(34\right)$

${\left(Y_1\right)}_i^{t_{j+1}}=1-{\left(Y_2\right)}_i^{t_{j+1}}---\left(35\right)$

步骤3.5)基于式(1)和式(2)构造渗流速度V₁和V₂的算法，即

$\begin{array}{c}{\left(V_1\right)}_i^{t_{j+1}}={\left(V_1\right)}_i^{t_j}-\\ [\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1{}_{}}^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_{i-1}^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1}^{t_j}{\left(I_1\right)}_{i-1}^{t_j}}{\left(V_1\right)}_{i-1}^{t_j}+\frac{{\left({\mathrm{\β}}_1\right)}_{i-1}^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1}^{t_j}}{\left({V_1}^2\right)}_{i-1}^{t_j}]\frac{x_i-x_i}{x_{i+1}-x_{i-1}}{h}_t\\ -[\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_i^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}{\left(I_1\right)}_i^{t_j}}{\left(V_1\right)}_i^{t_j}+\frac{{\left({\mathrm{\β}}_1\right)}_i^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}}{\left({V_1}^2\right)}_i^{t_j}]\frac{x_{i+1}-x_i}{x_{i+1}-x_{i-1}}{h}_t\end{array}---\left(36\right)$

$\begin{array}{c}{\left(V_2\right)}_i^{t_{j+1}}={\left(V_2\right)}_i^{t_j}\\ -[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1{}_{}}^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_{i-1}^{t_j}+\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1}^{t_j}{\left(I_{2\mathrm{ee}}\right)}_{i-1}^{t_j}}{\left(V_{2e}^n\right)}_{i-1}^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_{i-1}^{t_j}}{{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1}^{t_j}}{\left(V_{2e}^2\right)}_{i-1}^{t_j}]\frac{x_i-x_i}{x_{i+1}-x_{i-1}}{h}_t\\ -[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_i^{t_j}+\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}{\left(I_{2\mathrm{ee}}\right)}_i^{t_j}}{\left(V_{2e}\right)}_i^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_i^{t_j}}{{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}}{\left(V_{2e}^2\right)}_i^{t_j}]\frac{x_{i+1}-x_i}{x_{i+1}-x_{i-1}}{h}_t\end{array}---\left(37\right)$

步骤4)根据节点孔隙度和沙体积分数，利用一元三点不等距Lagrange插住法分别计算出各节点处水相流度、非Darcy流β因子、加速度系数、沙相有效流度、非Darcy流β因子、加速度系数-沙体积分数关系式系数，并利用渗透性参量-沙体积分数关系式计算节点上渗透性参量的值，进入下一步循环。

步骤4.1)单元划分

沿破碎岩样高度方向建立坐标轴Ox，坐标轴的正向与渗流速度方向相同，将破碎岩样均分为N个单元，单元长度为节点k坐标为

$x_n\left(k\right)=(k-1)h_e=\frac{k-1}{N}h,k=\mathrm{1,2},...,N+1.---\left(38\right)$

单元k的左、右端点坐标分别为x_k和x_k+1。

步骤4.2)渗透参量-体积分数回归系数的确定

步骤4.2.1)根据试验结果，给出m_p级孔隙度下的渗透参量-体积分数回归系数，孔隙度分别记为φ¹、φ²、……、对于的渗透参量-体积分数回归系数构成的矩阵分别记为

$\underset{\‾}{A^1}=\left[\begin{array}{ccc}{A}_{10}^1& A_{11}^1& A_{12}^1\\ A_{20}^1& A_{21}^1& A_{22}^1\\ A_{30}^1& A_{31}^1& A_{32}^1\\ A_{40}^1& A_{41}^1& A_{42}^1\\ A_{50}^1& A_{51}^1& A_{52}^1\\ A_{60}^1& A_{61}^1& A_{62}^1\end{array}\right],\underset{\‾}{A^1}=\left[\begin{array}{ccc}{A}_{10}^2& A_{11}^2& A_{12}^2\\ A_{20}^2& A_{21}^2& A_{22}^2\\ A_{30}^2& A_{31}^2& A_{32}^2\\ A_{40}^2& A_{41}^2& A_{42}^2\\ A_{50}^2& A_{51}^2& A_{52}^2\\ A_{60}^2& A_{61}^2& A_{62}^2\end{array}\right],......,\underset{\‾}{A^{m_p}}=\left[\begin{array}{ccc}{A}_{10}^{m_p}& A_{11}^{m_p}& A_{12}^{m_p}\\ A_{20}^{m_p}& A_{21}^{m_p}& A_{22}^{m_p}\\ A_{30}^{m_p}& A_{31}^{m_p}& A_{32}^{m_p}\\ A_{40}^{m_p}& A_{41}^{m_p}& A_{42}^{m_p}\\ A_{50}^{m_p}& A_{51}^{m_p}& A_{52}^{m_p}\\ A_{60}^{m_p}& A_{61}^{m_p}& A_{62}^{m_p}\end{array}\right]$

步骤4.2.2)将A ¹ 、A ² 、……、组装为分块矩阵

${\underset{\‾}{A}}^p=\left[\begin{array}{c}{\underset{\‾}{A}}^1\\ {\underset{\‾}{A}}^2\\ ...\\ {\underset{\‾}{A}}^{m_p}\end{array}\right]---\left(11\right)$

步骤4.3)给出节点物理量的初始值

步骤4.3.1)孔隙度的初始值

${{\mathrm{\φ}}_n\left(k\right)|}_{t=0}={\mathrm{\φ}}_{\mathrm{initial}}=\frac{{\mathrm{\π\α}}_{\mathrm{rock}}^{2}h-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}}{{\mathrm{\π\α}}_{\mathrm{rock}}^{2}h},k=\mathrm{1,2},....N+1---\left(39\right)$

步骤4.3.2)沙体积分数的初始值

$Y_{2n}\left(k\right)=\frac{\frac{M_{\mathrm{sand}}}{m_{\mathrm{sand}}}}{{\mathrm{\π\α}}_{\mathrm{rock}}^{2}h-\frac{M_{\mathrm{rock}}}{m_{\mathrm{rock}}}},k=\mathrm{1,2},....N+1---\left(40\right)$

步骤4.3.3)水体积分数的初始值

Y_1n(k)＝1-Y_1n(k)，k＝1，2，…，N+1   (41)

步骤4.3.4)压力的初始值

$p_n\left(k\right)=p_0+\frac{x_k}{H}(p_1-p_0),k=\mathrm{1,2},...,N+1---\left(42\right)$

步骤4.3.5)水渗流速度的初始值

$V_{1n}\left(k\right)=\frac{p_0-p_1}{h}{I}_1,k=\mathrm{1,2},...,N+1---\left(43\right)$

步骤4.3.6)沙渗流速度的初始值

$V_{2n}\left(k\right)={\left(\frac{p_0-p_1}{h}{I}_1\right)}^{\frac{1}{n}},k=\mathrm{1,2},...,N+1---\left(44\right)$

步骤4.4初始时刻的渗透性参量计算

步骤4.4.1)利用试验结果得到的m_p级孔隙度φ⁽ⁱ⁾(i＝1，2，…，n_f)下的渗透性参量系数

$a_{\mathrm{jk}}^{\left(i\right)}(i=\mathrm{1,2},...,m_p;j=\mathrm{1,2},...,6;k=\mathrm{0,1,2})$

通过一元三点插值法计算初始孔隙度φ_initial下渗透性参量系数

a_jk(j＝1，2，…，6；k＝1，2，3)；

步骤4.4.2)根据式(9)计算节点上流度、非Darcy流β因子和加速度系数；

步骤4.5)单元物理量计算

步骤4.5.1)单元k上孔隙度对坐标的偏导数按式(24)计算，即

${\left(\frac{\∂\mathrm{\φ}}{\∂x}\right)}_k=\frac{{\mathrm{\φ}}_{k+1}-{\mathrm{\φ}}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(24\right)$

步骤4.5.2)体积分数对坐标的偏导数和按式(25)和式(26)计算，即

${\left(\frac{{\∂Y}_1}{\∂x}\right)}_k=\frac{{\left(Y_1\right)}_{k+1}-{\left(Y_1\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(25\right)$

${\left(\frac{{\∂Y}_2}{\∂x}\right)}_k=\frac{{\left(Y_2\right)}_{k+1}-{\left(Y_2\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(26\right)$

步骤4.5.3)渗流速度对坐标的偏导数和按式(27)和式(28)计算，即

${\left(\frac{{\∂V}_1}{\∂x}\right)}_k=\frac{{\left(Y_1\right)}_{k+1}-{\left(Y_1\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(27\right)$

${\left(\frac{{\∂V}_2}{\∂x}\right)}_k=\frac{{\left(Y_2\right)}_{k+1}-{\left(Y_2\right)}_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(28\right)$

步骤4.5.4)压力对坐标的偏导数按式(29)计算，即

${\left(\frac{\∂p}{\∂x}\right)}_k=\frac{p_{k+1}-p_k}{x_{k+1}-x_k},k=\mathrm{1,2},...,N---\left(29\right)$

步骤4.6)节点压力的计算

步骤4.6.1)内部节点的压力按式(32)计算，即

$p_i^{t_{j+1}}=p_i^{t_j}-\frac{1}{2}[\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{{\∂V}_1}{\∂x}+\frac{{\∂V}_2}{\∂x})}_{i-1}^{t_j}+\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{{\∂V}_1}{\∂x}+\frac{{\∂V}_2}{\∂x})}_i^{t_j}]h_t,i=\mathrm{2,3}...,N---\left(32\right)$

步骤4.6.2)边界节点1的压力按式(45)计算，即

$p_1^{t_{j+1}}=p_1^{t_j}-\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{{\∂V}_1}{\∂x}+\frac{{\∂V}_2}{\∂x})}_1^{t_j}{h}_t---\left(45\right)$

步骤4.6.3)边界节点N+1的压力按式(46)计算，即

$p_{N+1}^{t_{j+1}}=p_{N+1}^{t_j}-\frac{1}{{\mathrm{\φ}}_{p_0}{c}_{\mathrm{\φ}}}{(\frac{{\∂V}_1}{\∂x}+\frac{{\∂V}_2}{\∂x})}_N^{t_j}{h}_t---\left(46\right)$

步骤4.7)节点孔隙度的计算，节点孔隙度按式(33)计算，即

${\mathrm{\φ}}_i^{t_{j+1}}={\mathrm{\φ}}_{p_0}[1+c_{\mathrm{\φ}}(p_i^{t_{j+1}}-p_0)],i=\mathrm{1,2},...,N+1.---\left(33\right)$

步骤4.8)节点体积分数Y₂的计算

步骤4.8.1)内部节点的体积分数Y₂按式(34)计算，即

$\begin{array}{c}{\left(Y_2\right)}_i^{t_{j+1}}={\left(Y_2\right)}_i^{t_j}+\frac{h_t}{2{\left({\mathrm{\φ}}_e\right)}_{i-1}^{t_j}}[{\left(Y_{2e}\right)}_{i-1}^{t_j}{\left(\frac{{\∂V}_1}{\∂x}\right)}_{i-1}^{t_j}-{\left(Y_{1e}\right)}_{i-1}^{t_j}{\left(\frac{{\∂V}_2}{\∂x}\right)}_{i-1}^{t_j}]\\ +\frac{h_t}{2{\left({\mathrm{\φ}}_e\right)}_i^{t_j}}[{\left(Y_{2e}\right)}_{i-1}^{t_j}{\left(\frac{{\∂V}_1}{\∂x}\right)}_i^{t_j}-{\left(Y_{1e}\right)}_{i-1}^{t_j}{\left(\frac{{\∂V}_2}{\∂x}\right)}_i^{t_j}]\end{array}---\left(34\right)$

步骤4.8.2)边界节点1的体积分数Y₂按式(47)计算，即

${\left(Y_2\right)}_1^{t_j}={\left(Y_2\right)}_1^{t_j}+\frac{1}{{\left({\mathrm{\φ}}_e\right)}_1^{t_j}}[{\left(Y_{2e}\right)}_1^{t_j}{\left(\frac{{\∂V}_1}{\∂x}\right)}_1^{t_j}-{\left(Y_{1e}\right)}_1^{t_j}{\left(\frac{{\∂V}_2}{\∂x}\right)}_1^{t_j}]h_t---\left(47\right)$

步骤4.8.3)边界节点N+1的体积分数Y₂按式(48)计算，即

${\left(Y_2\right)}_{N+1}^{t_j}={\left(Y_2\right)}_N^{t_j}+\frac{1}{{\left({\mathrm{\φ}}_e\right)}_N^{t_j}}[{\left(Y_{2e}\right)}_1^{t_j}{\left(\frac{{\∂V}_1}{\∂x}\right)}_N^{t_j}-{\left(Y_{1e}\right)}_1^{t_j}{\left(\frac{{\∂V}_2}{\∂x}\right)}_N^{t_j}]h_t---\left(48\right)$

步骤4.9)节点体积分数Y₁的计算

节点体积分数Y₁按式(35)计算，即

${\left(Y_1\right)}_i^{t_{j+1}}=1-{\left(Y_2\right)}_i^{t_{j+1}},i=\mathrm{1,2},...,N+1---\left(35\right)$

步骤4.10)节点渗流速度V₁的计算

步骤4.10.1)内部节点的渗流速度V₁按式(36)计算，即

$\begin{array}{c}{\left(V_1\right)}_i^{t_{j+1}}={\left(V_1\right)}_i^{t_j}-\\ [\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1}^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_{i-1}^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1}^{t_j}{\left(I_1\right)}_{i-1}^{t_j}}{\left(V_1\right)}_{i-1}^{t_j}+\frac{{\left({\mathrm{\β}}_1\right)}_{i-1}^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_{i-1}^{t_j}}{\left({V_1}^2\right)}_{i-1}^{t_j}]\frac{x_i-x_{i-1}}{x_{i+1}-x_{i-1}}{h}_t\\ -[\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_i^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}{\left(I_1\right)}_i^{t_j}}{\left(V_1\right)}_i^{t_j}+\frac{{\left({\mathrm{\β}}_1\right)}_i^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_i^{t_j}}{\left({V_1}^2\right)}_i^{t_j}]\frac{x_{i+1}-x_{i-1}}{x_{i+1}-x_{i-1}}{h}_t,\\ i=\mathrm{2,3},...,N\end{array}---\left(36\right)$

步骤4.10.2)边界节点1渗流速度V₁按式(49)计算，即

$\begin{array}{c}{\left(V_1\right)}_1^{t_{j+1}}={\left(V_1\right)}_1^{t_j}\\ -[\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_1^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_1^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_1^{t_j}{\left(I_{1e}\right)}_1^{t_j}}{\left(V_{1e}\right)}_1^{t_j}+\frac{{\left({\mathrm{\β}}_{1e}\right)}_1^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_1^{t_j}}{\left(V_{1e}^2\right)}_1^{t_j}]h_t\end{array}---\left(49\right)$

步骤4.10.3)边界节点N+1渗流速度V₁按式(50)计算，即

$\begin{array}{c}{\left(V_1\right)}_{N+1}^{t_{j+1}}={\left(V_1\right)}_{N+1}^{t_j}\\ -[\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_N^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_1^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_N^{t_j}{\left(I_{1e}\right)}_N^{t_j}}{\left(V_{1e}\right)}_N^{t_j}+\frac{{\left({\mathrm{\β}}_{1e}\right)}_N^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_N^{t_j}}{\left(V_{1e}^2\right)}_N^{t_j}]h_t\end{array}---\left(50\right)$

步骤4.11)节点渗流速度V₂的计算

步骤4.11.1)内部节点的渗流速度V₂按式(37)计算，即

$\begin{array}{c}{\left(V_2\right)}_i^{t_{j+1}}={\left(V_2\right)}_i^{t_j}\\ -[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1}^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_{i-1}^{t_j}+\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1}^{t_j}{\left(I_{2\mathrm{ee}}\right)}_{i-1}^{t_j}}{\left(V_{2e}^n\right)}_{i-1}^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_{i-1}^{t_j}}{{\left({\stackrel{\‾}{c}}_a^2\right)}_{i-1}^{t_j}}{\left(V_{2e}^2\right)}_{i-1}^{t_j}]\frac{x_i-x_i}{x_{i+1}-x_{i-1}}{h}_t\\ -[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_i^{t_j}+\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}{\left(I_{2\mathrm{ee}}\right)}_i^{t_j}}{\left(V_{2e}\right)}_i^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_i^{t_j}}{{\left({\stackrel{\‾}{c}}_a^2\right)}_i^{t_j}}{\left(V_{2e}^2\right)}_i^{t_j}]\frac{x_{i+1}-x_i}{x_{i+1}-x_{i-1}}{h}_t,\\ i=\mathrm{2,3},...,N\end{array}---\left(37\right)$

步骤4.11.2)边界节点1渗流速度V₂按式(51)计算，即

${\left(V_2\right)}_1^{t_{j+1}}={\left(V_2\right)}_1^{t_j}-\begin{array}{c}[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_1^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_1^{t_j}+\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_1^{t_j}{\left(I_{2\mathrm{ee}}\right)}_1^{t_j}}{\left(V_{2e}\right)}_1^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_1^{t_j}}{{\left({\stackrel{\‾}{c}}_a^2\right)}_1^{t_j}}{\left(V_{2e}^2\right)}_1^{t_j}]h_t\end{array}---\left(51\right)$

步骤4.11.3)边界节点N+1渗流速度V₂按式(52)计算，即

$\begin{array}{c}{\left(V_2\right)}_{N+1}^{t_{j+1}}={\left(V_2\right)}_{N+1}^{t_j}\\ -[\frac{1}{m_2{\left({\stackrel{\‾}{c}}_a^2\right)}_N^{t_j}}{\left(\frac{\∂p}{\∂x}\right)}_1^{t_j}+\frac{1}{m_1{\left({\stackrel{\‾}{c}}_a^1\right)}_N^{t_j}{\left(I_{2\mathrm{ee}}\right)}_N^{t_j}}{\left(V_{2e}\right)}_N^{t_j}+\frac{{\left({\mathrm{\β}}_{2e}\right)}_N^{t_j}}{{\left({\stackrel{\‾}{c}}_a^1\right)}_N^{t_j}}{\left(V_{2e}^2\right)}_N^{t_j}]h_t\end{array}---\left(52\right)$

本发明的有益效果是：一种水驱沙过程中破碎岩石渗流场计算方法，(1)考虑了水驱沙过程中，由于沙质量流失水引起的沙体积分数Y₂的变化；(2)考虑了沙体积分数Y₂对水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2的影响；(3)根据试验得到m_p级孔隙度下水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2与沙子体积分数Y₂的回归系数利用一元三点不等距差值可以计算出各节点渗透性参量；(4)在算法设计中，水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2随孔隙度和沙体积分数Y₂变化，时变渗流场利用异步迭代的方法计算，节省储存空间，运算速度快，收敛性和数值稳定性好，计算结果精确。

附图说明

下面结合附图及实施例对本发明作进一步说明。

图1是变量间的关系框图；

图2是破碎岩石孔隙中水驱沙过程渗流场计算流程图；

图3是体积分数Y₂在不同时刻的分布；

图4是渗流速度V₂在不同时刻的分布；

图5是渗流速度V₁在不同时刻的分布；

图6是流度I₁在不同时刻的分布；

图7是非Darcy流β因子β₁在不同时刻的分布；

图8是加速度系数c_a1在不同时刻的分布；

图9是流度I_e2在不同时刻的分布；

图10是非Darcy流β因子β₂在不同时刻的分布；

图11是加速度系数c_a2在不同时刻的分布。

具体实施方式

变量间的关系框图是建立水驱沙过程中渗流场算法的基础。根据式(1)～(4)、式(5)、式(7)～(8)，构造出变量间的关系框图，如图1。

具体实施步骤如下：

1.试验数据处理

将某种粒径的破碎岩石与某种粒径的沙混合后，放入渗透仪进行压缩。根据式(17)，可以计算出破碎岩石孔隙度。设置m_p级孔隙度φ¹、φ²、……、分别进行每级孔隙度下破碎岩石孔隙中水驱沙试验，用式(9)拟合水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2与沙子体积分数Y₂之间的关系，将回归系数组装为分块矩阵

${\underset{\‾}{A}}^p=\left[\begin{array}{c}{\underset{\‾}{A}}^1\\ {\underset{\‾}{A}}^2\\ ...\\ {\underset{\‾}{A}}^{m_p}\end{array}\right]$

其中，

$\underset{\‾}{A^1}=\left[\begin{array}{ccc}{A}_{10}^1& A_{11}^1& A_{12}^1\\ A_{20}^1& A_{21}^1& A_{22}^1\\ A_{30}^1& A_{31}^1& A_{32}^1\\ A_{40}^1& A_{41}^1& A_{42}^1\\ A_{50}^1& A_{51}^1& A_{52}^1\\ A_{60}^1& A_{61}^1& A_{62}^1\end{array}\right],\underset{\‾}{A^2}\text{=}\left[\begin{array}{ccc}{A}_{10}^2& A_{11}^2& A_{12}^2\\ A_{20}^2& A_{21}^2& A_{22}^2\\ A_{30}^2& A_{31}^2& A_{32}^2\\ A_{40}^2& A_{41}^2& A_{42}^2\\ A_{50}^2& A_{51}^2& A_{52}^2\\ A_{60}^2& A_{61}^2& A_{62}^2\end{array}\right],......,\underset{\‾}{A^{m_p}}=\left[\begin{array}{ccc}{A}_{10}^{m_p}& A_{11}^{m_p}& A_{12}^{m_p}\\ A_{20}^{m_p}& A_{21}^{m_p}& A_{22}^{m_p}\\ A_{30}^{m_p}& A_{31}^{m_p}& A_{32}^{m_p}\\ A_{40}^{m_p}& A_{41}^{m_p}& A_{42}^{m_p}\\ A_{50}^{m_p}& A_{51}^{m_p}& A_{52}^{m_p}\\ A_{60}^{m_p}& A_{61}^{m_p}& A_{62}^{m_p}\end{array}\right]$

在本例中，孔隙度的级数m_p＝7，各级孔隙度的数值分别为φ¹＝0.562，φ²＝0.568，φ³＝0.574，φ⁴＝0.580，φ⁵＝0.586，φ⁶＝0.591、φ⁷＝0.597，相应的渗透参量-体积分数回归系数矩阵

$\underset{\‾}{A^1}=\left[\begin{array}{ccc}6.38\×{10}^{-13}& -2.23\×{10}^{-12}& 1.71\×{10}^{-12}\\ 4.77\×{10}^{+00}& 7.49\×{10}^{+00}& -5.85\×{10}^{+00}\\ 5.60\×{10}^{+00}& 5.41\×{10}^{+00}& -4.22\×{10}^{+00}\\ 7.48\×{10}^{-12}& -2.28\×{10}^{-11}& 1.58E\×{10}^{-11}\\ 8.19\×{10}^{+00}& 5.99\×{10}^{-01}& 3.88\×{10}^{+00}\\ 8.10\×{10}^{+00}& 4.12\×{10}^{-01}& 2.67\×{10}^{+00}\end{array}\right],$

$\underset{\‾}{A^2}=\left[\begin{array}{ccc}-1.42\×{10}^{-12}& 4.20\×{10}^{-11}& -2.39\×{10}^{-11}\\ 5.97\×{10}^{+00}& -1.14\×{10}^{+01}& 7.76\×{10}^{+00}\\ 6.46\×{10}^{+00}& -8.25\×{10}^{+00}& 5.60\×{10}^{+00}\\ 6.76\×{10}^{-12}& 5.02\×{10}^{-12}& -7.46\×{10}^{-12}\\ 8.23\×{10}^{+00}& -3.17\×{10}^{+00}& 4.27\×{10}^{+00}\\ 8.13\×{10}^{+00}& -2.18\×{10}^{+00}& 2.94\×{10}^{+00}\end{array}\right],$

$\underset{\‾}{A^3}=\left[\begin{array}{ccc}1.50\×{10}^{-12}& 4.25\×{10}^{-12}& -1.49\×{10}^{-12}\\ 4.06\×{10}^{+00}& -1.80\×{10}^{+00}& 9.84\×{10}^{-01}\\ 5.08\×{10}^{+00}& -1.30\×{10}^{+00}& 7.11\×{10}^{-01}\\ 1.48\×{10}^{-11}& 2.49\×{10}^{-11}& -4.17\×{10}^{-11}\\ 7.66\×{10}^{+00}& -5.10\×{10}^{+00}& 7.97\×{10}^{+00}\\ 7.74\×{10}^{+00}& -3.50\×{10}^{+00}& 5.48\×{10}^{+00}\end{array}\right]$

$\underset{\‾}{A^4}=\left[\begin{array}{ccc}4.93\×{10}^{-12}& 2.43\×{10}^{-11}& -1.95\×{10}^{-11}\\ 3.48\×{10}^{+00}& -4.18\×{10}^{+00}& 3.33\×{10}^{-00}\\ 4.67\×{10}^{+00}& -3.02\×{10}^{+00}& 2.40\×{10}^{+00}\\ 5.09\×{10}^{-11}& 1.32\×{10}^{-10}& -1.95\×{10}^{-10}\\ 6.99\×{10}^{+00}& -6.54\×{10}^{+00}& 9.50\×{10}^{+00}\\ 7.27\×{10}^{+00}& -4.50\×{10}^{+00}& 6.53\×{10}^{+00}\end{array}\right]$

$\underset{\‾}{A^5}=\left[\begin{array}{ccc}2.37\×{10}^{-11}& -1.28\×{10}^{-12}& 1.04\×{10}^{-12}\\ 1.85\×{10}^{+00}& 4.31\×{10}^{-02}& -3.48\×{10}^{-02}\\ 3.49\×{10}^{+00}& 3.11\×{10}^{-02}& -2.51\×{10}^{-02}\\ 2.61\×{10}^{-11}& 5.78\×{10}^{-10}& -8.92\×{10}^{-10}\\ 5.82\×{10}^{+00}& -6.07\×{10}^{+00}& 9.08\×{10}^{+00}\\ 6.47\×{10}^{+00}& -4.17\×{10}^{+00}& 6.25\×{10}^{+00}\end{array}\right]$

$\underset{\‾}{A^6}=\left[\begin{array}{ccc}2.41\×{10}^{-11}& -3.43\×{10}^{-14}& 1.58\×{10}^{-13}\\ 1.84\×{10}^{+00}& 1.10\×{10}^{-03}& -5.10\×{10}^{-03}\\ 3.48\×{10}^{+00}& 8.11\×{10}^{-04}& -3.70\×{10}^{-03}\\ 3.39\×{10}^{-10}& 7.52\×{10}^{-10}& -1.16\×{10}^{-09}\\ 5.65\×{10}^{+00}& -6.07\×{10}^{+00}& 9.08\×{10}^{+01}\\ 6.35\×{10}^{+00}& -4.17\×{10}^{+00}& 6.25\×{10}^{+00}\end{array}\right]$

$\underset{\‾}{A^7}=\left[\begin{array}{ccc}2.11\×{10}^{-11}& -1.25\×{10}^{-12}& 1.10\×{10}^{-12}\\ 1.94\×{10}^{+00}& 4.66\×{10}^{-02}& -4.10\×{10}^{-02}\\ 3.55\×{10}^{+00}& 3.36\×{10}^{-02}& -2.96\×{10}^{-02}\\ 3.85\×{10}^{-10}& 6.16\×{10}^{-10}& -1.07\×{10}^{-09}\\ 5.10\×{10}^{+00}& -6.44\×{10}^{+00}& 9.92\×{10}^{+00}\\ 5.84\×{10}^{+00}& -4.65\×{10}^{+00}& 7.17\×{10}^{+00}\end{array}\right]$

2.输入有关参数

输入有关参数，包括以下几类：

(1)破碎岩体粒径、沙粒径；

(2)水的质量密度m₁、动力粘度μ₁；

(3)沙的质量密度m₂、稠度系数C和幂指数n；

(4)初始孔隙度φ_initial；沙体积分数Y₂的初始值

(5)岩样半径a_rock、高度h、顶部压力p_top、底部压力p_base，孔隙压缩系数c_φ；

(6)单元个数N、时间步长h_t，积分时间t_total。

本例中，这些参数分别为

(1)破碎岩体粒径＝15～15mm，沙粒径＝0.074～0.25mm；

(2)水的质量密度m₁＝1000(kg/m³)，动力粘度μ₁＝1.01×10^-3(Pa·s)；

(3)沙的质量密度m₂＝2670kg/m³、稠度系数C＝3.26×10^-2(Pa·s^0.19)，幂指数n＝0.19；

(4)初始孔隙度φ_initial＝0.585，沙体积分数Y₂的初始值

(5)岩样半径a_rock＝0.05(m)，高度h＝0.6(m)，顶部压力p_top＝1.2(MPa)，底部压力p_base＝0，孔隙压缩系数c_φ＝0.326×10^-8；

(6)单元个数N＝10，时间步长h_t＝0.1×10^-3(s)，积分时间t_total＝0.36×10⁺⁰⁴(s)。

3.渗流场计算过程与结果

根据图2，利用Fortran语言编制破碎岩石孔隙中水驱沙过程渗流场计算程序。在该程序中输入孔隙度φⁱ(i＝1，2，…，7)、A ⁱ (i＝1，2，…，7)及岩样质量密度等参量的数值后，便可计算水驱沙过程的时变渗流场，并输出各时刻、各节点的体积分数Y₁和Y₂、渗流速度V₁和V₂、孔隙度φ和压力p、水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2。图3给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻体积分数Y₂在区间[0，h]上的分布情况。图4给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻渗流速度V₂在区间[0，h]上的分布情况。图5给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻渗流速度V₁在区间[0，h]上的分布情况。图6给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻流度I₁在区间[0，h]上的分布情况。图7给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻非Darcy流β因子β₁在区间[0，h]上的分布情况。图8给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻加速度系数c_a1在区间[0，h]上的分布情况。图9给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻流度I_e2在区间[0，h]上的分布情况。图10给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻非Darcy流β因子β₂在区间[0，h]上的分布情况。图11给出了t＝0(s)、t＝300(s)、t＝600(s)、……，t＝3600(s)时刻加速度系数c_a2在区间[0，h]上的分布情况。

Claims (6)

1.一种水驱沙过程中破碎岩石渗流场计算方法，在水驱沙过程中破碎岩石孔隙被水和沙两种介质占据，根据水、沙渗透过程中的质量守恒方程、动量守恒方程、孔隙压缩方程以及渗透性参量-沙体积分数关系式，建立了渗流过程中物理量关系框图；根据物理量关系框图设计渗流场计算方法，其特征在于，渗流场的计算是通过以下步骤实现的： 

步骤1)计算得到的若干级孔隙度下的渗透性参量-沙体积分数关系式，利用一元三点不等距Lagrange插住法分别计算出各节点在当前孔隙度下水相流度I₁、非Darcy流β因子β₁、加速度系数c_a1、沙相有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2、渗透性参量-沙体积分数Y₂关系式系数； 

步骤2)根据破碎岩石的质量、高度、半径和真密度，计算出初始时刻(t＝0)破碎试验的孔隙度，根据沙子的质量、真密度计算出初始时刻沙的体积分数，并计算出水的体积分数；设定渗透条件，写出节点压力和渗流速度的分布； 

步骤3)根据节点物理量计算单元上孔隙度对坐标的偏导数、压力对坐标的偏导数、体积分数对坐标的偏导数、渗流速度对坐标的偏导数；通过节点附近区域平衡方程的积分，分别计算出压力、孔隙度、体积分数和渗流速度对时间的偏导数，在通过时间步长上的积分，计算出下一时刻节点的压力、孔隙度、体积分数和渗流速度； 

步骤4)根据节点孔隙度和沙体积分数，利用一元三点不等距 Lagrange插住法分别计算出各节点处水相流度、非Darcy流β因子、加速度系数、沙相有效流度、非Darcy流β因子、加速度系数-沙体积分数关系式系数，并利用渗透性参量-沙体积分数关系式计算节点上渗透性参量的值，进入下一步循环。 

2.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，所述质量守恒方程、动量守恒方程和孔隙压缩方程的计算方法为： 

a.分别用下标1和2表示水相和沙相介质的物理量，破碎岩石孔隙中质量守恒方程为 

和 

其中，t为时间，x为空间坐标，p为水压，V渗流速度，m为质量密度，I为流度，I_e为有效流度，β为非Darcy流β因子，c_a为加速度系数，n为粘度指数； 

b.假设在破碎岩石孔隙中没有空气，即孔隙体积完全被水和沙占据，记水和沙的体积分数分别Y₁和Y₂，则有如下关系 

Y₁+Y₂＝1            (3) 

破碎岩石孔隙中水、沙的质量守恒方程分别为 

和 

将式(4)和式(5)相加，并利用式(3)，可以得到 

c.破碎岩石的孔隙具有压缩性，即孔隙度随压力变化，亦即孔隙压缩方程为 

由式(4)、式(5)和式(7)，可以得到 

式(1)～(3)、式(5)～(6)和式(8)构成六个变量p,φ,V₁,V₂,Y₁,Y₂的封闭方程组，此方程组中水的质量密度m₁、沙的质量密度m₂、破碎岩石的孔隙压缩系数c_φ和参考压力p₀下的孔隙度可近似当做常量来处理，而水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2则随孔隙度φ和体积分数Y₂变化。 

3.根据权利要求1或2所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，计算出各节点在当前孔隙度下的渗透性参量以及渗透性参量-沙体积分数关系的具体步骤是： 

步骤1.1)计算渗透性参量—沙体积分数关系 

对于某一特定孔隙度的破碎岩石，通过水驱沙试验和线性回归，得到水的流度I₁、非Darcy流β因子β₁、加速度系数c_a1，沙的有效流度I_2e、非Darcy流β因子β₂、加速度系数c_a2与沙子体积分数Y₂的关系， 即 

I₁＝A₁₀+A₁₁Y₂+A₁₂Y₂ ²                         (9a) 

lgβ₁＝A₂₀+A₂₁Y₂+A₂₂Y₂ ²              (9b) 

lgc_a1＝A₃₀+A₃₁Y₂+A₃₂Y₂ ²              (9c) 

I_2e＝A₄₀+A₄₁Y₂+A₄₂Y₂ ²               (9d) 

lgβ₂＝A₅₀+A₅₁Y₂+A₅₂Y₂ ²                 (9e) 

lgc_a2＝A₆₀+A₆₁Y₂+A₆₂Y₂ ²             (9f) 

将回归系数A_ij(i＝1,2,…,6；j＝0,1,2)组装成矩阵 

步骤1.2)计算渗透性参量—体积分数回归系数三维数组的构造对同一粒径破碎岩石设置m_p个不同的孔隙度 通过水驱沙试验得到m_p个不同孔隙度下的回归系数矩阵

将组装成三维矩阵A ^b，其中A ^b的第i层为A ⁱ，为了读数方便，也可将组装成二维分块矩阵A ^p，即 

或 

步骤1.3)渗透性参量的计算 

在水驱沙过程中，孔隙度随时间变化，当孔隙度介于φ₁和之间时，应用一元三点Lagrange插值的方法计算系数A_ij(i＝1,2,…,6；j＝0,1,2)，即 

再根据式(9)计算水相和沙相的渗透性参量。 

4.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步骤2的具体计算方法为： 

步骤2.1)设定边界条件和初始条件 

根据试验条件，岩样上端x＝0的压力保持恒定，即 

p|_x＝0＝p₀          (14) 

下端与大气连通，故 

p|_x＝h＝p₁＝0          (15) 

岩样上端x＝0没有沙子流入，故渗流速度 

V₂|_x＝0＝0          (16) 

步骤2.2)根据破碎岩石的质量M_rock、高度h、半径a_rock和破碎前的真密度m_rock，可以计算出初始时刻(t＝0)破碎试验的孔隙度 

步骤2.3)根据沙子的质量M_sand和真密度m_sand，可以计算出初始时刻沙的体积分数，即 

步骤2.4)利用式(3)可以计算出初始时刻水的体积分数，即 

步骤2.5)初始压力呈分布 

步骤2.6)初始时刻水渗流速度V₁

步骤2.7)沙渗流速度为 

。

5.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步骤3的具体计算方法为： 

步骤3.1)水驱沙过程中水沙两相流动的一维流场算法的核心是将物理量分为两类：单元物理量和节点物理量，将区间[0,h]均分为N个单元，单元长度为节点k坐标为 

单元k的左、右端点坐标分别为x_k和

节点物理量包括孔隙度φ_n(k)、体积分数Y_1n(k)、Y_2n(k)、渗流速度V_1n(k)、V_2n(k)、压力p_n(k)，k＝1,2,…,N+1，单元物理量包括孔隙度梯度体积分数梯度渗流速度梯度和压力梯度 k＝1,2,…,N； 

单元k上物理量对坐标偏导数的差分格式分别为 

步骤3.2)基于式(8)构造时刻t^j节点i上压力对时间的偏导数，即在区间上对式(8)积分，得到 

对于等距节点，式(21)可以简化为 

将式(22)对时间积分一步，得到 

步骤3.3)基于式(7)构造节点孔隙度的算法，即 

步骤3.4)基于式(5)和式(3)构造体积分数Y₂和Y₁的算法，即 

步骤3.5)基于式(1)和式(2)构造渗流速度V₁和V₂的算法，即 

。

6.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步骤4的具体计算方法为： 

步骤4.1)单元划分 

沿破碎岩样高度方向建立坐标轴Ox，坐标轴的正向与渗流速度方向相同，将破碎岩样均分为N个单元，单元长度为节点k坐 标为 

单元k的左、右端点坐标分别为x_k和x_k+1。 

步骤4.2)渗透参量-体积分数回归系数的确定 

步骤4.2.1)根据试验结果，给出m_p级孔隙度下的渗透参量-体积分数回归系数，孔隙度分别记为对于的渗透参量-体积分数回归系数构成的矩阵分别记为 

步骤4.2.2)将组装为分块矩阵 

步骤4.3)给出节点物理量的初始值 

步骤4.3.1)孔隙度的初始值 

步骤4.3.2)沙体积分数的初始值 

步骤4.3.3)水体积分数的初始值 

Y_1n(k)＝1-Y_1n(k),k＝1,2,…,N+1           (41) 

步骤4.3.4)压力的初始值 

步骤4.3.5)水渗流速度的初始值 

步骤4.3.6)沙渗流速度的初始值 

步骤4.4初始时刻的渗透性参量计算 

步骤4.4.1)利用试验结果得到的m_p级孔隙度φ⁽ⁱ⁾(i＝1,2,…,n_f)下的渗透性参量系数 

通过一元三点插值法计算初始孔隙度φ_initial下渗透性参量系数 

a_jk(j＝1,2,…,6；k＝1,2,3)； 

步骤4.4.2)根据式(9)计算节点上流度、非Darcy流β因子和加速度系数； 

步骤4.5)单元物理量计算 

步骤4.5.1)单元k上孔隙度对坐标的偏导数按式(24)计算，即 

步骤4.5.2)体积分数对坐标的偏导数和按式(25)和式(26)计算，即 

步骤4.5.3)渗流速度对坐标的偏导数和按式(27)和式(28)计算，即 

步骤4.5.4)压力对坐标的偏导数按式(29)计算，即 

步骤4.6)节点压力的计算 

步骤4.6.1)内部节点的压力按式(32)计算，即 

步骤4.6.2)边界节点1的压力按式(45)计算，即 

步骤4.6.3)边界节点N+1的压力按式(46)计算，即 

步骤4.7)节点孔隙度的计算，节点孔隙度按式(33)计算，即 

步骤4.8)节点体积分数Y₂的计算 

步骤4.8.1)内部节点的体积分数Y₂按式(34)计算，即 

步骤4.8.2)边界节点1的体积分数Y₂按式(47)计算，即 

步骤4.8.3)边界节点N+1的体积分数Y₂按式(48)计算，即 

步骤4.9)节点体积分数Y₁的计算 

节点体积分数Y₁按式(35)计算，即 

步骤4.10)节点渗流速度V₁的计算 

步骤4.10.1)内部节点的渗流速度V₁按式(36)计算，即 

i＝2,3,…,N                (36) 

步骤4.10.2)边界节点1渗流速度V₁按式(49)计算，即 

步骤4.10.3)边界节点N+1渗流速度V₁按式(50)计算，即 

步骤4.11)节点渗流速度V2的计算 

步骤4.11.1)内部节点的渗流速度V₂按式(37)计算，即 

i＝2,3,…,N              (37) 

步骤4.11.2)边界节点1渗流速度V₂按式(51)计算，即 

步骤4.11.3)边界节点N+1渗流速度V₂按式(52)计算，即 

。