本发明实施例涉及一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法，包括步骤S10，获取地层参数、施工参数、完井参数；步骤S20：设定初始的缝压和缝宽；步骤S30：建立耦合井筒流动的水平井分段多簇压裂裂缝扩展计算模型；步骤S40：建立时间步递增模型；步骤S50：建立压裂缝内支撑剂运移计算模型；步骤S60：判断该时刻是否加入支撑剂，若是，则采用高阶加权本质不震荡有限差分方法对压裂裂缝扩展过程中支撑剂运移进行计算，得到缝内支撑剂分布；若不是，则进入下一步；步骤S70：判断时间步是否达到施工时间，若没有则返回步骤S4，否则结束计算。
水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模
方法
技术领域
[0001] 本发明实施例涉及油气田开发技术领域，特别涉及一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法。
背景技术
[0002] 非常规油气藏的开发依赖于水平井分段多簇压裂的发展和应用。水平井分段多簇压裂通过“分簇射孔”技术实现一段多裂缝扩展，从而进一步提高施工效率、节约施工成本。多簇压裂后在地层会形成具有一定导流能力的填砂裂缝，裂缝中的支撑剂分布形态将直接影响裂缝导流能力和油气井产能，故对压裂缝中的支撑剂运移规律进行研究十分必要。然而，目前对于水平井分段多簇压裂裂缝扩展、多裂缝非均匀扩展过程缝内支撑剂分布规律等问题研究并不充分，特别是水平井分段多簇压裂裂缝扩展与支撑剂运移一体化计算方法研究鲜有报道。
发明内容
[0003] 本发明实施方式的目的在于提供一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法，旨在提供水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法。
[0004] 为解决上述技术问题，本发明的实施方式提供了一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法,包括：
[0005] 步骤S10，获取地层参数、施工参数、完井参数；
[0006] 步骤S20：设定初始的缝压和缝宽；
[0007] 步骤S30：建立耦合井筒流动的水平井分段多簇压裂裂缝扩展计算模型；
[0008] 步骤S40：建立时间递增模型；
[0009] 步骤S50：建立压裂缝内支撑剂运移计算模型；
[0010] 步骤S60：判断该时刻是否加入支撑剂，若是，则采用高阶加权本质不震荡有限差分方法对压裂裂缝扩展过程中支撑剂运移进行计算，得到缝内支撑剂分布；若不是，则进入下一步；
[0011] 步骤S70：判断时间是否达到施工时间，若没有则返回步骤S4，否则结束计算。
[0012] 优选地，所述步骤S10包括：
[0013] 获取地应力分布、岩石力学参数、注入程序、加砂程序、液体粘度、排量、注入时间、分簇数量、裂缝间距、射孔数、射孔直径。
[0014] 优选地，所述步骤S200的具体过程如下：
[0015] 假设每个单元的宽度为常数，根据叠加原理，得到初始t时刻裂缝内压力与宽度的离散方程为
[0016]
[0017] 其中，N为单元总数；
[0018] p(x,t)为缝内流体压力，MPa；
[0019] σh(x)为最小主应力，MPa；
[0020] c(x,ξi)为平面三维裂缝核函数，MPa/m；
[0021] w(ξi)为裂缝宽度，m；
[0022] ξi为裂缝面网格点的坐标；
[0023] x为计算点的坐标；
[0024] 通过求解式(1)可得到初始时刻的缝压与缝宽；
[0025] 其中式(1)中的矩阵形式为
[0026] p‑σh＝Cw (2)
[0027] 其中，p为缝内压力向量，MPa；
[0028] σh为最小主应力向量，MPa；
[0029] C为影响系数矩阵，MPa/m；
[0030] w为缝宽向量。
[0031] 优选地，所述步骤S300包括：
[0032] 建立井筒流动模型；
[0033] 建立裂缝宽度与压力的流固耦合方程；
[0034] 建立裂缝动边界。
[0035] 优选地，所述井筒流动模型的建立过程如下：
[0036] 井筒到各簇裂缝的分流满足质量守恒和压力连续条件，式(3)和式(4)构成井筒流动模型，式(3)和式(4)如下：
[0037]
[0038] pw＝pp,k+pc,k+pin,k (4)
[0039] 式(2)中
[0040] 式(2)中
[0041] 其中，Qt为一段开启Nf条裂缝的情况，注入总流量，Qi为第i条分支裂缝流量；
[0042] pw为井底压力,MPa；
[0043] k＝1，2，......，Nf；
[0044] pp,k为k裂缝的射孔摩阻，MPa；
[0045] pc,k为井口到k裂缝的井筒流动摩阻，MPa；
[0046] pin,k为k裂缝的入口压力，MPa；
[0047] Qk为k裂缝的入口流量，m3/s；
[0048] ρ为液体密度，kg/m3；
[0049] nk为k射孔簇的射孔数量；
[0050] dk为k射孔簇的射孔直径，mm；
[0051] K为射孔磨蚀修正系数，无因次；
[0052] fc为沿程摩阻系数，无因次；
[0053] Dw为压裂管柱内径，m；lk为井口到k裂缝的管柱长度，m；
[0054] ε为压裂管柱的内壁粗糙度，m；
[0055] Vw为井筒内液体流速，m/s；
[0056] Re为雷诺数，Re＝DwρVw/μ；
[0057] μ为液体粘度，mPa·s；
[0058] ρ为液体密度，kg/m3。
[0059] 优选地，所述裂缝宽度与压力的流固耦合方程的建立过程如下：
[0060] 岩石变形通过边界元方法计算，计算式如下：
[0061] p‑σh＝Cw (7)
[0062] 缝内流动满足质量守恒和层流方程，对缝内流动方程的空间项进行有限体积离散，得到流动方程的一阶微分方程形式：
[0063]
[0064] 将式(7)带入式(8)，得到所述裂缝宽度与压力的流固耦合方程：
[0065]
[0066] 井筒模型得到各簇裂缝进液流量后，求解式(9)可得到更新单元宽度和压力分布；同时式(9)得到的井底压力与井筒模型进行对比，直至收敛；
[0067] 其中，p为缝内压力向量，MPa；σh为最小主应力向量，MPa；C为影响系数矩阵，MPa/m；w为缝宽向量，m；θ为系数，0≤θ≤1；w0、p0分别为上一步宽度和压力分布；A(w)为流动方程的系数矩阵；S为源汇项。
[0068] 优选地，所述裂缝动边界的建立过程：
[0069] 建立尖端扩展判断准则，过程如下：
[0070]
[0071] 其中，d为距尖端的距离，m；
[0072] K'＝4(2/π)0.5KIc，MPa·m0.5；
[0073] KIc为I型断裂韧性，MPa·m0.5；
[0074] E'为平面应变杨氏模量，E'＝E/(1‑v2)，MPa；
[0075] v‑岩石泊松比，无因次；
[0076] w为缝宽向量；
[0077] 式(6)可得到单元发生扩展的临界宽度，通过比较当前时刻尖端单元的宽度是否达到临界宽度，若达到则增加单元，否则单元数量不增加。
[0078] 优选地，所述建立时间递增模型包括：
[0079] 时间步计算公式为：
[0080]
[0081] 时间递增为：
[0082] tm＝t+Δt (12)
[0083] 其中，λ为松弛因子，λ∈(0,1]；
[0084] 为所有待检查单元的最大扩展速度，m/s；
[0085] lf为裂缝长度，m；
[0086] ε＝0.08；
[0087] smax＝301；
[0088] ΔtE为显式算法的时间步长，s；
[0089] Δx为x方向单元大小，m；
[0090] Δy为y方向单元大小，m。
[0091] 优选地，所述压裂缝内支撑剂运移计算模型的建立过程包括：
[0092] 裂缝内携砂液流动和支撑剂运移的质量守恒方程为
[0093]
[0094]
[0095]
[0096]
[0097] 其中，qs、qp分别为携砂液和支撑剂的缝内流量，m3/s；
[0098] q1为滤失项；
[0099] p为封内流体压力,MPa；
[0100] 为归一化的加砂浓度，无因次；
[0101] 为归一化的支撑剂浓度，
[0102] φ为支撑剂体积浓度，无因次；
[0103] φm为最大容许浓度，无因次；
[0104] Qk为k裂缝的入口流量，m3/s；
[0105] a为支撑剂半径，m；
[0106] g为重力加速度，m/s2；
[0107] ρp、ρf分别为支撑剂和压裂液密度，kg/m3；
[0108] ez为z轴方向单位向量；
[0109] 为描述携砂液流动的无因次函数；
[0110] 分别为描述支撑剂对流和沉降的无因次函数；
[0111] x为缝内任一点的坐标向量；
[0112] xin,k为第k裂缝的入口坐标向量；
[0113] δ(x‑xin,k)为Delta函数，输入值为0时取1，其他值为0；
[0114] B为阻塞函数，描述支撑剂在缝内砂堵现象。
[0115] 优选地，所述步骤S70包括：
[0116] 判断时间递增是否达到施工时间，若没有则返回步骤S40，否则结束计算。
[0117] 本发明了一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法，首先建立耦合井筒流动的水平井分段多簇压裂裂缝扩展模型，然后建立支撑剂运移模型，并将其与裂缝扩展模型进行耦合，接着采用加权本质不震荡有限差分方法对耦合方程进行求解，最后得到缝内支撑剂分布。
[0118] 进一步地，本发明提出了裂缝扩展与支撑剂运移一体化计算方法，其中裂缝扩展是平面三维多裂缝模型，支撑剂运移通过高阶加权本质不震荡方法求解，可解决支撑剂不连续分布问题。裂缝扩展模型为平面三维模型，精度高，而且缝内支撑剂运移可以自动实现砂堵、沉降、流态转变等，并基于高精度方法求解支撑剂分布问题，方法更准确可靠。该方法可应用于水力压裂加砂设计，是提高压裂设计效果和油井增产的有效设计工具。
附图说明
[0119] 一个或多个实施例通过与之对应的附图中的图片进行示例性说明，这些示例性说明并不构成对实施例的限定，附图中具有相同参考数字标号的元件表示为类似的元件，除非有特别申明，附图中的图不构成比例限制。
[0120] 图1为压裂井与光纤监测井的物理模型示意图；
[0121] 图2为裂缝扩展形态的示意图；
[0122] 图3为支撑剂分布的示意图。
[0123] 本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0124] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0125] 需要说明，若本发明实施例中有涉及方向性指示(诸如上、下、左、右、前、后……)，则该方向性指示仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。
[0126] 另外，若本发明实施例中有涉及“第一”、“第二”等的描述，则该“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。另外，各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0127] 本发明提供一种水平井分段多簇压裂裂缝扩展与支撑剂运移一体化的建模方法,包括：
[0128] 步骤S10，获取地层参数、施工参数、完井参数；
[0129] 具体地，所述步骤S10包括：
[0130] 获取地应力分布、岩石力学参数、注入程序、加砂程序、液体粘度、排量、注入时间、分簇数量、裂缝间距、射孔数、射孔直径。
[0131] 在本实施例中，具体参数取值参见表1，在其他实施例中，可以根据实际情况而定。
[0132] 表1基本参数表
[0133]
[0134]
[0135] 步骤S20：设定初始的缝压和缝宽；
[0136] 具体地，所述步骤S200的具体过程如下：
[0137] 假设每个单元的宽度为常数，根据叠加原理，得到初始t时刻裂缝内压力与宽度的离散方程为
[0138]
[0139] 其中，N为单元总数；
[0140] p(x,t)为缝内流体压力，MPa；
[0141] σh(x)为最小主应力，MPa；
[0142] c(x,ξi)为平面三维裂缝核函数，MPa/m；
[0143] w(ξi)为裂缝宽度，m；
[0144] ξi为裂缝面网格点的坐标；
[0145] x为计算点的坐标；
[0146] 通过求解式(1)可得到初始时刻的缝压与缝宽；
[0147] 其中式(1)中的矩阵形式为
[0148] p‑σh＝Cw (2)
[0149] 其中，p为缝内压力向量，MPa；
[0150] σh为最小主应力向量，MPa；
[0151] C为影响系数矩阵，MPa/m；
[0152] w为缝宽向量，MPa。
[0153] 步骤S30：建立耦合井筒流动的水平井分段多簇压裂裂缝扩展计算模型；
[0154] 具体地，所述步骤S300包括：
[0155] 步骤S301：建立井筒流动模型；
[0156] 具体地，所述步骤S301所述井筒流动模型的建立过程如下：
[0157] 井筒到各簇裂缝的分流满足质量守恒和压力连续条件，式(3)和式(4)构成井筒流动模型，一段开启Nf条裂缝的情况，注入总流量等于各分支流量之和，式(3)和式(4)如下：
[0158]
[0159] 每条裂缝构成的分支回路具有相同压力降，
[0160] pw＝pp,k+pc,k+pin,k (4)
[0161] 式(2)中
[0162] 式(2)中
[0163] 其中，Qt为一段开启Nf条裂缝的情况，注入总流量，Qi为第i条分支裂缝流量；
[0164] pw为井底压力,MPa；
[0165] k＝1，2，......，Nf；
[0166] pp,k为k裂缝的射孔摩阻，MPa；
[0167] pc,k为井口到k裂缝的井筒流动摩阻，MPa；
[0168] pin,k为k裂缝的入口压力，MPa；
[0169] Qk为k裂缝的入口流量，m3/s；
[0170] ρ为液体密度，kg/m3；
[0171] nk为k射孔簇的射孔数量；
[0172] dk为k射孔簇的射孔直径，mm；
[0173] K为射孔磨蚀修正系数，无因次；
[0174] fc为沿程摩阻系数，无因次；
[0175] Dw为压裂管柱内径，m；lk为井口到k裂缝的管柱长度，m；
[0176] ε为压裂管柱的内壁粗糙度，m；
[0177] Vw为井筒内液体流速，m/s；
[0178] Re为雷诺数，Re＝DwρVw/μ；
[0179] μ为液体粘度，mPa·s；
[0180] ρ为液体密度，kg/m3。
[0181] 步骤S302：建立裂缝宽度与压力的流固耦合方程；
[0182] 具体地，所述裂缝宽度与压力的流固耦合方程的建立过程如下：
[0183] 将裂缝所在平面离散为一系列矩形单元，岩石变形通过边界元方法计算，计算式如下：
[0184] p‑σh＝Cw (7)
[0185] 缝内流动满足质量守恒和层流方程，对缝内流动方程的空间项进行有限体积离散，得到流动方程的一阶微分方程形式：
[0186]
[0187] 将式(7)带入式(8)，得到所述裂缝宽度与压力的流固耦合方程：
[0188]
[0189] 井筒模型得到各簇裂缝进液流量后，求解式(9)可得到更新单元宽度和压力分布；同时式(9)得到的井底压力与井筒模型进行对比，直至收敛；
[0190] 其中，p为缝内压力向量，MPa；σh为最小主应力向量，MPa；C为影响系数矩阵，MPa/m；w为缝宽向量，m；θ为系数，0≤θ≤1；w0、p0分别为上一步宽度和压力分布；A(w)为流动方程的系数矩阵；S为源汇项。
[0191] 步骤S303：裂缝动边界
[0192] 具体地，建立尖端扩展判断准则，过程如下：
[0193]
[0194] 其中，d为距尖端的距离，m；
[0195] K'＝4(2/π)0.5KIc，MPa·m0.5；
[0196] KIc为I型断裂韧性，MPa·m0.5；
[0197] E'为平面应变杨氏模量，E'＝E/(1‑v2)，MPa；
[0198] v‑岩石泊松比，无因次；
[0199] w为缝宽向量；
[0200] 式(6)可得到单元发生扩展的临界宽度，通过比较当前时刻尖端单元的宽度是否达到临界宽度，若达到则增加单元，否则单元数量不增加。
[0201] 通过步骤S30建立耦合了井筒流动模型的流固耦合方程式(9)，求解流固耦合方程可得到更新单元的宽度和压力分布以及井底压力；得到的井底压力与井筒模型进行对比，直至收敛。
[0202] 步骤S40：建立时间递增模型；
[0203] 具体地，所述建立时间递增模型包括：
[0204] 时间步计算公式为：
[0205]
[0206] 时间递增为：
[0207] tm＝t+Δt (12)
[0208] 其中，λ为松弛因子，λ∈(0,1]；
[0209] lf为裂缝长度，m，可以根据裂缝面积进行近似计算；
[0210] 为了避免裂缝长度太小时时间步过大，可以选择ε＝0.08；
[0211] 为了避免积分步数太大而可能出现数值震荡，smax＝301，最大积分步数为301；
[0212] 为所有待检查单元的最大扩展速度，m/s；
[0213] ΔtE为显式算法的时间步长，s；
[0214] Δx为x方向单元大小，m；
[0215] Δy为y方向单元大小，m。
[0216] 步骤S50：建立压裂缝内支撑剂运移计算模型；
[0217] 具体地，所述压裂缝内支撑剂运移计算模型的建立过程包括：
[0218] 考虑支撑剂后，缝内流动为携砂液流动。支撑剂视为携砂液的组分，裂缝内携砂液流动和支撑剂运移的质量守恒方程为
[0219]
[0220]
[0221]
[0222]
[0223] 其中，w为缝宽向量，m；
[0224] qs、qp分别为携砂液和支撑剂的缝内流量，m3/s；
[0225] q1为滤失项；
[0226] Qk为第k裂缝的入口流量，m3/s；
[0227] p为缝内流体压力；
[0228] 为归一化的加砂浓度，无因次；
[0229] 为归一化的支撑剂浓度，
[0230] φ为支撑剂体积浓度，无因次；
[0231] φm为最大容许浓度，无因次；
[0232] a为支撑剂半径，m；
[0233] g为重力加速度，m/s2；
[0234] μ为压裂液粘度，mPa.s；
[0235] ρp、ρf分别为支撑剂和压裂液密度，kg/m3；
[0236] ez为z轴方向单位向量；
[0237] 为描述携砂液流动的无因次函数；
[0238] 分别为描述支撑剂对流和沉降的无因次函数；
[0239] x为缝内任一点的坐标向量；
[0240] xin,k为第k裂缝的入口坐标向量；
[0241] δ(x‑xin,k)为Delta函数，输入值为0时取1，其他值为0；
[0242] B为阻塞函数，描述支撑剂在缝内砂堵现象。
[0243] 步骤S60：判断该时刻是否加入支撑剂，若是，则采用高阶加权本质不震荡(WENO，Weighted EssentiallyNon‑Oscillatory)有限差分方法对压裂裂缝扩展过程中支撑剂运移进行计算，得到缝内支撑剂分布；若不是，则进入下一步；
[0244] 所述步骤S60包括：
[0245] 步骤S601：判断该时刻是否加入了支撑剂。因为支撑剂运移方程为非线性双曲方程，其数值求解过程中会出现间断而导致震荡误差，为了解决这一问题并同时要满足特征线稳定性要求，加上本发明中模型为结构化网格，于是采用加权本质不震荡有限差分方法求解支撑剂运移方程。另外，本发明中支撑剂模型与裂缝扩展模型采用相同网格，而裂缝扩展模型的空间步长一般为数米级别，因此支撑剂运移的计算需高阶计算格式。
[0246] 具体的计算过程：调用式(13)与固体方程进行耦合得到裂缝扩展流固耦合方程，该方程通过显式高效方法求解。求解裂缝扩展模型，得到裂缝宽度、压力等结果后，缝内携砂液的流场即可确定；根据携砂液流场，求解支撑剂运移的连续性方程(14)，即可得到支撑剂分布。
[0247] 步骤S602：判断无支撑剂加入，则进入下一步。
[0248] 请参阅图2和图3，通过步骤S60判断该时刻是否加入支撑剂，若是，则调用步骤S50中的携砂液流动方程式(13)与固体方程进行耦合得到裂缝扩展流固耦合方程，该方程通过显式高效方法求解。求解裂缝扩展模型，得到裂缝宽度、压力等结果后，缝内携砂液的流场即可确定；根据携砂液流场，求解支撑剂运移的连续性方程(14)，即可得到支撑剂分布如图2。
[0249] 步骤S70：判断时间是否达到施工时间，若没有则返回步骤S4，否则结束计算。
[0250] 具体地，所述步骤S70包括：
[0251] 判断时间递增是否达到施工时间，若没有则返回步骤S40，否则结束计算。
[0252] 传统方法缝内支撑剂运移通常为静态裂缝，即将裂缝扩展模型与支撑剂运移模型割裂开来，并没有一体化同步模拟。由于支撑剂运移是非连续问题，本发明采用加权本质不震荡方法，可以实现对支撑剂浓度不连续的准确捕捉，实现支撑剂运移的准确求解。
[0253] 以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。
