基于有限元法与密度相关NorSand临界状态本构土模型的板沉陷及原位bevameter剪切试验数值模拟

《Journal of Terramechanics》：Numerical modeling of plate sinkage and in-situ bevameter shear tests using FEM and the density-dependent NorSand critical state constitutive soil model

　　

•包含弯曲元件、三轴和 Bevameter 测试的开放获取多尺度数据集。
•创新性地将临界状态 NorSand 有限元模型应用于模拟现场 Bevameter 测试。
•实现了板沉陷和剪切的仅向前实验室到现场的预测，误差为 6%–16%。
•开发了用于 NorSand 硬化率标定的自动化贝叶斯优化方案。
•构建了土壤剪切刚度和硬化参数的密度相关元模型。

引言
Bekker-Wong 模型是地形力学中预测软土上车辆牵引力的既定框架和基石。该模型依赖于从现场板沉陷和现场 Bevameter 剪切测试中得出的经验土壤参数（例如 kc, k?, n, c, ?, K）。然而，基于 Bevameter 的土壤参数估算方法因实验测试装置因素而存在变异性（He et al., 2020; Apfelbeck et al., 2011）。从 Bevameter 现场剪切测试中估算内聚力（c）和摩擦角（?）等剪切强度参数尤为困难。例如，研究表明，对于同一土壤，Mohr-Coulomb 总应力摩擦角根据剪切机制、工具形状和估算方法的不同，范围从 17°到 38°，这对基于 Bevameter 的土壤表征以及依赖 Bevameter 测试的 Bekker-Wong 土轮模型的可靠性提出了挑战（Kruger et al., 2023）。同样，板沉陷 Bevameter 测试拟合的经验模型并不代表真正独特的固有土壤参数，而是反映了系统响应变量（He et al., 2019）。此外，影响现场剪切测试的重要因素，如推土力和侧向履带阻力，在实验中仍难以量化和测量（Kruger et al., 2023）。经验板沉陷和剪切测试模型也不适合推广到新条件；土壤状态、含水量、工具几何形状或边界条件的任何改变都需要重新标定参数。将垂直平板应力直接映射到圆形或柔性轮几何形状上也存在显著的认知模型不确定性，尤其是对于较小的车轮，这进一步阻碍了使用经验 Bekker-Wong 方法进行的预测努力（Meirion-Griffith and Spenko, 2011）。这些问题对行星漫游车、农业机械和军用车辆机动性应用具有广泛的影响。

基于物理的土壤本构模型，特别是基于临界状态土力学（CSSM）原则的模型，代表了在大地力学应用中取代经验模型的有前景的途径。这些模型提供了对材料宏观行为的见解，而非系统级响应，并依赖于真正不变的土壤参数，如临界状态摩擦比，这些参数不是几何形状、应力路径或土壤状态的函数（Wood, 1991）。为了模拟一般载荷、变化土壤状态和耦合水力学效应下的土壤行为，需要控制宏观土壤行为的本构模型（Karstunen and Amavasai, 2017）。具体而言，如有效应力框架，能够明确解耦土骨架的力学响应和孔隙流体压力，是捕捉不同排水条件下真实土壤行为的基本要求（Schofield and Wroth, 1968）。然而，将这些新兴的基于物理的本构模型直接应用于完整的轮胎 - 土壤相互作用引入了额外的复杂性，因为这涉及额外的子系统，如柔性轮胎（Sandu et al., 2019）。因此，建议首先针对 Bevameter 现场板沉陷和剪切测试验证这些方法，然后再扩展到更复杂的系统。

已有几项研究利用基于物理的土壤模型模拟了板沉陷测试，但没有已发表的工作模拟了 Bevameter 现场平移或旋转剪切测试。基于物理的板沉陷模拟利用了各种数值方法，包括有限元方法（FEM）、离散元方法（DEM），以及最近的光滑粒子流体动力学（SPH）（Hemmat et al., 2012, Xia et al., 2014, Bahrami et al., 2020）。
DEM 是一种广泛用于土壤建模的微观力学方法，擅长捕捉大变形、粒子相互作用和不连续行为。然而，其显式接触公式，加上需要非常小的时间步长和大量的粒子数，使其计算量很大（O'Sullivan, 2011）。DEM 从涌现的粒子相互作用中推导材料行为；至关重要的是，这种方法并不严格强制宏观屈服面，也不固有地嵌入 CSSM 的核心公理，如状态相关膨胀或明确规定的渐近收敛到临界状态（Chang and Yin, 2010）。然而，众所周知，大多数土壤，包括砂土、粉土和未固结粘土，当充分剪切时，倾向于在应力和体积空间中趋向唯一的渐近临界状态（Schofield and Wroth, 1968, Wood, 1991）。因此，DEM 难以重现某些宏观响应，如损伤、峰后软化（分叉）和压力相关的体积应变演化（Kuhn and Chang, 2006, Xiong and Wang, 2023, Aboul Hosn et al., 2017）。准确追踪孔隙比等宏观变量仍然具有挑战性；像粒子升标和简化几何形状这样的常见做法往往会扭曲孔隙比和结构演化，导致非物理膨胀（Zeraati-Shamsabadi and Sadrekarimi, 2025）。基于干燥接触力学的经典 DEM 缺乏固有的孔隙压力公式，这是真实土壤行为的关键方面，使其不适合模拟水力学耦合的饱和或非饱和土壤，除非与 DEM-CFD 框架耦合（Liu et al., 2015）。通过采用明确嵌入并标定这些 CSSM 公理和宏观土壤响应的方法来解决这一差距，是本研究的核心目标。

SPH 和物理信息神经网络（PINNs）是大地力学和地形力学中有前景的新兴方法。SPH 的无网格特性使其擅长处理大变形；然而，它缺乏强制严格屈服面的稳健公式，难以遵循宏观应力 - 应变路径，并且不强制临界状态原则，如膨胀或收缩（Islam et al., 2020）。虽然 PINNs 在求解一般偏微分方程约束问题上很强大，但在大地力学和岩土工程建模中仍处于早期开发阶段，实验研究仅模拟了基本土壤行为（Lu and Mei, 2022）。PINNs 在 Bevameter 测试中的应用尚未被探索。尽管 PINNs 可以在稀疏甚至没有数据的情况下运行，但边界或载荷条件的每次变化都需要重新训练神经网络（Wang et al., 2025）。因此，PINNs 表现出有限的边界条件泛化能力，并产生高昂的计算训练成本（Wong et al., 2025）。由于这些限制，本研究决定不采用上述基于物理的方法。

在基于连续介质的 FEM 方法中，存在各种本构公式。Mohr-Coulomb 模型因简单而广受欢迎，但它假设线性弹性后接完全塑性，忽略了非线性应变硬化/软化行为（Yu, 2007）。Mohr-Coulomb 参数也不是独特的固有土壤基本属性（Craig, 2004）。同样，Drucker-Prager 模型也是一个简单的弹塑性模型。帽型硬化变体在压缩区域引入了一些非线性；然而，它仍然省略了临界状态框架，使用固定的屈服面而没有膨胀或土壤密度依赖性，并且需要针对不同土壤状态使用不同参数（Sinha et al., 2010）。Cam-Clay 和修正 Cam-Clay 模型结合了更先进的 CSSM 概念，更好地代表了体积变化对剪切强度的影响。然而，它们假设孔隙比 e 和平均有效应力 p'之间存在唯一关系，这对于砂土是不成立的（Shuttle and Jefferies, 2010）。砂土根据初始密度表现出无限条正常固结线（NCLs），这种行为无法被 Cam-Clay 等固定屈服面模型捕捉。历史上，从实验室标定数据进行的现场尺度预测不可靠（Shoop, 1993, Sion et al., 2022）。这被假设为源于简化的、非状态相关的模型，这些模型不遵守 CSSM 公理。

为了克服这些限制，本研究采用在有效应力框架下制定的 NorSand 本构模型，明确分离骨架应力和孔隙压力效应，并使用大变形 FEM 求解器实现（Jefferies, 1993）。NorSand 通过引入状态参数 ψ=e?ec 扩展并泛化了 Cam-Clay，该参数表示 CSSM 框架中当前孔隙比与临界状态孔隙比之间的距离（Jefferies, 1993, Shuttle and Jefferies, 2010）。这种状态相关公式允许 NorSand 的屈服面和剪切强度随加载而演化，捕捉收缩（超临界，ψ>0）和膨胀（亚临界，ψ<0）行为。动态屈服面由密度相关的动员摩擦比 ηi 控制，该比值适应当前土壤密度状态，不同于传统模型中固定的临界状态摩擦比 Mtc（Shuttle and Jefferies, 2010）。这使得单一标定的本构模型能够表示从松散、易液化状态到致密、膨胀状态的土壤响应。内部而言，NorSand 明确建模膨胀并追踪孔隙比，这是许多传统本构模型中不存在但对于砂土至关重要的能力。

所提出的方法基于 CSSM 理论，采用实验室标定的 NorSand 模型及其不变土壤参数，旨在提供一种机制的、与边界条件无关的预测，预测现场板沉陷和现场剪切响应，涵盖不同的土壤密度和动员摩擦比，并具有固有的体积收缩和膨胀追踪与建模能力。这与上述讨论的传统非状态相关或固定屈服面本构模型（Mohr-Coulomb, Drucker-Prager 和 Cam-Clay），以及缺乏严格 CSSM 强制执行的 DEM/SPH 方法形成对比。本工作通过以下关键创新贡献推进了地形力学的发展：
1. 一个关于 Cullinan 砂的综合多尺度实验数据集，包括 6 个小应变弯曲元件剪切波测量、18 个宏观实验室三轴（CU/CD）测试，以及使用专门建造的 4 自由度 Bevameter 在五种土壤密度状态下进行的 52 个现场尺度现场板沉陷和剪切测试。一项针对 30 篇 Bevameter 相关出版物（编纂自 Kruger (2022) 的参考文献列表和 Wong (2026) 的综述）的调查表明，原始力 - 位移 Bevameter 数据很少与实验室测试数据一起发表，并且未识别出具有稳定存档链接的公开可访问数据集。该数据集弥补了这一差距，提供了用于多尺度标定和验证的配对现场和实验室数据。
2. 将密度相关的临界状态本构土壤模型（NorSand）创新性地应用于现场 Bevameter 工具 - 土壤相互作用，建立了一个基于 CSSM 的框架，作为 Bekker 方法中高度可变、几何依赖的经验系数的替代方案，从系统级响应转向与几何形状和边界条件无关的本构建模。与先前基于物理的地形力学方法不同，它严格强制宏观状态和密度相关的 CSSM 公理，包括膨胀、动态屈服面和渐近收敛到临界状态。
3. 对 Bevameter 板沉陷和现场剪切测试进行仅向前的实验室到现场预测建模和验证，通过有限元求解的前向本构土壤模型（NorSand）实现了 6%–16% 的现场板沉陷和剪切响应相对误差，仅通过实验室数据标定，并针对现场测量值（未见过数据）进行验证，没有任何事后调整或反向拟合。这与典型的 Bevameter 实践形成对比，后者是从现场测试中反向拟合土壤参数。
4. 引入了一种新的 NorSand 硬化率自动化标定程序，使用贝叶斯优化和 RBF 代理模型，整合了数据驱动的循环 FEM 模拟与快速、近最优、具有在线超参数优化的 RBF 代理模型，取代了 NorSand 标定中通常使用的“手动”调整标定。
5. 通过引入新的经验关系，将小应变剪切模量 G0、大应变剪切模量 G 和塑性硬化率 H 与可测量的土壤量如孔隙比 e 和平均有效应力 p'联系起来，为 NorSand 增强了剪切刚度和硬化率的密度和应力相关经验元模型，扩展了 Shuttle 和 Jefferies (2010) 的元模型（模型的层次模型）概念。这些先前在 NorSand 中被视为固定的参数，被制定为状态相关关系，增强了预测准确性。至关重要的是，这些额外关系使得能够直接从可测量的现场土壤属性中估算现场剪切刚度和硬化率。

包含实验室和现场测试的多尺度 Cullinan 砂数据集以开放获取形式提供：https://doi.org/10.5281/zenodo.17584462。这旨在支持广泛大地力学建模方法的开放研究和发展，包括数据驱动和基于物理的方法，如 FEM、DEM、SPH、PINNs 和经典经验模型。此外，NorSand 标定的贝叶斯优化库已开源发布 https://github.com/raykruger/BayesNorSand。

NorSand 本构土壤模型由一组耦合的非线性代数方程和普通微分方程组成，控制本构应力 - 应变关系和塑性流动运动学；这些需要数值积分求解，因为不存在闭合表达式（Jefferies and Been, 2016）。本研究使用 NorSand-M 版本 6（Shuttle and Jefferies, 2010），基于 Jefferies (1993) 的原始公式构建。控制本构方程总结在表 1 中，该表定义了所有相关量，包括不变材料参数、演化状态变量和积分期间使用的辅助隐藏状态变量。
NorSand 依赖于八个规范的土壤参数（G, ν, Γ, λ, Mtc, N, χ, H）（Shuttle and Jefferies, 2010）。所有参数遵循标准的 CSSM 术语；有关详细阐述，请参阅 Jefferies (1993) 和 Shuttle and Jefferies (2010)。

在本文的工作中，公式被扩展了七个额外的数据驱动元模型土壤参数：两个用于硬化率H，三个用于小应变弹性G0，两个用于大应变剪切衰减G/G0，总共得到十五个需要通过实验室内标定的土壤参数。每个参数都有明确的物理意义，共同构成了一组独特的土壤参数，能够捕捉土壤力学行为的不同方面（Shuttle和Jefferies，2010）。为了完整捕捉应力路径和加载历史，公式（1）中的数值实现追踪了总共45个（38个标准参数＋7个额外参数）耦合的土壤状态量，每个有限元积分点处都有。除了固定的材料参数外，这些包括：（i）连续介质应力和塑性应变张量；（ii）可观测的演化状态变量（如密度状态ψ、平均有效应力p′和动员摩擦比η）；（iii）隐藏的内部状态变量（如镜像压力pi′、镜像摩擦比Mi、镜像状态参数ψi和剪胀性Dp），这些变量不可直接观测，但对于追踪屈服面和孔隙比演化至关重要。这将土壤行为构建为一个非线性状态空间系统，其中隐藏的内部状态控制可观测的宏观应力-应变和体积响应。状态向量在空间上变化并随载荷增量动态演化。xi={e,ec,p′,pi′,q,σij,?ijp,??qp,??vp,ξ3,tc,ξ3,te,??3,??1,ψ,ψi,η,ηc,θ,M,Mtc,Mi,Mi,tc,Mi,te,H,N,χ,χtc,K,λ,Γ,G,ν,Dp,Dtcp,Dtep,a,b,c}

为了进一步实现高水平的宏观土壤建模保真度，NorSand依赖于以下关键本构机制：
1. **内部状态变量**：状态参数ψ=e?ec量化了相对于临界状态的土壤密度，ψi针对应力进行调整，Mi通过状态依赖性修正摩擦比以控制剪切强度。
2. **临界状态**：临界状态线（CSL），ec=Γ?λln(p′)，将孔隙比与平均有效应力联系起来，而临界状态摩擦比Mtc在三维应力状态下随Lode角θ变化。
3. **屈服面**：随动员摩擦比η=q/p′动态演化，该比值通过镜像条件Mi依赖于状态参数ψi。其最大范围受由χ和ψi定义的内部帽面限制，从而实现对剪切强度的状态依赖性控制。
4. **硬化规律**：模量H控制屈服面随塑性应变率演化的速率，实现应力状态依赖性的硬化/软化。
5. **应力-剪胀关系**：剪切应变率和体积应变率通过应力比和状态参数联系起来，其中N和χtc控制状态依赖性的剪胀和收缩。
6. **弹性**：剪切模量G和体积模量K通过泊松比ν相关联。

与上述描述的高保真度、状态依赖性NorSand公式相比，后者追踪38-45个空间变化的土壤状态量，而经验Bekker型模型（仍广泛用于板沉陷和原位剪切试验）依赖于少量随土壤应力、应变、孔隙比、加载历史、水-力耦合和试验几何形状变化的曲线拟合系数：Bekker-Wong参数：{kc,k?,n,c,?,K}。这些参数不是状态依赖性的，且与空间分辨的应力-应变场和边界条件断开，限制了其在一般加载条件下捕捉土壤变形或识别固有材料特性的能力。此外，常用的表观土壤强度参数如内聚力c和摩擦角?并非真正的材料常数，而是随密度、含水率、孔隙压力和加载历史变化的响应变量（Hvorslev，1960；Schofield和Wroth，1968）。存在更基本的底层土壤参数，如Hvorslev将内聚力分解为真实ct(e)、结构cs(e)和毛细cc分量，这些分量是破坏时孔隙比e的函数。类似地，摩擦角可以用临界状态摩擦比Mtc表示，从中可以推导出动员和峰值摩擦比ηi作为状态和孔隙比的函数。这种可变性和非不变性在Bevameter试验中由于非标准化设备、边界条件效应和试验装置敏感性而进一步放大（Kruger等人，2023）。

NorSand通过使用更基本的底层临界状态摩擦比Mtc解决了表观土壤参数（如总摩擦角?）的局限性，对于给定土壤，Mtc保持恒定，而动员峰值摩擦比ηmax随当前密度状态ψ演化。随着土壤压实，ηmax增加，使模型能够捕捉密度依赖性的峰值土壤强度和随后的峰值后向临界状态的软化。Mtc也不受水-力耦合或加载历史的影响，反映了真正的固有土壤特性，这是Bekker型公式在其当前形式下无法达到的物理保真度水平。

最后，为了解决任意边界值问题的NorSand本构关系，采用了非线性、大应变、准静态隐式（即误差控制）有限元公式。域在积分点处离散化并定义本构行为。牛顿法利用从NorSand变量的偏导数组装的一致切线矩阵来最小化内力和外力之间的残差。应力回映算法强制屈服面约束，而采用大应变理论处理几何非线性的增量载荷步进法。带有额外应变能残差的线搜索增强了收敛稳健性。NorSand的动态屈服面和软化特性使标准隐式求解器的收敛变得具有挑战性，需要先进的收敛控制策略。

**章节摘要**

**土壤特性**：本研究所用的土壤是一种天然细石英砂，称为Cullinan砂，平均粒径为149μm，颗粒呈棱角形。它来自南非比勒陀利亚以东50公里处的一个砂矿。Cullinan砂是比勒陀利亚大学用于岩土工程、土力学和离心机实验的标准化测试介质，配有室内和室外土壤箱。因此，它在以往的研究中被广泛使用（Archer和……）。

**NorSand参数估计结果**：表3总结了标定的NorSand参数，按本构类别（弹性、临界状态、塑性和状态变量）组织。表4列出了用于泛化表3中关键参数的相关元模型的功能形式和系数。

**实验室三轴试验预测（训练数据）**：CU三轴标定试验的有限元模拟结果如图9所示。该图比较了NorSand有限元预测与两种松散（左……）排水三轴试验的实验结果。

**结论与建议**：本研究建立了一个基于第一性原理的物理基础本构建模框架，用于利用NorSand临界状态土壤模型模拟原位板沉陷和剪切试验。作为经验Bekker-Wong方法的替代方案，本研究采用基于临界状态土壤力学（CSSM）的方法论，明确捕捉密度依赖性塑性、体积耦合和动态演化屈服面——这些现象在当前经验模型中并未明确表示。值得注意的是，该模型……

**CRediT作者贡献声明**：R. Kruger：撰写-原稿、可视化、验证、软件、方法论、调查、正式分析、数据管理、概念化。P.A. Crous：撰写-审阅与编辑、数据管理。Y. Trinidad：撰写-审阅与编辑、监督、项目管理。M.Z. Tekeste：撰写-审阅与编辑。

**利益冲突声明**：作者声明不存在任何可能影响本文所报告工作的已知竞争性经济利益或个人关系。

R. Kruger | P.A. Crous | Y. Trinidad | M.Z. Tekeste