两相两点物质点法一维固结验证：饱和与非饱和响应

一维固结是检验流固耦合算法是否同时满足水力扩散和骨架变形规律的基础算例。只比较孔隙水压力并不足以证明模型正确：压力可能看似符合解析解，而沉降量仍会因为耦合系数、积分体积或基准应力处理不一致而偏离。

本文采用两相两点物质点法（two-phase two-point MPM），分别以饱和 L8 和非饱和 L08U 工况验证平均孔压固结度、沉降固结度及沿深度的压力剖面。新版 L08U 结果使用变形一致解析解，替代了早期仅适用于压力子问题的冻结孔隙率参考；因此，本文只将新版结果作为非饱和固结通过的依据。

一、物理问题与验证目标

模型为高$H=10\ \mathrm{m}$、截面$1\ \mathrm{m}\times1\ \mathrm{m}$的土柱，取$z=0$为不透水底面、$z=H$为排水顶面。侧面限制法向位移且不透水，底面限制竖向位移且不透水，顶面施加相对于初始平衡状态的$1\ \mathrm{kPa}$压缩荷载。

两相两点离散使用两套独立物质点：固相点携带骨架质量、应力和运动状态，液相点携带孔压、饱和度和液相状态。两相在背景网格上交换动量并求解水力—力学耦合。算例采用 10 层均匀网格，每层每相布置 8 个点，共 80 个固相点和 80 个液相点。

本文检查三个互相独立的量：

1. 平均超静孔压定义的固结度$U_p$；
2. 顶面沉降定义的固结度$U_s$；
3. 不同无量纲时间下的归一化超静孔压剖面$u/u_0$。

只有三者同时满足误差门限，才能认为压力扩散与变形耦合均得到验证。

二、控制方程与变形一致耦合

在忽略惯性的一维小应变固结极限下，两相模型可归结为骨架平衡、液相质量守恒和 Darcy 渗流。三维形式下的固相与液相动量方程可写为

$$(1-n)\rho_s\boldsymbol{a}_s =\nabla\cdot\boldsymbol{\sigma} +(1-n)\rho_s\boldsymbol{g} +\boldsymbol{f}^{\mathrm{drag}}, \tag{1}\label{eq:solid-momentum}$$

$$nS_r\rho_w\boldsymbol{a}_w =-\nabla(nS_rp_w) +nS_r\rho_w\boldsymbol{g} -\boldsymbol{f}^{\mathrm{drag}}, \tag{2}\label{eq:water-momentum}$$

其中，$n$为孔隙率，$S_r$为饱和度，$\rho_s$和$\rho_w$分别为固相与液相真实密度，$p_w$为孔隙水压力，$\boldsymbol{f}^{\mathrm{drag}}$为两相拖曳力。低速渗流时，拖曳关系退化为 Darcy 定律。

采用 Bishop 有效应力表示总应力：

$$\boldsymbol{\sigma} =\boldsymbol{\sigma}' -\alpha p_w\boldsymbol{I}, \tag{3}\label{eq:bishop-stress}$$

其中$\alpha$由饱和度和孔隙率决定。饱和极限下$S_r=1$、$\alpha=1$，式 \eqref{eq:bishop-stress} 退化为经典有效应力关系。

非饱和状态的变形一致线性化

早期 L08U 对照把土—水特征曲线视为只依赖压力，即冻结孔隙率后仅保留$\partial S_r/\partial p_w$。这种处理可以近似描述压力扩散，却遗漏了变形引起的孔隙率变化，因而不能作为沉降验证的解析参考。

新版参考同时保留

$$S_{r,p}=\left.\frac{\partial S_r}{\partial p_w}\right|_n, \qquad S_{r,n}=\left.\frac{\partial S_r}{\partial n}\right|_{p_w}. \tag{4}\label{eq:saturation-tangents}$$

由此得到压力耦合系数$\beta$、质量—体应变耦合系数$\gamma$和 Bishop 应力的体应变切线$\alpha_\varepsilon$：

$$\beta =\frac{\partial(\alpha p_w)}{\partial p_w} =\alpha+p_w(1-n)S_{r,p}, \tag{5}\label{eq:beta}$$

$$\gamma =1+\frac{n(1-n)}{S_r}S_{r,n}, \tag{6}\label{eq:gamma}$$

$$\alpha_\varepsilon =(1-n)\left[1-S_r+(1-n)S_{r,n}\right]. \tag{7}\label{eq:alpha-epsilon}$$

在基准吸力$p_0$附近，侧限有效模量相应修正为

$$M_{\mathrm{eff}}=M-p_0\alpha_\varepsilon. \tag{8}\label{eq:effective-modulus}$$

压力增量与体积应变增量必须在同一线性化状态下求解。相应的增量形式为

$$C_p\Delta p_w+\gamma\Delta\varepsilon_v +\Delta t\,\mathcal{H}(\Delta p_w)=R_p, \tag{9}\label{eq:mass-balance}$$

$$\Delta\boldsymbol{\sigma} =\mathbb{C}:\Delta\boldsymbol{\varepsilon} -\beta\Delta p_w\boldsymbol{I} -p_0\alpha_\varepsilon\Delta\varepsilon_v\boldsymbol{I}. \tag{10}\label{eq:stress-increment}$$

其中，$C_p$为储水系数，$\mathcal{H}$为离散渗流算子，$R_p$为质量平衡残差。饱和极限满足$S_{r,n}=0$、$\beta=\gamma=1$、$\alpha_\varepsilon=0$，可严格恢复经典饱和固结方程。

本算例采用小应变参考构形。固相与液相积分权重均固定在参考构形，避免基准吸力应力随当前体积被重复缩放。该处理是新版沉降曲线能够与解析解一致的关键物理契约，不应直接外推到大变形更新拉格朗日问题。

三、算例参数与解析基准

饱和 L8 与非饱和 L08U 保持几何、网格、弹性参数、初始孔隙率、饱和渗透系数和附加荷载一致。两者仅在初始水力状态及其相应的非饱和耦合系数上不同。

参数	符号	数值
土柱高度	$H$	$10\ \mathrm{m}$
弹性模量	$E$	$50\ \mathrm{MPa}$
泊松比	$\nu$	$0.3$
初始孔隙率	$n_0$	$0.3$
水的体积模量	$K_w$	$500\ \mathrm{MPa}$
饱和渗透系数	$k_s$	$1.0\times10^{-3}\ \mathrm{m/s}$
固相真实密度	$\rho_s$	$2600\ \mathrm{kg/m^3}$
液相真实密度	$\rho_w$	$1000\ \mathrm{kg/m^3}$
附加压缩荷载	$q$	$1\ \mathrm{kPa}$

饱和 L8

线弹性材料的侧限压缩模量为

$$M =\frac{E(1-\nu)}{(1+\nu)(1-2\nu)} =67.307692\ \mathrm{MPa}. \tag{11}\label{eq:constrained-modulus}$$

考虑水的可压缩性后，固结系数为

$$c_v =\frac{k_s} {\rho_wg\left(M^{-1}+n_0K_w^{-1}\right)} =6.594802\ \mathrm{m^2/s}. \tag{12}\label{eq:saturated-cv}$$

最终沉降量为

$$s_\infty=\frac{qH}{M} =1.485714\times10^{-4}\ \mathrm{m}. \tag{13}\label{eq:saturated-settlement}$$

无量纲时间因子定义为

$$T_v=\frac{c_vt}{H^2}. \tag{14}\label{eq:time-factor}$$

对于顶面排水、底面不透水的土柱，超静孔压和平均固结度采用 Terzaghi 单向排水 Fourier 级数作为解析基准。

非饱和 L08U

L08U 在$p_0=-400\ \mathrm{kPa}$的基准吸力状态激活同一土—水特征曲线。初始饱和度为$S_{r0}=0.931035$，相对渗透率为$k_{r0}=0.550073$。为了满足 Bishop 牵引平衡，加载后的初始孔压为$-398.799828\ \mathrm{kPa}$，对应初始超静孔压$u_0=1.200172\ \mathrm{kPa}$，并非直接复制饱和算例的$1\ \mathrm{kPa}$。

基准状态下

$$\beta_0=0.832637,\qquad \gamma_0=0.892770,\qquad \alpha_\varepsilon=-0.184672. \tag{15}\label{eq:unsaturated-couplings}$$

由式 \eqref{eq:effective-modulus} 得$M_{\mathrm{eff}}=67.233823\ \mathrm{MPa}$，变形一致固结系数为$c_u=0.377083\ \mathrm{m^2/s}$，最终沉降量为$1.487347\times10^{-4}\ \mathrm{m}$。非饱和工况的固结速度明显低于饱和工况，主要原因是基质吸力降低了相对渗透率并改变了储水与应变耦合。

四、饱和固结结果

饱和 L8 在$T_v=0.05$、$0.20$、$0.50$和$0.80$四个目标时刻的结果如下。$U_p$由平均超静孔压计算，$U_s=s/s_\infty$。

$T_v$	$t$（s）	$U_p$数值解	$U_p$解析解	$U_s$数值解	$U_s$解析解
0.05	0.7885	0.25637	0.25231	0.25642	0.25231
0.20	3.0327	0.50322	0.50409	0.50321	0.50409
0.50	7.5817	0.76285	0.76395	0.76281	0.76395
0.80	12.1308	0.88645	0.88740	0.88648	0.88740

除最早时刻受 10 层离散和输出采样影响较明显外，孔压与沉降两条曲线几乎重合。四个目标时刻的归一化孔压剖面均保持单调，并满足顶面排水、底面零通量的边界特征。

五、非饱和固结的新验证结果

旧 L08U 结果只使用冻结孔隙率的压力参考，且液相积分体积会随密度和饱和度变化而错误缩放基准吸力应力。因此，旧结果即使压力趋势合理，也不能证明沉降通过。新版结果同时修正了变形一致解析基准和小应变参考构形积分权重，未采用压力裁剪、人工阻尼或沉降回写。

新版 10 层算例在时间步尺度系数为$0.2$时得到：

$T_v$	$U_p$数值解	$U_p$解析解	$U_s$数值解	$U_s$解析解	压力剖面 RMS 误差	压力剖面最大误差
0.05	0.25310	0.25183	0.25188	0.25183	0.00344	0.00624
0.20	0.50553	0.50412	0.50485	0.50412	0.00210	0.00270
0.50	0.76535	0.76399	0.76525	0.76399	0.00180	0.00259
0.80	0.88833	0.88743	0.88828	0.88743	0.00113	0.00160

四个目标时刻的最大$U_p$绝对误差为$0.00141$，最大$U_s$绝对误差为$0.00126$；压力剖面的 RMS 与最大误差均低于$0.01$和$0.025$的门限。时间步尺度从$0.2$依次减小至$0.1$、$0.05$，以及网格从 10 层加密至 20 层时，结果均满足预设收敛门限。

下图上排分别比较饱和 L8 和非饱和 L08U 的$U_p$、$U_s$时间历程，下排比较四个目标时刻的压力剖面。空心点为 MPM 结果，连续线为各自的解析解。非饱和算例中，孔压与沉降均与变形一致解析解重合，说明新版结果已经从“压力子问题可用”提升为完整的一维非饱和固结验证通过。

六、数值稳定性与适用范围

1. 压力稳定化：低阶等阶插值在近不可压缩条件下可能产生体积自锁和棋盘状孔压。算例使用压力投影稳定化与质量加权压力恢复，并检查同一轴向平面上的粒子离散度和交替重构模态。
2. 时间与空间收敛：单个时间步或单一网格上的拟合不足以证明方法可靠。至少应同时比较时间步加密、网格加密、平均固结度和完整压力剖面。
3. 初始平衡：非饱和工况的附加荷载必须建立在基准吸力应力平衡之上。直接把饱和算例的初始超孔压复制到非饱和状态，会破坏初始牵引平衡。
4. 参考构形：本文结论适用于标准 MPM、线弹性、小应变和一维单向排水。大变形问题需要完整的更新构形线性化与独立验证，不能直接沿用本文的固定参考积分权重。
5. 验证边界：本文未验证塑性、GIMP、三维排水、受控通量入渗或显式两惯性声学响应。上述能力应分别使用与其物理机制相匹配的基准算例。

综合孔压、沉降、压力剖面以及时间和网格收敛结果，可以确认两相两点 MPM 在本文限定范围内同时通过饱和与非饱和一维固结验证。

本文含有 AI 生成内容