AI 搭的流固耦合：PFC × Python 管域水力压裂，从踩坑到沉淀成一个技能
先说结论：我用 AI 结对编程，把一套"PFC2D 颗粒力学 + Python 管域流体"的水力压裂双向耦合从零搭到能出论文级结果——11843 颗粒、400 次交换、60 分钟跑完，破裂压力 28.04 MPa（Hubbert–Willis 理论 24 MPa），质量守恒残差 3.6×10⁻⁹。更关键的是，整个工作流最后沉淀成了一个可复用的 AI 技能（skill），下次再做水力压裂算例，一句话就能调起来。
这篇文章讲三件事：功能怎么实现的、两个数值坑怎么填的、后处理怎么看，以及"沉淀成技能"这一步为什么值得单独说。
算例速览
几何：88×88 mm 试样，中心 r=2.5 mm 圆孔井筒
颗粒：11843 颗（rd 0.35–0.50 mm），平行胶结 29268 键
应力：σ_v 8 MPa / σ_h 4 MPa（裂缝走垂向）
注液：q=1e-7 m²/s，恒速注入到破裂
1. 一、为什么是"管域模型"
PFC 原生 CFD 模块只有 PFC3D 有、而且不含流场求解器（只交换数据），2D 走不通。水力压裂关心的是孔隙压力场驱动裂缝张开这个准静态过程，要的是"每个孔隙单元一个压力、每条接触一根管"的粗粒化网络——这正是 Cundall 管域模型（domain/pipe network）的主场：
域（domain） = 颗粒接触图（平面图）的一个面。中心圆孔是一个大多边形，颗粒之间的三角形孔隙是一个个小多边形——流体场的基本单元就是它们 管（pipe） = 一个接触，平行板泊肃叶：q = (a³/12μL)(p_j − p_i) 开度律：完好胶结 a = a₀（0.1 μm）；断键后 a = a_res + 张开量（跳到几 μm~几百 μm）
域体积 = 多边形面积（鞋带公式）− 各颗粒扇形面积。建网后拓扑冻结，每步只更新几何与开度。这套东西 PFC 里没有，得自己写——这就是 Python 的活。
难点在时间推进。裂纹一张开，开度涨三个量级，电导率涨九个量级，显式格式的稳定步长掉到纳秒级，直接不可用。所以压力场用全隐式后向欧拉：每个交换解一次 17427 阶稀疏线性方程组（scipy 直接法），步长无限制。
2. 二、耦合架构：FISH socket 交错，一个交换 = 一次松弛 + 一次流体步
机械和流体各干各的，靠 TCP socket 握手：
压力反力按几何精确积分：域 d 作用在顶点颗粒 i 上的净力 F = −p·r_i·(sinβ−sinα, cosα−cosβ)，每个域的力之和恒为零（闭合边界法向积分抵消，单测验证过）。
时间账：60 分钟里 97% 花在 PFC 机械求解（400 cycle × 9 s），Python 流体全套只占 0.25 s/交换——隐式解 0.16 s + 几何/力/socket 0.09 s。想提速找 PFC 多核或降 M，优化 Python 端没意义。
3. 三、两个数值坑，以及 AI 是怎么帮我定位的
3.1 坑 1：井筒压力 ±4 MPa 方波锯齿
第一版跑出来，p_well 严格交替 +3.88 / −3.80 MPa，锯齿幅值精确等于 K_f·rel_cap/2 = 2 GPa × 0.002 / 2。这不是随机噪声，是周期-2 极限环。
根因是交错耦合的固有滞后：pⁿ 产生的流体压力，要等 PFC 松弛 400 cycle 之后才反映到几何上，所以体积增量 dV 相对压力恒滞后一步，dV_k = χ·(p_{k−1} − p_{k−2})。只把 dV 当显式源项时，回路增益 G = K_f·χ/V 实测 1.3 > 1，必然振荡，被限幅钳成方波。
修法是把机械柔度 χ 计入流体电容：V_cap = V + K_f·χ，增益降为 G/(1+G) < 1，振荡按构造几何衰减。χ 用 dV_k/dp_{k−1} 在线估计 + 低通滤波——分子分母含同一个滞后压力增量，与阻尼强度无关，不会自激。修复后 dp 变号次数从 91/99 步降到 4/99 步。
3.2 坑 2：外边界不透水，全试样憋压
第一版把 222 根边界管直接禁用（无排水），结果整个试样压力场均匀抬升、没有远场漏失，破裂压力虚高 10%。
修法是加一个远场储层域：体积取试样孔隙体积的 100 倍、压力恒 0，边界管全部接到它上面，等价于常压排水边界。妙处在于它仍是"普通域"，漏失量自动进守恒账，一行记账代码都不用改。
一个诚实的细节：χ 修复的灵感来自我对数据的"较真"——锯齿幅值恰好等于理论值 K_f·rel_cap/2，说明是限幅饱和而非噪声。AI 先给了"循环内记账间隙随断键增长"的归因，被相关系数证伪（corr(Δrem,Δnbrk)=0.50 < corr(Δrem,KE)=0.86）后重做回归才挖出真因。人机配合里，人负责不信"听起来合理"的故事。
4. 四、结果：破裂压力 28 MPa，裂缝走垂向，账能对上
图 1：初始状态。左：11843 颗粒密实堆积，中心红圈是墙环法压出的平整圆孔（r=2.5 mm）；右：29268 条胶结网络，0 损伤。
造孔手法参考了预制裂隙技能的核心思路：先插 36 段圆环墙 → 平衡压平孔壁 → 删环内颗粒 → 加胶结（墙撑着孔不被焊死）→ 伺服 → 删墙。直接删球再胶结会留颗粒级锯齿孔壁，还会扯断跨孔键产生假损伤。
图 2：井筒压力曲线。线性升压 6.5 s → 冲顶 28.04 MPa 破裂 → 跌到 15 MPa → 稳态渗流。原始曲线与配对均值完全重合 = 锯齿已消除。
破裂压力怎么读：取全曲线峰值，不是首次破键。第 84 次交换井壁就掉过 1 根键（p=13.6 MPa），那只是局部起裂；真正的宏观破裂是第 215 次交换压力冲顶、张开管雪崩那一下。Hubbert–Willis 理论值 (3σ_h − σ_v + T₀) = 24 MPa，实测高 17%——纯圆孔没有预制缺陷的尖端应力集中，起裂更难，方向符合预期。
图 3：孔隙压力场，流体场基本单元直接可视化。
四联图：(a) 全部 17426 个域多边形填色 + 轮廓，井筒红框；(b) 压力场以圆孔为中心向外衰减；© 井壁放大（±8 mm），看圆孔与三角孔隙共边拼接的几何关系；(d) 38 根张开管 = 裂缝网络。注意一个坑：拓扑冻结的多边形必须配建网时刻坐标画，配末态坐标孔壁会自交成星形（孔被压裂撑大了几毫米，旧环跟不上）。
图 4：流体体积账。注入 = 留存 + 漏失，末态去向：压缩 65%、膨胀 35.5%、漏失 0%。
这张图是这次放大算例最有信息量的后处理：裂缝垂向扩展到 ±28 mm 就停了，没打到 ±44 mm 边界，所以漏失几乎为零、注入流体全留在试样里。对比之前 40 mm 小模型（裂缝直接贯通边界，漏失 47%）——同一套代码，尺寸一放大，边界效应自己就显现出来了。守恒残差 3.6×10⁻⁹，账目机器精度闭合。
图 5：裂缝演化。破裂瞬间（6.5 s）断键与张开管雪崩式增长，之后随压力回落进入阶梯式慢扩展。
图 6：压力剖面。过井筒做水平/垂直两条剖线：垂向（裂缝向）压力延伸到 ±28 mm 还有值，水平向 2 mm 外就快速衰减——裂缝就是各向异性渗流通道，这张图把它定量画出来了。
图 7：开度场。左：对数直方图双峰，胶结峰 2.9 万根 @0.1 μm、张开峰 38 根 @10~160 μm；右：张开管空间分布，清楚看到从井筒上下端起裂的垂向裂缝带。
七张图 + KPI 卡片（破裂压力/偏差/断键数/漏失占比/守恒残差/墙钟）自动汇总成一个自包含 HTML 文件（base64 内嵌，无外部依赖，双击即开、直接转发），这是后处理的最终交付物。
5. 五、沉淀成技能：从"跑通一次"到"一句话复现"
算例跑完不是终点。我把整条链打包成了一个 AI 技能 pfc-hf-pipe：
SKILL.md 里最值钱的是两类内容：
踩坑清单（全是实测教训）——FISH 变量名大小写不敏感（hf_M 和 hf_m 是同一个变量，撞过一次名）；ball-facet 接触绝不能在插墙期间软化（颗粒会顶穿盒墙）；万级球伺服步数要加倍；def 体内调函数不带 @……这些坑一个能让调试浪费半天。
验收门 + 回归基准——每次跑完核对五条：守恒残差 <1e-6、升压段 dp 同号爬升（出现 ±K_f·rel_cap/2 方波 = χ 被关）、裂缝垂直最小主应力、破裂压力与 Hubbert–Willis 同量级、孔壁顶点数 ≈ 墙环接触数。再对一组基准数字：88 mm 模型 28.04 MPa / 44 断键 / 漏失 0，40 mm 模型 23.65 MPa / 漏失 47%。下次任何人（或任何 AI）改坏了什么，跑一遍就知道。
为什么值得单独做这一步：代码写完放在算例目录里，三个月后没人记得哪些参数动过、哪些是坑。技能目录 = 代码 + 说明书 + 验收标准 + 回归基准四件套，AI 下次接手时读 SKILL.md 就能直接干活，不用重新"考古"。
6. 六、局限与下一步
准静态交错耦合：一个交换 = M 个机械 cycle 松弛 + 一个流体步，流体时间与物理时间解耦，适合升压-破裂，不适合研究瞬态渗流速率 管域拓扑冻结：新裂缝沿旧三角网格走，网格分辨率（粒径）限制裂缝路径精度 2D only；3D 管域要换四面体面遍历，工作量不小 下一步：延长注入时间看裂缝能否贯通边界进入长期漏失；PFC 多核 + 降 M 提速；以及把"破裂前/冲顶/贯通"三时刻压力场做成动画
参考资料：Cundall (1994) 管域模型原始思想；Hubbert & Willis (1957) 破裂压力准则；Hazzard & Young 的 DEM 水力压裂系列工作。
