一种水驱沙过程中破碎岩石渗流场计算方法
【专利摘要】本发明涉及一种水驱沙过程中破碎岩石渗流场计算方法，根据水、沙渗透过程中的质量守恒方程、动量守恒方程、孔隙压缩方程以及渗透性参量-沙体积分数关系式，建立了渗流过程中物理量关系框图；根据物理量关系框图设计渗流场计算方法，该方法考虑了水驱沙过程中，由于沙质量流失水引起的沙体积分数的变化；考虑了沙体积分数对水和沙的流度、非Darcy流β因子、加速度系数的影响；可计算出各节点渗透性参量；在算法设计中，水和沙的流度、非Darcy流β因子、加速度系数，随孔隙度和沙体积分数变化，时变渗流场利用异步迭代的方法计算，节省储存空间，运算速度快，收敛性和数值稳定性好，计算结果精确。
【专利说明】一种水驱沙过程中破碎岩石渗流场计算方法
【技术领域】
[0001] 本发明涉及煤矿地质灾害与防治领域，具体涉及一种水驱沙过程中破碎岩石渗流 场计算方法。
【背景技术】
[0002] 在毛乌素沙漠，煤层具有埋深小、地表沉沙厚的特点。在采煤过程中，地表的砂层 随地表水沿着煤层覆岩向采空区运移。为了研究煤层开采过程中由于覆岩破坏造成突水溃 沙的机理、防治突水溃沙的发生或减轻突水溃沙的危害，根据水驱沙试验得到的水相和沙 相渗透率、非Darcy流3因子、加速度系数随破碎岩石孔隙中沙子体积分数变化规律，建立 一种水驱沙过程中破碎岩石渗流场计算方法。
[0003] 目前，人们描述和模拟水驱沙过程中破碎岩石孔隙中水、沙两相介质的流动主要 参照油藏工程中水驱开发理论和方法。由于沙子本质上属于非连续介质，沙子在破碎岩 石孔隙中的迁移与油在油藏中渗流的运动机理、运动方式和阻力特性等方面存在显著的差 异。特别是，水和沙之间不存在分界面。因此，参照水驱油方法计算水驱沙过程中破碎岩石 的渗流场这种做法欠妥，构建一种有别于水驱油渗流场计算方法的、水沙在破碎岩石中水 沙两相渗流计算方法既是必要的。
【发明内容】
[0004] 本发明的目的是提供一种水驱沙过程中破碎岩石渗流场计算方法，根据水驱沙试 验得到的水相和沙相渗透率、非Darcy流P因子、加速度系数随破碎岩石孔隙中沙子体积 分数变化关系，建立水驱沙过程中破碎岩体中水沙渗流场的计算方法
[0005] 本发明是通过以下技术方案实现的：一种水驱沙过程中破碎岩石渗流场计算方 法，在水驱沙过程中破碎岩石孔隙被水和沙两种介质占据，根据水、沙渗透过程中的质量守 恒方程、动量守恒方程、孔隙压缩方程以及渗透性参量-沙体积分数关系式，建立渗流过程 中物理量关系框图；根据物理量关系框图设计渗流场计算方法。
[0006] 本发明的基本原理为：
[0007] 分别用下标1和2表示水相和沙相介质的物理量，破碎岩石孔隙中质量守恒方程 为
【权利要求】
1. 一种水驱沙过程中破碎岩石渗流场计算方法，在水驱沙过程中破碎岩石孔隙被水和 沙两种介质占据，根据水、沙渗透过程中的质量守恒方程、动量守恒方程、孔隙压缩方程以 及渗透性参量-沙体积分数关系式，建立了渗流过程中物理量关系框图；根据物理量关系 框图设计渗流场计算方法，其特征在于，渗流场的计算是通过以下步骤实现的： 步骤1)计算得到的若干级孔隙度下的渗透性参量-沙体积分数关系式，利用一元三点 不等距Lagrange插住法分别计算出各节点在当前孔隙度下水相流度I1、非Darcy流β因 子、加速度系数C al、沙相有效流度I2e、非Darcy流β因子β2、加速度系数ca2、渗透性参 量-沙体积分数Y 2关系式系数； 步骤2)根据破碎岩石的质量、高度、半径和真密度，计算出初始时刻（t = 0)破碎试 验的孔隙度，根据沙子的质量、真密度计算出初始时刻沙的体积分数，并计算出水的体积分 数；设定渗透条件，写出节点压力和渗流速度的分布； 步骤3)根据节点物理量计算单元上孔隙度对坐标的偏导数、压力对坐标的偏导数、体 积分数对坐标的偏导数、渗流速度对坐标的偏导数；通过节点附近区域平衡方程的积分，分 别计算出压力、孔隙度、体积分数和渗流速度对时间的偏导数，在通过时间步长上的积分， 计算出下一时刻节点的压力、孔隙度、体积分数和渗流速度； 步骤4)根据节点孔隙度和沙体积分数，利用一元三点不等距Lagrange插住法分别计 算出各节点处水相流度、非Darcy流β因子、加速度系数、沙相有效流度、非Darcy流β因 子、加速度系数-沙体积分数关系式系数，并利用渗透性参量-沙体积分数关系式计算节点 上渗透性参量的值，进入下一步循环。
2. 根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，所 述质量守恒方程、动量守恒方程和孔隙压缩方程的计算方法为： a. 分别用下标1和2表示水相和沙相介质的物理量，破碎岩石孔隙中质量守恒方程为
其中，t为时间，X为空间坐标,p为水压，V渗流速度,m为质量密度，I为流度，为有 效流度，β为非Darcy流β因子，Ca为加速度系数，η为粘度指数； b. 假设在破碎岩石孔隙中没有空气，即孔隙体积完全被水和沙占据，记水和沙的体积 分数分别Y1和Y2，则有如下关系 Y^Y2 = 1 (3) 破碎岩石孔隙中水、沙的质量守恒方程分别为
将式⑷和式（5)相加，并利用式（3)，可以得到
C.破碎岩石的孔隙具有压缩性，即孔隙度随压力变化，亦即孔隙压缩方程为
由式（4)、式（5)和式（7)，可以得到
式⑴?（3)、式（5)?（6)和式⑶构成六个变量p，（^VpLYpY2的封闭方程组， 此方程组中水的质量密度Hl1、沙的质量密度m2、破碎岩石的孔隙压缩系数C41和参考压力P tl 下的孔隙度4可近似当做常量来处理，而水的流度I1、非Darcy流β因子P1、加速度系数 Cal，沙的有效流度I2e、非Darcy流β因子β2、加速度系数c a2则随孔隙度Φ和体积分数Y2 变化。
3.根据权利要求1或2所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是， 计算出各节点在当前孔隙度下的渗透性参量以及渗透性参量-沙体积分数关系的具体步 骤是： 步骤1. 1)计算渗透性参量一沙体积分数关系 对于某一特定孔隙度的破碎岩石，通过水驱沙试验和线性回归，得到水的流度I1、非 Darcy流β因子P1、加速度系数Cal，沙的有效流度I2e、非Darcy流β因子β 2、加速度系 数ca2与沙子体积分数Y2的关系，即
将回归系数Aij (i = 1，2,…，6 ; j = 0, 1，2)组装成矩阵
步骤1. 2)计算渗透性参量一体积分数回归系数三维数组的构造对同一粒径破碎岩石 设置mp个不同的孔隙度H…，<，祕<4 <色，···，，之<<)，通过水驱沙试验得到mp个 不同孔隙度下的回归系数矩阵; 将左，…，▲?组装成三维矩阵#，其中Ab的第i层为Ai，为了读数方便，也可将
4.￥，···.屋?组装成二维分块矩阵Ap，即
步骤I. 3)渗透性参量的计算 在水驱沙过程中，孔隙度随时间变化，当孔隙度介于^和< 之间时，应用一元三点 Lagrange插值的方法计算系数Ai1. (i = 1，2,…，6 ;j = 0, 1，2)，即
再根据式（9)计算水相和沙相的渗透性参量。
4.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步 骤2的具体计算方法为： 步骤2. 1)设定边界条件和初始条件 根据试验条件，岩样上端X = 0的压力保持恒定，即 Plx = O = Po (14) 下端与大气连通，故 plx = h = Pi = 〇 (15) 岩样上端X = 0没有沙子流入，故渗流速度 V丄= 〇 = 〇 (16) 步骤2. 2)根据破碎岩石的质量Mrodt、高度h、半径arodt和破碎前的真密度Hirodt,可以计 算出初始时刻（t = 0)破碎试验的孔隙度
步骤2. 3)根据沙子的质量Msand和真密度msand，可以计算出初始时刻沙的体积分数，即
步骤2. 4)利用式（3)可以计算出初始时刻水的体积分数，即
步骤2. 5)初始压力呈分布
步骤2. 6)初始时刻水渗流速度V1
步骤2. 7)沙渗流速度为
5.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步 骤3的具体计算方法为： 步骤3. 1)水驱沙过程中水沙两相流动的一维流场算法的核心是将物理量分为两类： 单元物理量和节点物理量,将区间[0, h]均分为N个单元,单元长度戈
_节点k坐标 为
单元k的左、右端点坐标分别为xk和\ ; 节点物理量包括孔隙度Φη00、体积分数Yln(k)、Y2n(k)、渗流速度V ln(k)、V2n(k)、压力 pn(k)，k = 1，2,…，N+1，单元物理量包括孔隙度梯择
、体积分数梯度_
，渗 流速度梯度
单元k上物理量对坐标偏导数的差分格式分别为
步骤3. 2)基于式⑶构造时刻P节点i上压力对时间的偏导数，即在区间
上对式⑶积分，得到
对于等距节点，式（21)可以简化为
将式（22)对时间积分一步，得到
步骤3. 3)基于式（7)构造节点孔隙度的算法，即
步骤3. 4)基于式（5)和式（3)构造体积分数Y2和Y1的算法，即
步骤3. 5)基于式（1)和式（2)构造渗流速度V1和V2的算法，即
6.根据权利要求1所述的一种水驱沙过程中破碎岩石渗流场计算方法，其特征是，步 骤4的具体计算方法为： 步骤4. 1)单元划分 沿破碎岩样高度方向建立坐标轴〇x，坐标轴的正向与渗流速度方向相同，将破碎岩样 均分为N个单元，单元长度，
，节点k坐标为
单元k的左、右端点坐标分别为xk和xk+1。 步骤4. 2)渗透参量-体积分数回归系数的确定 步骤4. 2. 1)根据试验结果，给出mp级孔隙度下的渗透参量-体积分数回归系数，孔隙 度分别记为彡'、#、……^，对于的渗透参量-体积分数回归系数构成的矩阵分别记 为
步骤4. 2. 2)将！、f、……、^组装为分块矩阵
步骤4. 3)给出节点物理量的初始值 步骤4. 3.1)孔隙度的初始值
步骤4. 3. 2)沙体积分数的初始值
步骤4. 3. 3)水体积分数的初始值 Yln(k) = 1-Yln(k)，k = 1，2, ...，N+1 (41) 步骤4.3.4)压力的初始值
步骤4. 3. 5)水渗流速度的初始值
步骤4. 3. 6)沙渗流速度的初始值
步骤4. 4初始时刻的渗透性参量计算 步骤4. 4. 1)利用试验结果得到的mp级孔隙度Φ (i) (i = 1，2,…，nf)下的渗透性参量 系数
通过一元三点插值法计算初始孔隙度下渗透性参量系数 aJk(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)渗流速度对坐标的偏导数
K式（27)和式（28)计算，即 V ^ K. ΛΧ1 fV
步骤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)节点体积分数Y2的计算 步骤4. 8. 1)内部节点的体积分数Y2按式（34)计算，即
步骤4. 8. 2)边界节点1的体积分数Y2按式（47)计算，即 、-/1 -
-步骤4. 8. 3)边界节点N+1的体积分数Y2按式（48)计算，即
步骤4. 9)节点体积分数Y1的计算 节点体积分数Y1按式（35)计算，即
步骤4. 10)节点渗流速度V1的计算 步骤4. 10. 1)内部节点的渗流速度V1按式（36)计算，即
步骤4. 10. 2)边界节点1渗流速度V1按式（49)计算，即
步骤4. 10. 3)边界节点N+1渗流速度V1按式（50)计算，即
步骤4. 11)节点渗流速度V2的计算 步骤4. 11. 1)内部节点的渗流速度V2按式（37)计算，即
步骤4. 11. 2)边界节点1渗流速度V2按式（51)计算，即
步骤4. 11. 3)边界节点N+1渗流速度V2按式（52)计算，即
【文档编号】G06F19/00GK104392131SQ201410680250
【公开日】2015年3月4日 申请日期:2014年11月24日 优先权日:2014年11月24日
【发明者】丁焕德, 杜锋, 陈占清, 李虎民, 石建祥, 惠学齐, 朱南京 申请人:中国矿业大学, 河南理工大学, 华电煤业集团有限公司