本发明公开了一种地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，具体步骤：S1、采用压裂模拟器Fracpro PT模拟水力压裂过程中，获取水力裂缝扩展过程中的缝长、缝宽、缝高；S2、采集岩样，将岩样按照步骤S1获取的缝长与缝高等比例缩小至岩心尺寸，并加工成试样，利用人工劈裂的方式获得具有粗糙裂缝面的试样，然后激光扫描获得数字化的裂缝表面粗糙形貌特征参数；S3、获得具有地层尺度的粗糙裂缝形貌特征参数的裂缝表面；S4、获得由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝空间；S5、建立用于数值模拟的复杂流动空间网格系统；S6、模拟支撑剂在地层尺度的粗糙裂缝内的运动过程，最终获得水力压裂结束后地层尺度的粗糙裂缝内的支撑剂输送特征。
1.一种地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，其特征在于，包括以下步骤：
S1、采用压裂模拟器Fracpro PT模拟水力压裂过程中，获取水力裂缝扩展过程中的缝长、缝宽和缝高；
S2、采集岩样，将岩样按照步骤S1获取的缝长与缝高等比例缩小至岩心尺寸，并加工成试样，利用人工劈裂的方式获得具有粗糙裂缝面的试样，然后激光扫描获得数字化的裂缝表面粗糙形貌特征参数；
S3、将步骤S2获得的两个数字化的裂缝表面粗糙形貌特征参数按照步骤S1获得的缝长与缝高等比例还原到地层尺度，获得两个具有地层尺度的粗糙裂缝形貌特征参数的裂缝表面；
S4、在三维空间内通过“不错位靠近”的方式获得由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝空间，其中两个裂缝面之间的距离为步骤S1获得的缝宽；
S5、建立地层尺度的粗糙裂缝复杂流动空间几何模型，并通过导入成熟的网格处理软件数值化，建立用于数值模拟的复杂流动空间网格系统；
S6、将步骤S5建立的复杂流动空间网格系统导入Open Foam开源流体力学计算软件，选择将固相支撑剂拟流体化的两流体方法，设置相应的边界条件，模拟支撑剂在地层尺度的粗糙裂缝内的运动过程，最终获得水力压裂结束后地层尺度的粗糙裂缝内的支撑剂输送特征；
步骤S6包括以下子步骤：S61、将导出的“.msh”格式文件导入Open Foam开源流体力学计算软件中；
S62、根据裂缝模型几何尺寸，将现场施工排量转换为线速度并作为模拟初始速度，同时将现场施工砂比或砂浓度转化为液、固两相体积分数并作为模型初始相分数场；
S63、使用广义颗粒温度Θ表征固相运动，广义颗粒温度Θ的计算公式如下：式中，Θt——颗粒温度，即颗粒平移速度波动产生的动能相关项；
Θc——固相构型温度，即固相形变速度波动产生的动能相关项；
Js——单位体积内因颗粒非弹性碰撞产生的能量耗散，J/(m·s)；
Πs——单位体积内因相间能量交换产生的能量耗散，J/(m·s)；
S64、计算固相剪切黏度μs，计算公式如下：式中，μs,k——颗粒传输和碰撞产生的剪切黏度，mPa·s；
μs,f——固相接触摩擦产生的剪切黏度，mPa·s；
μs——考虑颗粒间流体影响的固相黏性系数，mPa·s；
Fr、r'、s'——摩擦模型系数，通过实验确定(默认值分别取0.05、2和5)；
S65、根据如下公式分别求解液、固两相的速度场和相分数场，最终通过固相分数的分布变化反映支撑剂在裂缝中的输运情况：式中：εf——固相体积分数，小数；
2.如权利要求1所述的地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，其特征在于，所述步骤S5包括以下两个子步骤：S51、将步骤S4获得的三维复杂裂缝空间的特征数据，导入成熟的网格处理平台ICEM，建立地层尺度的粗糙裂缝复杂流动空间几何模型；
S52、继续在ICEM软件平台中，将地层尺度的粗糙裂缝复杂流动空间几何模型进行合理的网格划分，建立用于数值模拟的复杂流动空间网格系统，导出为“.msh”的格式文件。
3.如权利要求1所述的地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，其特征在于，步骤S1具体操作是：根据现场实际参数，将目标压裂井层的井筒资料、地质资料、压裂液资料、施工参数输入到压裂模拟器Fracpro PT中，模拟水力压裂过程中储层水力裂缝延伸扩展过程，获取水力裂缝扩展过程中的缝长、缝宽和缝高。
4.如权利要求1所述的地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，其特征在于，步骤S2，采集的岩样是目的层的井下岩心或油气储层同层位的露头岩样。
5.如权利要求4所述的地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，其特征在于，步骤S2中，岩样等比例缩小的具体公式如下：式中：Lformation——通过模拟获得的水力裂缝扩展过程中的缝长；
Hformation——通过模拟获得的水力裂缝扩展过程中的缝高；
一种地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法
技术领域
[0001] 本发明涉及油气田开发工程领域，尤其是水力压裂改造过程中的一种地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法。
背景技术
[0002] 自美国于1947年首次开展水力压裂作业开始，水力压裂技术逐渐发展成为众多储层增产改造的“破冰利刃”，尤其在最近十数年对致密砂岩、页岩及煤层等非常规储层的有
效改造和经济开发中，该项技术已经无可替代。水力压裂施工过程中，通过向地层注入高压
流体实现压开裂缝，沟通储层并为油气流动提供高速通道；施工结束后，缝内压力降低导致
裂缝壁面逐渐闭合，需要依靠早期混在压裂液中的高强度支撑剂支撑裂缝。支撑剂在裂缝
中的充填使裂缝保持较高的导流能力，即保证了“油气高速公路”畅通，例如卢聪,郭建春,
王文耀,邓燕，刘登峰.支撑剂嵌入及对裂缝导流能力损害的实验[J].天然气工业,2008,28
(2):99‑101.。因此，支撑剂在缝内的运移规律与分布形态对施工效果的影响至关重要。
[0003] 缝内支撑剂运移规律与分布形态研究分为理论分析、物模实验和数值模拟三大类。由于大多现存理论已经滞后于现场应用，物模实验和数值模拟成为主要研究手段并共
同推进理论发展。实验研究表明缝内支撑剂运动受混砂液黏度、排量、颗粒浓度以及裂缝宽
度的影响。地层岩石在高压流体作用下受拉张或剪切破坏，产生粗糙裂缝壁面、缝宽非均匀
分布，而实验装置大多采用光滑壁面、均匀缝宽分布。虽然粘合岩板(例如LiuY J,Sharma M
M.Effect of Fracture Width and Fluid Rheology on Proppant Settling and
Retardation[C].SPE 96208,2005.)和3D打印技术(例如专利CN201910322572.4)被引入用
于考虑粗糙裂缝壁面对支撑剂运动的影响，受实验规模、成本和研究条件限制，需要大量结
合数值模拟工作。
[0004] 而且，现有技术忽略了地层尺度的裂缝粗糙面对支撑剂运动影响，容易导致支撑剂运移规律和堆积形态认识出现偏差，影响非常规油气储层水力压裂设计和效果评价。
发明内容
[0005] 本发明的目的是针对现有技术忽略了地层尺度的裂缝粗糙面对支撑剂运动影响，容易导致支撑剂运移规律和堆积形态认识出现偏差的问题，提供一种地层尺度下的粗糙裂
缝内支撑剂输送数值模拟方法。
[0006] 本发明提供的地层尺度下的粗糙裂缝内支撑剂输送数值模拟方法，包括以下步骤：
[0007] S1、根据现场实际参数，将目标压裂井层的井筒资料、地质资料、压裂液资料、施工参数输入到压裂模拟器Fracpro PT中，模拟水力压裂过程中储层水力裂缝延伸扩展过程，
获取水力裂缝扩展过程中的缝长、缝宽、缝高。
[0008] S2、采集目的层的井下岩心或油气储层同层位的露头岩样，按照步骤S1获得的缝长与缝高比例等比例缩小到岩心尺寸，并加工成试样；参照专利ZL 201810203404.9中描述
的方法，利用人工劈裂的方式将试样剖开，获得具有两个粗糙裂缝面的试样；并采用三维激
光扫描方法获得数字化的裂缝表面粗糙形貌特征参数。具体按比例缩小的公式如下；
[0009]
[0010] 式中：Lformation——通过模拟获得的水力裂缝扩展过程中的缝长；
[0011] Hformation——通过模拟获得的水力裂缝扩展过程中的缝高；
[0012] Lcore——待加工岩心的目标试样长度；
[0013] Hcore——待加工岩心的目标试样宽度。
[0014] S3、将S2获得的两个数字化的裂缝表面粗糙形貌特征参数，按照步骤S1获得的缝长与缝高比例，等比例放大还原到地层尺度，获得两个具有地层尺度的粗糙裂缝形貌特征
参数的裂缝表面。具体方法是：(1)将步骤S2获得的两个数字化的裂缝表面粗糙形貌特征参
数，沿着裂缝面长度、宽度、高度方向建立三维直角坐标系，并建立裂缝面特征参数矩阵；
(2)按照长度方向对应缝长、宽度方向对应缝高的方式，以步骤S1中缩小比例等比例还原到
地层尺度，获得两个具有地层尺度的粗糙裂缝形貌特征参数的裂缝表面。
[0015] S4、将步骤S3获得的两个地层尺度的粗糙裂缝形貌特征参数的裂缝表面，在三维空间内通过“不错位靠近”的方式获得由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝
空间，其中两个裂缝面之间的距离为步骤S1获得的缝宽。
[0016] S5、将步骤S4所获得的由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝空间的特征数据，建立地层尺度的粗糙裂缝复杂流动空间几何模型，并通过导入成熟的网格处理
软件数值化，建立用于数值模拟的复杂流动空间网格系统。包括下列步骤：
[0017] (1)将所获得的由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝空间特征数据，导入成熟的网格处理平台ICEM，建立地层尺度的粗糙裂缝复杂流动空间几何模型；
[0018] (2)继续在ICEM软件平台中，将地层尺度的粗糙裂缝复杂流动空间几何模型进行合理的网格划分，建立用于数值模拟的复杂流动空间网格系统，导出为“.msh”的格式文件。
[0019] S6、将步骤S5获得的复杂流动空间网格系统导入Open Foam开源流体力学计算软件，选择将固相支撑剂拟流体化的两流体(Two Fluid Method，TFM)方法，设置相应的边界
条件，模拟支撑剂在地层尺度的粗糙裂缝内的运动过程，最终获得水力压裂结束后地层尺
度的粗糙裂缝内的支撑剂输送特征。具体操作步骤如下：
[0020] (1)将导出复杂流动空间网格系统的“.msh”格式文件导入Open Foam开源流体力学计算软件中；
[0021] (2)根据裂缝模型几何尺寸，将现场施工排量转换为线速度并作为模拟初始速度，同时将现场施工砂比或砂浓度转化为液、固两相体积分数并作为模型初始相分数场；
[0022] (3)传统的固相拟流体化(视为连续相)方法仅使用颗粒温度(即颗粒平移速度波动产生的动能相关项)Θt表征固相运动，在此基础上引入固相构型温度Θc(即固相形变速
度波动产生的动能相关项)以捕捉中‑高固相浓度条件下支撑剂的运动过程，通过公式(2)
计算广义颗粒温度Θ：
[0023]
[0024] 式中，εs——固相体积分数，小数；
[0025] ρs——支撑剂密度，kg/m3；
[0026] κs——颗粒能量扩散系数，J·s/m3；
[0027] Js——单位体积内因颗粒非弹性碰撞产生的能量耗散，J/(m3·s)；
[0028] Πs——单位体积内因相间能量交换(相互作用)产生的能量耗散，J/(m3·s)。
[0029] us——颗粒平均速度，m/s；
[0030] ξs——固相应力张量，Pa；
[0031] β——常数，通过实验获得；
[0032] D——常数，与颗粒内摩擦角有关；
[0033] t——系综平均的时间范围，s；
[0034] d——颗粒粒径，m。
[0035] 其中，输入模型的初始颗粒温度Θt通过公式(3)计算获得，颗粒瞬时速度νs和颗粒平均速度us通过PIV大型平板实验测试获得(详见专利ZL201910087439.5中的方法)：
[0036]
[0037] 式中，cs——颗粒平移速度波动，m/s(<>表示系综平均)。
[0038] 固相应力张量ξs通过公式(4)计算：
[0039]
[0040] 式中，ps——固相压力，Pa；
[0041] μb——固相体积黏度，mPa·s；
[0042] I——单位张量；
[0043] τs——固相剪切应力，Pa；
[0044] ep——颗粒碰撞恢复系数，通过实验获得；
[0045] μs——固相剪切黏度，mPa·s；
[0046] us——颗粒平均速度，m/s；
[0047] T——转置符号；
[0048] 根据颗粒动理论，固相压力ps和固相体积黏度μb通过公式(5)计算：
[0049]
[0050] 式中：εs,max——固相最大堆积分数，小数；
[0051] g0——颗粒接触径向分布函数；
[0052] εs——固相体积分数，小数；
[0053] μ——平均固相黏度，mPa·s；
[0054] ρs——支撑剂密度，kg/m3；
[0055] d——颗粒粒径，m。
[0056] (4)基于土力学准静态摩擦理论(Srivastava A,Sundaresan S.Analysis of a Frictional–kinetic Model for Gas–particle Flow[J].PowderTechnology,2003,72‑
85.)，固相剪切黏度μs考虑了颗粒传动、碰撞和摩擦三方面作用，其中颗粒传动和碰撞的影
响考虑在μs,k中，摩擦的影响考虑在μs,f中，由公式(6)计算：
[0057]
[0058] 式中：μs,k——颗粒传输和碰撞产生的剪切黏度，mPa·s；
[0059] μs,f——固相接触摩擦产生的剪切黏度，mPa·s；
[0060] α——固相黏性常数，取1.6；
[0061] μs*——考虑颗粒间流体影响的固相黏性系数，mPa·s；
[0062] g0——颗粒接触径向分布函数；
[0063] εs——固相体积分数，小数；
[0064] μb——固相体积黏度，mPa·s；
[0065] ρs——支撑剂密度，kg/m3；
[0066] μ——平均固相黏度，mPa·s；
[0067] βc——固、液动量交换系数；
[0068] ps——固相压力，Pa；
[0069] Ss——固相剪切形变量，s‑1；
[0070] d——颗粒粒径，m；
[0071] pc——固相临界压力，Pa；
[0072] n——屈服曲线形状因子；
[0073] Fr、r'、s'——摩擦模型系数，通过实验确定(默认值分别取0.05、2和5)；
[0074] φ——颗粒内摩擦角，通过实验测定；
[0075] ——产生摩擦作用的最小固相体积分数，取50％；
[0076] ——固相最大堆积体积分数，取63％；
[0077] us——颗粒平均速度，m/s。
[0078] (5)完成对固相支撑剂拟流体化的方程后，通过公式(7)所示的动量方程分别求解液、固两相的速度场和相分数场，最终通过固相分数的分布变化反映支撑剂在裂缝中的输
运情况：
[0079]
[0080] 式中：εf——液相体积分数，小数；
[0081] ρf——液相密度，kg/m3；
[0082] uf,i、uf,j——液相流速，m/s；其中，i、j是张量表示法；
[0083] xj——位移，m(张量表示方法)；
[0084] g——重力加速度，取9.8m/s2；
[0085] ξ——液相应力张量，Pa；
[0086] nfi'——单位面积流‑固相互作用，N/m3；
[0087] ρs——支撑剂密度，kg/m3；
[0088] εs——固相体积分数，小数；
[0089] us,i、us,j——固相速度，m/s，其中，i、j是张量表示法；
[0090] ξs——固相应力张量，Pa。
[0091] 与现有技术相比，本发明的有益之处在于：
[0092] 本发明的方法采用数值模拟手段研究水力压裂过程中支撑剂在粗糙裂缝中的运移规律和分布形态。该方法原理可靠，能考虑真实水力裂缝形态和裂缝中支撑剂的运动特
征，为非常规油气储层水力压裂设计和效果评价提供指导。
[0093] 本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
[0094] 图1是本发明中地层尺度的具有粗糙裂缝形貌特征参数的裂缝表面。
[0095] 图2是本发明中裂缝面“不错位靠近”方式示意图。
[0096] 图3是本发明中地层尺度的粗糙裂缝复杂流动空间网格。
[0097] 图4是本发明中地层尺度的粗糙裂缝液相和固相流场示意图。
[0098] 图5是本发明中地层尺度的粗糙裂缝内支撑剂分布模拟结果。
具体实施方式
[0099] 以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0100] 本发明的方法应用在具体的压裂井层案例中，具体步骤如下：
[0101] (1)根据现场实际参数，将目标压裂井层的井筒资料、地质资料、压裂液资料、施工参数输入到压裂模拟器Fracpro PT中，模拟水力压裂过程中储层水力裂缝延伸扩展过程，
获得水力裂缝扩展过程中的缝长为80m，缝宽为5mm，缝高为20m。
[0102] (2)采集油气储层同层位的露头岩样，将岩样按照水力裂缝扩展过程中的缝长与缝高的比例4：1(缝长：缝高＝80：20)等比例缩小，通过加工切割成长度为8cm，宽度为2cm的
试样。
[0103] (3)参照专利ZL 201810203404.9中描述的方法，利用人工劈裂的方式将试样剖开，获得具有两个粗糙裂缝面的试样，并采用三维激光扫描方法获得数字化的裂缝表面粗
糙形貌特征参数。
[0104] (4)将数字化的裂缝表面粗糙形貌特征参数，按照8cm：2cm＝80m：20m，即1000倍的放大比例还原到地层尺度，获得两个具有地层尺度的粗糙裂缝形貌特征参数的裂缝表面，
如图1所示。
[0105] (5)将获得的两个地层尺度的粗糙裂缝形貌特征参数的裂缝表面，在三维空间内通过“不错位靠近”的方式(如图2所示)获得由两个地层尺度的粗糙裂缝面组成的三维复杂
裂缝空间，其中两个裂缝面之间的距离为水力裂缝扩展过程中的缝宽，即为5mm。
[0106] (6)将获得的由两个地层尺度的粗糙裂缝面组成的三维复杂裂缝空间特征数据，建立地层尺度的粗糙裂缝复杂流动空间几何模型，并通过导入成熟的网格处理软件数值
化，建立用于数值模拟的复杂流动空间网格系统，如图3所示，并导出为“.msh”的格式文件。
[0107] (7)将导出复杂流动空间网格系统的“.msh”格式文件导入Open Foam开源流体力学计算软件中，选择将固相支撑剂拟流体化的两流体(Two Fluid Method，TFM)方法，模拟
粗糙裂缝内支撑剂的输送。
[0108] (8)并根据裂缝模型几何尺寸，将现场施工排量转换为线速度并作为模拟初始速度，同时将现场施工砂比或砂浓度转化为液、固两相体积分数并作为模型初始相分数场，具
体边界设置参数如表1所示。
[0109] (9)输入表2中的模拟基础数据，模拟支撑剂在地层尺度的粗糙裂缝内的运动过程，分别计算地层尺度的粗糙裂缝液相流场和固相流场。由于本案例中液相流场和固相流
场的速度场差异不大，因此用一个流场示意图(图4)展示。图中箭头的长度不同代表不同的
速度，箭头越长，表示速度越大。最终通过固相分数的分布变化反映支撑剂在裂缝中的输运
情况，获得水力压裂结束后地层尺度的粗糙裂缝内的支撑剂输送特征(如图5所示)。
[0110] 表1模型边界条件设置
[0111]
[0112]
[0113] 表2模型参数设置
[0114]
[0115] 以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人
员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰
为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质
对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
