裂缝性气藏分支水平井试井模型Python复现与参数敏感性分析

简介一份面向石油工程与试井分析专业人员的裂缝性气藏分支水平井试井模型技术资料系统覆盖从理论建模、Laplace变换求解、Stehfest数值反演到典型曲线绘制的完整流程。内容基于学术论文复现穿插可运行的Python代码及逐段解释涉及参数敏感性分析、九阶段渗流特征识别、分支几何参数对压力响应的影响规律等核心主题并给出分支参数优化、试井测试设计与数据解释流程的现场建议。研究明确指出分支长度越短或间距越大时第二拟径向流段更明显的结论有助于工程人员预判不同布井条件下的压力动态特征。完整成果整理在1个docx文档中约60KB以理论推导、代码片段和应用说明为主结构紧凑。已有73人学习使用适合具备石油工程或地质学基础、希望深入理解复杂气藏压力动态并掌握数值求解方法的开发工程师、科研人员及在校研究生。1. 裂缝性气藏分支水平井这套试井模型的 Python 复现能解决什么问题一口鱼骨型分支水平井的关井压力恢复数据摆在眼前用常规均质径向流模型硬拟合曲线挠度对不上九个流动阶段只认得出两三个。问题多半不在井而在模型天然裂缝性气藏本质是裂缝-基质双孔隙介质分支井又叠加了井储、表皮和分支干扰解析解绕不开 Laplace 变换和 Stehfest 数值反演。这套代码把论文里的数学模型完整落成了可运行的 Python 类从参数初始化、拉普拉斯空间解到典型曲线绘制与参数敏感性分析一条龙走通还附了流动阶段自动识别和分支长度、间距的影响分析。适合搞试井解释的现场工程师、做开发方案的研究人员以及正在复现论文的在校生——照着跑一遍就能自己动手改参数验证结论。2. 双孔隙介质与鱼骨型分支井建模八个储层参数如何进入无量纲方程2.1 双孔隙介质基本假设裂缝给通道基质给储量天然裂缝性气藏和均质砂岩最本质的区别是储层被天然裂缝切割成了两套系统裂缝系统渗透率高但孔隙体积小是主要的流动通道基质系统渗透率低但储集空间大是主要的储集体。基质里的气靠压差慢慢补给裂缝这个补给过程就是窜流。代码里的两个核心参数正是描述这套关系的。储容比 omega 是裂缝系统弹性储容占总储容的份额典型值在 0.010.1 之间代码里给的是 0.1偏大但便于看出窜流特征窜流系数 lamda 描述基质向裂缝补给的强度量级一般在 1e-5 到 1e-7 之间代码默认 1e-5对应窜流段出现在中后期。这两个参数加在一起决定了双孔隙典型曲线的“凹槽”位置和深度——凹槽越深说明基质补给越滞后储容比越小。参数代码变量典型量级物理含义分支长度L200800 m单支水平段长度分支间距d100400 m相邻分支垂直距离裂缝渗透率kf150 mD主流动通道渗透能力基质渗透率km0.010.5 mD基质供气能力储容比omega0.010.1裂缝储容占比窜流系数lamda1e-51e-7基质向裂缝补给强度井储系数C0.011 m³/MPa井筒存储效应分支数量n_branches36鱼骨分支数换到工程单位时要注意代码注释里井储系数写的是 m³/Pa实际现场习惯用 m³/MPa这两个单位差六个数量级直接决定早期井储段能不能拟合上。后文避坑章节会专门说这件事。2.2 控制方程与无量纲化为什么必须走 Laplace 变换双孔隙介质的渗流由两套方程耦合描述。裂缝系统方程是裂缝渗透率项加窜流源项等于裂缝储集项kf/μ · (∂²pf/∂x² ∂²pf/∂y² ∂²pf/∂z²) λ(pm − pf) φf·ct·∂pf/∂t基质系统方程则是负窜流项方向相反km/μ · (∂²pm/∂x² ∂²pm/∂y² ∂²pm/∂z²) − λ(pm − pf) φm·ct·∂pm/∂t这里 pf 和 pm 分别是裂缝、基质压力λ 是无量纲窜流系数φf、φm 是两套孔隙度ct 是总压缩系数。注意方程里已经做了拟稳态窜流假设即基质内部的压力梯度不显式求解只用一个窜流项描述基质和裂缝之间的交换。为什么这套论文不直接在时间域解偏微分方程因为鱼骨型分支井的边界条件太复杂——每个分支都有产量贡献分支之间还有压力干扰直接数值求解九段流动会非常麻烦。Laplace 变换把时间导数变成代数因子 s偏微分方程变成关于 s 的常微分方程边界条件在拉普拉斯域里处理起来简单得多解出显式表达式后再通过数值反演回到时间域。这就是整套代码的技术主线。2.3 离散单元法解析解的交叉验证手段论文的主解是拉普拉斯解析解但我一般会建议手里同时留一套离散单元法的验证代码用来检查解析解有没有量纲错。下面这段是原附录里的离散单元法实现把每支水平井分成若干段用影响系数构建线性方程组求解流量分布import numpy as np from scipy.sparse import lil_matrix from scipy.sparse.linalg import spsolve class DiscretizedModel: def __init__(self, n_segments20, n_branches4): self.n_segments n_segments self.n_branches n_branches self.total_elements n_branches * n_segments self.A lil_matrix((self.total_elements, self.total_elements)) self.b np.zeros(self.total_elements) def build_system(self, reservoir_params): 构建线性方程组系统reservoir_params 包含储层参数 for i in range(self.total_elements): # 自影响项本段对自身压降的贡献 self.A[i, i] -np.log( reservoir_params[rw] / reservoir_params[L_segment] ) / (2 * np.pi * reservoir_params[kf]) # 交叉影响项其他段对本段的压降贡献 for j in range(i 1, self.total_elements): distance self._calculate_distance(i, j, reservoir_params) self.A[i, j] 1 / (4 * np.pi * distance) self.A[j, i] self.A[i, j] # 右端项总产量均摊到每个离散单元 self.b[i] -reservoir_params[total_rate] / self.total_elements def solve_flow_distribution(self): 求解流量分布返回各段流量系数 A_csr self.A.tocsr() return spsolve(A_csr, self.b) def _calculate_distance(self, i, j, params): 计算离散单元间距同一分支用段长累加不同分支按几何关系近似 branch_i i // self.n_segments seg_i i % self.n_segments branch_j j // self.n_segments seg_j j % self.n_segments if branch_i branch_j: return params[L_segment] * abs(seg_i - seg_j) else: angle np.pi / 4 # 假设分支间夹角 45 度 return params[d] * np.sin(angle) params[L_segment] * (seg_i seg_j)逻辑上这段代码做的是把鱼骨分支井离散成若干源段每个段对井底的压降贡献分成自影响和交叉影响两部分自影响是段本身在井底产生的压降交叉影响是其他段通过地层传导过来的干扰。A 矩阵对角线是自影响非对角线是交叉影响。参数说明n_segments 控制离散精度20 段每支已经能较好逼近连续分支L_segment 是每段长度由分支长度除以段数得到夹角按 45 度近似如果现场知道实际分支夹角改这一行就行。求解出的流量分布如果和解析解的趋势对得上说明主模型的空间离散合理。3. 用 Python 复现 Laplace 解析解模型类设计、单位换算与 f(s) 的物理含义3.1 模型类初始化单位换算是第一个血泪坑复现试井模型最容易翻车的不是公式推导而是单位换算。原文代码里把这个问题处理得比较规范但初读时很容易忽略。先看初始化部分class FracturedGasReservoirModel: 裂缝性气藏分支水平井试井模型 L : 分支长度 (m) d : 分支间距 (m) kf : 裂缝渗透率 (mD) km : 基质渗透率 (mD) mu : 气体粘度 (cp) ct : 总压缩系数 (1/psi) rw : 井筒半径 (m) h : 储层厚度 (m) q : 产量 (m3/d) C : 井储系数 (m3/Pa) omega : 裂缝储容比 lamda : 窜流系数 n_branches : 分支数量 def __init__(self, L500, d200, kf10, km0.1, phi_f0.02, phi_m0.1, mu0.02, ct2e-5, rw0.1, h20, q1000, C0.1, omega0.1, lamda1e-5, n_branches4): self.L L self.d d self.kf kf * 9.86923e-16 # mD - m^2 self.km km * 9.86923e-16 # mD - m^2 self.phi_f phi_f self.phi_m phi_m self.mu mu * 0.001 # cp - Pa*s self.ct ct / 6894.76 # 1/psi - 1/Pa self.rw rw self.h h self.q q / 86400 # m3/d - m3/s self.C C self.omega omega self.lamda lamda self.n_branches n_branches四处换算逐一说明。渗透率乘 9.86923e-16 是把毫达西换成平方米这是渗流力学的标准换算因子粘度乘 0.001 是厘泊换帕秒压缩系数除以 6894.76 是每磅每平方英寸的倒数换成每帕产量除以 86400 是日产换秒产。漏掉任何一处典型曲线的横纵轴就会整体平移一到两个数量级拟合直接失败。特别提醒井储系数 C注释里写的是 m³/Pa但工程现场更常用 m³/MPa如果你现场拿到的 C 是 0.1 m³/MPa直接塞进这个类无量纲井储系数 CD 会差六个数量级早期段完全对不上。用之前先统一成 Pa 或 MPa 体系这是我的血泪经验。3.2 拉普拉斯空间的解逐项拆解核心公式核心方法是 laplace_space_solution它把论文里的解析解从纸面搬进了代码。逐项看def laplace_space_solution(self, s): 计算拉普拉斯空间的解s 为拉普拉斯变量返回无量纲压力 # 窜流项组合参数 eta eta self.lamda * self.omega * (1 - self.omega) / s # f 是双孔隙介质在拉普拉斯域内的等效修正项 f (self.omega * (1 - self.omega) * s self.lamda / (3 * eta * np.tanh(np.sqrt(eta)) - self.lamda 1)) # 裂缝系统径向流参数rD 是取点处的无量纲半径 alpha np.sqrt(s * self.phi_f * self.mu * self.ct / self.kf) rD self.rw * alpha K0 kv(0, rD) I0 iv(0, rD) # 分支井影响函数分支越多、越长井底压降越被摊薄 branch_effect self.n_branches * self.L / (self.d * np.sqrt(np.pi * s)) # 无量纲井储系数注意 h 和 rw 参与归一化 CD self.C / (2 * np.pi * self.phi_f * self.ct * self.h * self.rw**2) skin 0 # 假设表皮系数为 0 # 拉普拉斯空间压力解含井储和表皮项 pD ((K0 skin * s * I0) / (s * (I0 CD * s * (K0 skin * s * I0))) * branch_effect / (1 f)) return pD第一处是 eta它把 lamda、omega 和拉普拉斯变量 s 组合在一起是后面 f 函数的核心参数。第二处 f 是整个双孔隙模型的关键修正项实空间里它对应基质系统对裂缝系统压力响应的缓慢补给后面单独展开。第三处 alpha 和 rD 把井底取点换算成了修正贝塞尔函数 kv、iv 的自变量K0 和 I0 分别是第二类和第一类零阶修正贝塞尔函数径向流解的标准形态。第四处 branch_effect 是分支井的近似影响项物理含义是分支越多、单支越长同样的总产量下井底压降被多条分支分摊压力降落变缓。3.3 f(s) 里 tanh 的几何含义为什么窜流项长这样初次看3 * eta * np.tanh(np.sqrt(eta))会觉得奇怪为什么要用双曲正切。这其实是拟稳态窜流假设在拉普拉斯域里的解析形态。实空间里基质被假设为无限长条状或板状基质内部压力呈抛物线分布基质向裂缝的补给通量随时间演化把这个演化过程做拉普拉斯变换tanh 项就自然浮现。系数 3 是形状因子归一化的结果——不同基质几何形状板状、球状、圆柱状会得到不同的形状因子3 对应的是板状基质块。工程上不需要深挖 tanh 的数学推导但要知道它对曲线形状的影响eta 很小时tanh(√eta) ≈ √etaf 函数退化成与 λ 直接相关的简单形式对应早期裂缝系统径向流eta 很大时tanh 项趋近 1f 里的窜流修正逐渐消失对应总系统拟径向流。中间那段由 1 到接近极值的过渡就是典型曲线上那个下凹的窜流段。这就是为什么 ω 和 λ 的组合决定了窜流段出现的位置。4. Stehfest 数值反演与典型曲线N 取 8 的经验与九阶段标注方法4.1 Stehfest 反演系数计算、N 值玄学与精度边界拉普拉斯空间解是在 s 域里的要画典型曲线必须回到时间域。Stehfest 算法是试井分析里最常用的数值反演方法思路很简单在若干离散的 s 点s i·ln2/t上对拉普拉斯解取样乘上权重系数 V 求和再乘 ln2/t 得到实空间近似解。权重系数只和 N取偶数有关和具体函数无关。from math import factorial def _stehfest_coefficients(self, N): 计算 Stehfest 系数N 必须为偶数常用 8 或 10 V np.zeros(N) for i in range(1, N 1): k_min int((i 1) / 2) k_max min(i, N // 2) sum_k 0.0 for k in range(k_min, k_max 1): num k ** (N // 2) * factorial(2 * k) den ( factorial(N // 2 - k) * factorial(k) * factorial(k - 1) * factorial(i - k) * factorial(2 * k - i) ) sum_k num / den V[i - 1] (-1) ** (i N // 2) * sum_k return V def stehfest_inversion(self, t, N8): Stehfest 数值反演t 为实空间时间点返回无量纲压力 ln2 np.log(2.0) V self._stehfest_coefficients(N) sum_pD 0.0 for i in range(1, N 1): s i * ln2 / t pD self.laplace_space_solution(s) sum_pD V[i - 1] * pD return ln2 / t * sum_pD注意我把np.math.factorial换成了math.factorial新版 NumPy 已经把np.math移除了直接跑原附录代码会报错这是一处必踩的坑。系数计算里的双层循环看起来繁琐但本质就是把组合数形式的权重按定义式展开N 越大系数绝对值越大浮点误差也越大所以一般取 8 到 12论文代码默认 8 是个稳妥折中——这就是 N 值的玄学取大了高频振荡取小了精度不够。4.2 典型曲线生成与九阶段标注有了反演函数生成典型曲线就是在一个对数均匀的时间网格上逐点反演。这里有两个细节值得注意时间跨度用 logspace因为压力响应对对数时间才分辨得开计算点数量 100 个是速度和精度的平衡再少曲线就毛糙了。import matplotlib.pyplot as plt import numpy as np def generate_typecurve(self, t_end1e6, n_points100): 生成典型曲线时间轴对数均匀分布 t np.logspace(-2, np.log10(t_end), n_points) pD np.zeros_like(t) for i in range(len(t)): pD[i] self.stehfest_inversion(t[i]) return t, pD # 创建模型实例并生成典型曲线 model FracturedGasReservoirModel(L500, d200, n_branches4) t, pD model.generate_typecurve() # 绘制 loglog 典型曲线 plt.figure(figsize(10, 6)) plt.loglog(t, pD, b-, linewidth2) plt.xlabel(无量纲时间, fontsize12) plt.ylabel(无量纲压力, fontsize12) plt.title(裂缝性气藏分支水平井典型曲线, fontsize14) plt.grid(True, whichboth, ls--) # 按论文的九个流动阶段标注 stages [ (1e-3, 1e-2, 纯井储段), (1e-2, 1e-1, 过渡段), (1e-1, 1e0, 第一径向流), (1e0, 1e1, 线性流), (1e1, 1e2, 第二拟径向流), (1e2, 1e3, 分支干扰), (1e3, 1e4, 第三拟径向流), (1e4, 1e5, 窜流段), (1e5, 1e6, 总系统径向流), ] for t1, t2, label in stages: idx (t t1) (t t2) if any(idx): mid len(pD[idx]) // 2 plt.text(np.sqrt(t1 * t2), pD[idx][mid], label, hacenter, vacenter, bboxdict(facecolorwhite, alpha0.8)) plt.show()阶段标注用的是每个区间几何中心np.sqrt(t1 * t2)做锚点再取该区间压力数组的中间值定垂直位置这样标签不会和曲线重叠。从实际效果看九段识别可以辅助肉眼判断对数坐标里径向流段因为压力对时间呈对数关系loglog 曲线斜率接近 0线性流段斜率接近 1/2井储段斜率接近 1。这套几何特征比单纯看形状可靠得多。顺带说明一个原附录代码的坑generate_typecurve 在原始代码里缩进脱出了类直接运行会报missing 1 required positional argument: self把它缩进回类定义内就好。这也是读者复现时的高频报错点。5. 复现与应用的避坑指南五个容易翻车的参数和代码细节5.1 无量纲时间与实际压力对不上曲线平移一到两个数量级现象生成的典型曲线形状对但横轴无量纲时间和现场实际压恢时间换算后总是错位导数曲线也整体偏移。原因无量纲时间 tD 和压力 pD 的定义依赖 kf、μ、ct、rw、h 等一组参数只要其中一个单位的换算不一致比如 ct 用 psi⁻¹ 却直接套进 SI 公式整个横轴就会按指数平移。试井模型对单位极其敏感这是最常见的翻车点。解决我一般先在模型类里加一个自检方法把无量纲时间按定义反算回实际时间再用已知的均质径向流解析解做一次基准测试。若基准测试吻合再跑分支井模型。反算公式t_real omega * mu * ct * rw² / kf * tD每个参数先确认单位再代入。5.2 Stehfest 反演结果出现高频锯齿或负压现象典型曲线早期段振荡剧烈甚至出现负的无量纲压力loglog 图上像锯齿。原因Stenfest 的 N 取太大比如 16 以上阶乘项数值溢出或浮点抵消导致系数不精确另一个原因是时间 t 极小时s i·ln2/t 很大rD 偏离核心区贝塞尔函数计算失稳。解决先把 N 降到 8 或 10这是文献里的常用经验值再把时间范围从 1e-6 开始而非更小。若仍振荡检查 kv(0) 是否出现除零或 NaN必要时对 s 极小的点做数值截断。5.3 敏感性分析循环里没重置实例属性曲线互相污染现象跑不同分支长度的敏感性分析第二条曲线形状诡异和单变参数应有的趋势完全不符。原因原附录代码里用同一个 model 实例循环改 L但 d、omega 等属性保留着上一次循环的值。比如先跑 L800 再跑 L500d 可能还是上一轮的残留值曲线自然乱套。解决每轮循环重新创建实例或至少把参与计算的属性全部显式重置。我的习惯是写一个reset()方法把 L、d、omega、lamda 等全部恢复默认值再改目标参数避免属性残留污染敏感性结论。5.4 井储系数单位混淆m³/Pa 与 m³/MPa 差六个数量级现象早期纯井储段的斜率特征出不来CD 计算值异常大或异常小曲线整体形状与理论不符。原因注释里写 m³/Pa但现场报表常给 m³/MPa直接代入后 CD 差 10⁶ 倍无量纲井储系数完全失真。解决入口处统一规定单位体系。我建议把构造函数的 C 参数明确标注为 m³/PaSI现场拿到 m³/MPa 时先除以 1e6 再传入并加一行注释防止下次再踩。同理ct 的 1/psi 与 1/Pa 也要统一。5.5 两个代码级小坑np.math 被移除、方法缩进出类现象np.math.factorial报 AttributeError或调用生成曲线时报missing 1 required positional argument: self。原因新版 NumPy 移除了np.math命名空间需要改用 Python 标准库的math.factorial另外原附录代码里 generate_typecurve 定义时缩进脱出了类体导致它变成游离函数而不是绑定方法。解决把from math import factorial替换到代码顶部所有np.math.factorial改为factorial把 generate_typecurve 整体缩进回类内部。这两处改完代码基本能一次跑通。6. 把模型用起来参数敏感性、试井设计与拟合验证的技巧6.1 分支几何参数的影响方向L 与 d 的相反作用论文最关键的结论是分支长度越短或分支间距越长第二拟径向流段越明显。从物理上理解第二拟径向流对应裂缝系统整体的拟径向流动此时各分支的压降漏斗开始重叠。分支间距拉大后分支间的干扰减弱裂缝系统需要更长时间才进入整体径向流这个阶段在曲线上就会被拉长分支长度缩短同理压降波及范围小分支之间更早解耦。反过来分支越长或间距越短各分支几乎等效成一条长水平井第二拟径向流被吞掉曲线直接跳向第三拟径向流。做敏感性分析时L 取 200、500、800 三档d 取 100、200、400 三档就能清楚看到这个趋势。6.2 现场试井设计的重点事项用这套模型指导现场测试我最关注三个点。一是压恢时长要覆盖到总系统拟径向流按默认参数算无量纲时间要到 1e51e6 对应实际关井 110 天测试设计时应按这个量级安排关井周期。二是压力计的分辨率必须能识别早期井储段纯井储段持续的无量纲时间极短实际时间可能只有几秒到几分钟采样建议采用对数均匀布点而不是等间隔采样。三是拟合顺序先固定 ω 和 λ 确定窜流段位置再反演 L 和 d 匹配第二拟径向流的出现时间最后用 C 微调早期段。这样分步拟合比一次性全参数搜索稳定得多。6.3 验证结论是否可靠的两个手段第一用离散单元法交叉验证解析解的流量分布两者趋势一致才说明空间离散没问题。第二把生成曲线的导数曲线同时画出用径向流段斜率为零、线性流段斜率为 1/2 的几何特征做硬校验——形状像不代表对斜率对得上才是真的对。从那以后我每次拿到这类双孔隙加多分支的复现任务都强制先跑一遍单位换算自检、实例重置检查和斜率校验再谈画曲线和出结论。这几个习惯帮我省掉了大半返工时间希望帮到你。本文还有配套的精品资源点击获取