夜雨聆风
PINN 求解孔隙弹性固结正问题：Terzaghi 一维流固耦合 PDE 的分步训练实践
本教程是我们AI4PDE&CFD暑期课程固体力学专题系列的一部分，与课程配套的 PINN 算法代码实验同步更新。本期围绕饱和多孔介质的一维固结问题，展示如何将 Biot 孔隙弹性控制方程嵌入神经网络损失函数，在 JAX/Equinox 框架下完成流固耦合偏微分方程的正问题求解。读者无需预先阅读同系列其他篇章，但若能具备弹性力学与渗流力学的基础，理解会更加顺畅。
相关代码教程：
https://github.com/xgxgnpu/Physics-informed-vibe-coding
摘要导读
Terzaghi 一维固结问题是孔隙弹性力学的经典基准算例：饱和土柱在顶面瞬时施加上覆荷载后，孔隙水压力自排水面向深处逐渐消散，土体骨架随之发生压缩沉降。该过程由压力扩散方程与准静态力平衡方程耦合控制，时间尺度长、空间梯度陡，对流固耦合 PINN 的训练策略与损失设计提出了较高要求。我们在厚度 的竖向土柱上，采用双网络 PINN 架构分别表征无量纲孔隙压力 与竖向位移 ，并实施"先流体后固体、冻结压力网络"的两阶段分步训练，在 GPU 上约 101 秒完成全部训练。
与代码内置的 Fourier 级数解析解对比，全域 相对误差为：压力 （0.62%），位移 （0.50%）。两网络合计可训练参数 10,162 个，配置点 20,000 个。以下表格汇总了本实验的核心指标，供快速查阅。
本教程依次介绍问题的物理背景与控制方程、PINN 网络维度与损失构造、分步训练策略的设计逻辑，以及基于我们实验数据的精度对比与结果分析。
一、引言
1.1 问题的物理背景
土体、储层岩石、生物组织等天然材料往往以多孔介质的形式存在：固相骨架构成孔隙网络，孔隙中充满流体。当外荷载作用于饱和多孔介质表面时，荷载最初主要由孔隙流体承担，随后流体自排水边界渗出，有效应力逐步转移至骨架，介质发生与时间相关的压缩变形。Terzaghi 于 1923 年建立的一维固结理论正是这一物理图像的数学表述，也是 Biot 三维孔隙弹性理论的简化起点。
从一维竖向土柱的理想化模型出发，我们可以将问题归结为在时空域 上求解耦合的压力场与位移场。底部为刚性不透水边界，顶部为排水面并承受恒定上覆荷载 。这一构型在实验室固结仪试验、软土地基沉降分析以及油气开采中的储层压实问题中反复出现，具有明确的工程对应关系。
1.2 正问题表述与 PINN 的适用性
本实验求解的是正问题：全部材料参数与边界条件已知，目标是预测任意 处的 与 。正问题在工程上对应于"给定土性参数与荷载历史，预报沉降与孔压时程"的前向分析需求。传统数值途径需要对时空域进行网格或差分离散，并单独处理流固耦合项的时间积分与混合边界格式。
PINN 以配置点采样替代空间网格，以自动微分替代手工差分格式，将 PDE 残差、初值条件与边界条件统一写入损失函数，由神经网络在训练过程中同时满足这些约束。对于 Terzaghi 问题，压力方程为抛物型扩散方程，位移方程为椭圆型平衡方程，两者通过 Biot 系数 耦合——这种"双方程、双物理场"的结构天然适合用两个独立子网络分别表征，再通过损失函数中的交叉导数项建立联系。我们在本实验中采用分步训练策略，先锁定压力场再求解位移场，以规避同时优化两个强耦合残差时常见的收敛困难。
二、方法原理与核心数学
2.1 有量纲控制方程
考虑竖向一维饱和多孔介质，设 为竖向坐标（向上为正）， 为时间。记孔隙流体压力为 ，竖向位移为 。在 Biot 孔隙弹性框架下，小变形、Darcy 渗流、流体可压缩性由 Biot 模量 表征，控制方程组为：
流体质量守恒（压力扩散）：
其中 为储水系数， 为排水体积模量， 为渗透率， 为流体动力粘度， 为 Biot 系数。
准静态力平衡（竖向）：
其中 为剪切模量， 为 Lamé 第一常数， 对应有效应力贡献。
有效应力原理 （压应力为正约定下）将上述两式联系起来：孔压消散改变有效应力，驱动骨架进一步压缩。
2.2 无量纲化
为减少参数量纲对神经网络优化的影响，我们引入特征尺度进行无量纲化。代码中采用的特征量为：
定义无量纲变量 ，，，，经代换后得到代码中实际求解的方程组：
方程 (1) 为标准热传导型扩散方程，方程 (2) 表明位移场的 Laplacian 由压力梯度驱动，即 Biot 耦合的数学体现。
2.3 边界条件与初值条件
顶部应力边界来自有效应力关系： 对应 的无量纲形式，减去 后等于 （压应力为负）。
2.4 解析解基准
精度评估采用代码内置的 Fourier 级数解析解（200 项截断）。压力场：
位移场由级数项 与 Terzaghi 解的解析积分项叠加得到，再除以 无量纲化。所有 PINN 预测值与解析解的对比数据均来自本实验输出的 predictions.npz，不包含任何外部文献数据。
三、PINN 架构与维度分析
3.1 双网络结构
我们定义 TerzaghiPINN 模块，包含两个结构相同但参数独立的 MLP 子网络：net_p 输出 ，net_v 输出 。每个 MLP 接受二维输入 ，经 4 层隐藏层（每层 40 神经元， 激活）后输出标量，权重采用 Xavier 均匀初始化。
维度变换过程可概括为：
net_p | ||||
net_v | ||||
| 10,162 |
参数量计算：首层 ，中间三层各 ，输出层 ，单侧合计 5,081。
选择 而非 ReLU 的原因在于 PDE 残差涉及二阶空间导数与一阶时间导数，要求激活函数在全域光滑可微；ReLU 的导数在零点不连续，不利于高阶自动微分的数值稳定性。
3.2 配置点采样
训练数据通过随机配置点生成，不含任何观测标签，全部约束来自物理方程与边界条件。
| 合计 | 20,000 |
无量纲时间域上界 ，对应物理时间约 15,000 s（约 4.2 小时），覆盖固结过程的主要阶段。
3.3 PDE 残差与损失函数
JAX 自动微分计算网络输出对输入的各阶导数，构造 PDE 残差。流体残差对应方程 (1)：
def_fluid_residual(model, y, t): pd_t = jax.grad(_pd, argnums=2)(model, y, t) # ∂p_d/∂t_d pd_yy = jax.grad(jax.grad(_pd, argnums=1), argnums=1)(model, y, t)return pd_t - pd_yy # 应趋于 0固体残差对应方程 (2)，其中 来自已训练（或冻结）的 net_p：
def_solid_residual(model, y, t): vd_yy = jax.grad(jax.grad(_vd, argnums=1), argnums=1)(model, y, t) pd_y = jax.grad(_pd, argnums=1)(model, y, t)return (1 + 2*NU_STAR) * vd_yy - BIOT_B * pd_y总损失由 PDE 残差均方值与边界/初值约束组成，IC 与 BC 项权重为 10，以强化硬约束的满足程度：
3.4 分步训练策略
流固耦合方程 (1)(2) 理论上可联立求解，但我们在实践中发现：若两个网络的参数同时更新，固体 PDE 残差中的 在压力场尚未收敛时会引入大量噪声梯度，导致位移网络训练不稳定。因此采用两阶段分步策略：
阶段一（流体）：仅优化 net_p，最小化 ，直至压力场在全域达到 量级的 相对误差。
阶段二（固体）：冻结 net_p 的全部梯度（对 net_p 参数返回零梯度），仅优化 net_v，最小化 。冻结机制如下：
grads_p = eqx.tree_at(lambda m: m.net_p, grads, jax.tree.map(jnp.zeros_like, eqx.filter(model.net_p, eqx.is_array)))
flowchart TB subgraph input [输入空间] yt["(y_d, t_d) in R2"] end subgraph phase1 [阶段一流体] netP["net_p: MLP 4x40"] pdeF["PDE: dp/dt - d2p/dy2 = 0"] end subgraph phase2 [阶段二固体] netV["net_v: MLP 4x40"] pdeS["PDE: d2v/dy2 - b dp/dy = 0"] freeze["net_p 冻结"] end yt --> netP --> pdeF netP --> freeze yt --> netV --> pdeS netP --> pdeS每个阶段内部采用"随机梯度预训练 + 拟牛顿精细收敛"的两段式优化：前段使用余弦衰减学习率（）进行 10,000 轮迭代，并施加全局梯度范数裁剪（阈值 1.0）；后段切换至拟牛顿法进行最多 3,000 步精细优化，以进一步压低 PDE 残差。
net_p | |||||
net_v | net_p |
值得注意的是，阶段一精细优化后压力 误差几乎不变（），说明随机梯度预训练已接近该网络架构下的误差下界；阶段二则仍有约 1.3% 的相对改善（），表明位移场的精细结构需要二阶优化才能充分捕获。
四、实验结果与分析
4.1 训练损失演化
流体阶段初始总损失为 18.77，其中初值条件损失 占主导——网络初始输出与 相差较大。经过约 2,000 轮迭代，总损失降至 量级，PDE 残差 稳定在 附近。固体阶段初始总损失 4.80，力平衡 PDE 残差 与顶部应力边界损失 共同贡献；至第 10,000 轮，总损失降至 ，PDE 残差仅为 。
4.2 相对误差演化
压力误差在阶段一从初始的 79.9% 单调下降，至第 10,000 轮达到 。阶段二启动后压力误差保持不变（net_p 已冻结），位移误差从阶段二第 1 轮的 97.0% 快速下降——第 500 轮即降至 6.75%，第 7,000 轮达到 ，最终精细优化至 。
4.3 时空剖面对比
上图在 共八个时刻，将 PINN 预测（虚线）与解析解（实线）沿深度方向对比。早期（）孔压梯度最为剧烈，PINN 在过渡区出现约 0.006 MPa 的局部误差；中期（）误差降至 0.0001–0.0005 MPa；晚期（）孔压接近全域消散，误差有所回升至 0.0006 MPa，可能与接近零解时的相对误差放大有关。
位移场在固结早期（）的剖面 误差偏大（最高约 35.3%），这是因为此时位移绝对值本身极小（量级 ），相对误差度量对微小偏差极为敏感。随固结推进，位移单调增大， 时剖面 误差已降至 ，最终时刻为 。
4.4 时空等值线分布
等值线图以 为横轴、物理深度 （m）为纵轴，清晰展示了固结的时空演化结构。压力场自顶部排水面（）向底部（）逐步消散，等值线呈典型的"扩散前沿"形态；PINN 预测与解析解在视觉上一致，误差主要集中在 的早期阶段和 的顶部过渡区，最大偏差约 。位移场随时间单调增大，等值线近乎水平，表明同一深度处的沉降速率在大部分时段内较为均匀 ；误差量级在 以下，分布均匀，未出现局部集中。
4.5 综合精度对比
| 合计 | 100.9 s |
上述结果表明，分步 PINN 在 Terzaghi 一维固结正问题上达到了亚 1% 的全域相对精度。压力场在阶段一即已收敛，位移场在阶段二经约 7,000 轮预训练后进入 误差平台，精细优化阶段进一步将其压低约 1.3%。
五、结论、常见问题与未来展望
5.1 结论
我们针对 Terzaghi 一维孔隙弹性固结正问题，构建了基于 JAX/Equinox 的双网络 PINN 求解器，将 Biot 耦合控制方程的无量纲形式直接嵌入损失函数，通过"先流体后固体"的分步训练策略，在 10,162 参数、20,000 配置点的设定下，于 100.9 秒内将压力与位移的全域 相对误差分别控制在 0.62% 与 0.50%。实验结果与代码内置 Fourier 级数解析解高度吻合，验证了 PINN 处理流固耦合抛物-椭圆型 PDE 系统的可行性。
分步训练的核心价值在于解耦两个物理场的优化竞争：压力扩散方程结构简单、收敛快，先行放入第二阶段后，net_p 已提供稳定的 驱动项，固体残差的优化曲面显著平滑。这一策略对更复杂的孔隙弹性问题（如二维渗流、非恒定荷载）具有直接的可迁移性。
5.2 常见问题与应对
5.3 未来展望
本实验仅涉及正问题求解，材料参数（、、 等）均已知。在暑期课程后续篇章中，我们计划将框架拓展至反问题——从稀疏的位移或孔压观测数据中辨识未知渗透率或 Biot 系数，届时只需将待求参数注册为可训练变量并加入数据损失项，无需改变网络架构。此外，二维平面应变固结、三相应力路径下的弹塑性孔隙介质、以及时变荷载 等方向，均可在此双网络 PINN 框架上逐步扩展。对于工程尺度的大时空域问题，配置点自适应加密与多尺度网络结构将是提升效率与精度的关键路径。
