动态步长迭代法在CSUH模型计算中的应用研究¹

摘 要：CSUH模型能够合理描述不同类型土的应力-应变特性，反演参数时需要大量迭代计算，耗时较长。为提高CSUH模型的计算效率，给出了弹塑性本构模型柔度矩阵和不同应力路径下应力应变增量关系式，避免了迭代过程中求解方程的步骤。还提出了一种分段调整迭代增量步长的动态迭代方法。该方法根据应力-应变曲线斜率的变化特征，在曲线斜率较大区域采用较小的迭代步长，在斜率较小区域采用较大的迭代步长，使总轴向应变的30%区间呈等差数列步长分布，70%区间呈均匀步长分布。通过与传统均匀步长迭代法对比，结果表面在保证计算精度的前提下动态步长迭代法可将迭代步数从3000步减少到50步，计算时间降低约47倍。最后，采用该方法对Fujinomori粘土、钙质砂和长河坝堆石料的CSUH模型参数进行了反演，预测结果表明CSUH模型能准确描述不同类型土体在不同应力路径下的力学行为。研究成果可为提高CSUH模型及其他弹塑性模型迭代效率提供借鉴。

关键词：CSUH模型；动态迭代步长；柔度矩阵；增量关系式；土体力学特性

1 引言

土的本构关系或称本构模型是以土的应力应变规律为研究对象，反映土的复杂应力应变规律的数学表达式，是岩土工程学科的重要理论基础，在大坝、铁路等工程中有重要应用^[1-4]。修正剑桥模型（Modified Cam-Clay, MCC）模型^{[5, 6]}作为经典的弹塑性本构模型，在描述正常固结粘土的力学特性方面取得了巨大成就。然而，MCC模型在描述超固结粘土时存在弹塑性转变处应力-应变曲线突变、预测峰值强度偏高等问题^[7]。为了更合理地描述超固结饱和粘土的应变软化和剪胀现象^[7]，Yao等^[8-10]在MCC模型基础上提出的统一硬化（Unified Hardenning model, UH）模型。UH模型仅修改了MCC模型的硬化参数，便可以统一地描述不同密实度粘土（正常固结和超固结粘土）的应变硬化和软化、剪缩和剪胀现象，具有广泛的适用性。工程中不仅有粘土，也有粒状土，比如砂土和堆石料、道砟等，但适用于粘土的UH模型并不能合理地描述粒状土的应力应变特性。因此，Yao等^[11]基于UH模型又提出了能同时描述粘土和粒状土的统一硬化（UH model for Clays and Sands, CSUH）模型，且CSUH模型可退化到UH和MCC模型，是MCC模型的理论提升和创新性发展^[12]。

将CSUH模型应用于真三轴试验预测和数值模拟前需要标定模型参数，反演方法标定参数时需要计算本构模型理论曲线值，这包含了大量数值计算。在弹塑性本构模型的数值积分方法研究中，Wu^[13]提出了一个统一的数值积分公式(UNIF)，该方法虽然包含了多种积分方案，但仅适用于完全塑性von Mises模型，难以推广到更复杂的本构模型。为了解决数值积分过程中应力点偏离屈服面的问题，Potts和Gens^[14]对5种屈服面漂移校正方法进行了评估，建议采用考虑弹性应变变化的校正方法。但该方法需要在每个积分子步后进行迭代校正，用于反演参数则计算效率较低。Sloan等^[15]提出了一种带有自动误差控制的显式积分算法，通过估计局部误差来自动划分应变增量子步，可用于复杂的弹塑性本构模型。然而，该方法在确定子步大小时仅考虑了局部误差，未充分考虑应力-应变曲线斜率变化对所需迭代步长的影响，导致某些区域积分精度偏低或计算效率不高。Yin等^[16]编写了包含多种本构模型的参数标定软件ErosOpt，标定参数时需调用类似ABAQUS的用户材料子程序（UMAT）计算本构模型理论曲线值。上述弹塑性本构模型的数值积分方法可用于高精度的UMAT子程序编写，但这种有限元方法需组装单元矩阵，每次迭代时都需要根据边界条件解方程组，浪费了大量算力。

在反演参数时，另一种计算本构模型理论曲线的方法是显式迭代^{[12, 17]}。目前常用的迭代方法包括前向欧拉法、四阶龙格库塔法等。其中，四阶龙格库塔法虽然具有较高的计算精度和稳定性，但计算量大；前向欧拉法实现简单，为保证计算精度需要采用较小的迭代步长。这些方法在实际应用中均存在计算效率不高的问题。

土的应力-应变曲线通常表现为非线性特征，其曲线斜率在不同阶段存在显著差异。现有求解模型理论曲线方法多为有限元方法或显式均匀步长迭代方法，未充分考虑应力-应变曲线斜率变化对所需迭代步长的影响。为提高CSUH模型的计算效率，本文提出了一种分段调整迭代增量步长的动态迭代方法。该方法根据应力-应变曲线斜率的变化特征，采用不同的迭代步长，在保证计算精度的同时显著减少迭代步数。

本文首先介绍了一般弹塑性本构模型的柔度矩阵和不同应力路径下应力应变增量迭代式，之后给出了CSUH模型的基本理论框架，然后详细阐述了所提出的动态迭代方法的实现过程。通过与传统均匀步长迭代方法的对比，验证了动态步长迭代法的有效性。最后，采用该方法对Fujinomori粘土、钙质砂和长河坝堆石料的CSUH模型参数进行了反演，预测结果表明CSUH模型能准确描述不同类型土体在不同应力路径下的力学行为，研究成果可为提高CSUH模型及其他弹塑性模型迭代效率提供借鉴。

2 弹塑性本构关系理论

2.1 一般弹塑性本构关系

土的弹塑性本构关系中弹性部分一般采用胡克定律，塑性部分由屈服条件、流动法则和硬化规律组成。屈服条件可以用屈服面方程f定量刻画，是确定开始产生塑性变形的应力条件，屈服面方程f是应力σ和硬化参量H的函数：

流动法则可以用塑性势面方程g定量刻画，是确定塑性应变增量方向的法则。塑性势面方程g一般可写为下面的形式：

式中，C为常数。

硬化规律实质是要确定硬化参量H，是确定塑性应变增量大小的规律，硬化参量H一般是塑性应变的函数：

式中，为塑性应变，不同的本构模型采用的塑性应变也不尽相同，比如修正剑桥模型采用塑性体应变$$作为硬化参数，此时。CSUH模型采用含应力比系数的塑性体应变作为硬化参数，此时，η表示应力比，$$和$$分别表示潜在破坏应力比和特征状态应力比，后面具体介绍。

2.2 弹塑性应力应变关系

经典的小变形弹塑性理论将应变ε分为弹性应变和塑性应变两部分，弹性应变按胡克定律计算，塑性应变按塑性理论计算。在采用增量法时有：

假设材料在弹性区间的本构关系为

式中，为弹性柔度矩阵，当弹性应变增量向量和应力增量向量分别为主应变增量和主应力增量时，为3阶柔度矩阵。当弹性应变增量向量和应力增量向量维度为6行1列时，为6阶柔度矩阵。

塑性应变增量不产生应力增量，由式可知总的应力增量和总的应变增量的关系为：

或者写为：

式中，为塑性柔度矩阵。采用增量塑性理论，对非关联流动法则可写为：

式中，是表示塑性应变增量大小的一个标量，称为塑性标量因子。

在塑性加载时必须保持增量一致性df=0，即

又因为硬化参数H往往视为塑性应变或累积塑性应变的函数，故有

将式代入：

将式代入式得到：

将式和式代入式：

式也可以写为弹塑性柔度矩阵形式的应力应变增量方程：

2.3 不同边界条件下的应力应变增量迭代式

使用主应力、主应变就可以模拟室内真三轴和常规三轴试验中常见变量。求解不同边界条件或称不同应力路径下本构模型应力应变增量迭代式的思路是求解由3阶柔度矩阵形式的应力、应变增量方程式、中主应力系数方程和应力路径方程组成的方程组，用最大主应变增量表示中主应变增量、最小主应变增量、最大主应力增量、中主应力增量和最小主应力增量。由式可得3阶柔度矩阵所列应力、应变方程：

当然，式也可以写为以下形式：

在试验中一般最大主应变增量

式中，b为常数，作为初始条件给出。那么在增量过程中同样保持此关系：

增量式与模型无关，不同的模型使用不同的柔度矩阵即可。加载时和卸载时分别取弹塑性和弹性柔度矩阵。下面给出一些常用应力路径的增量式。

（1）等向压缩（ISO）应力路径

在等向压缩应力路径中，通过施加静水压力使试样各向均匀压缩。三个主应力方向应力增量相等，应变增量也相等。在此使用三个应力增量相等的关系，可知中主应力系数b=0，，再结合式可得到含五个方程共六个变量的方程组：

式中，中主应力系数b是根据实验可知的已知量；柔度矩阵元素是本构模型给出的；由此可知

若剪切初始状态（初始应力p₀和初始孔隙比e₀）已知，MCC和CSUH模型根据式就可以进行迭代计算了。

（2）侧限压缩（K0）应力路径

侧限压缩试验就是将土样放在限制土样侧向变形的铁环内进行排水压缩的试验。由此可知其边界条件为，由此可推测侧向应力，即中主应力系数b=0。类似于等向压缩应力路径列方程组并解方程组可得：

（3）等排水剪切应力路径

等剪切应力路径就是在剪切过程中保持

常规三轴固结排水剪切（CD）应力路径的试验就是最常见的等

（4）围压不变的不排水剪切应力路径

不排水剪切应力路径比较特殊，是体积应变控制路径。其边界条件是体变增量，类似于等向压缩应力路径列方程组并解解方程组得：

式中A、B、和D分表示：

A、B、和D并没有物理意义，仅用作缩短式长度。式中b=0时退化为常规三轴固结不排水（CU）应力路径。

（5）等p排水剪切（P0）应力路径

在剪切过程中保持有效平均主应力p不变的剪切应力路径是等p排水剪切应力路径。在此应力路径下，类似于等向压缩应力路径列方程组并解方程组可得：

至此，常见的5种应力路径的应力应变增量迭代式已全部给出。

3. CSUH模型

3.1三维化模型框架

姚仰平等^[11]基于MCC模型提出了能同时描述超固结粘土和砂土的统一硬化模型（CSUH）。CSUH模型能很好地描述粘土、砂土、堆石料等多种类型土在不同初始条件下多种应力路径中的应力应变特性，是MCC模型的理论提升和创新性发展^[12]。

要将本构模型应用于真三轴室内试验预测或工程实践中，需要将本构模型三维化。变换应力（TS）方法是Yao等^[18]提出的一种本构模型三维化方法，其主应力形式如下：

式中，、p、q分别表示真实应力空间中主应力、平均主应力和广义剪应力；、分别表示变换应力空间中主应力和广义剪应力，根据强度准则不同取不同值；真实应力空间中p与变换应力空间中相等。

由于该方法可以基于不同的强度准则(如SMP准则、Mohr-Coulomb准则、Lade准则、广义非线性强度准则等)进行构建，因此能够适用于具有不同强度特性的材料。图1是基于SMP强度准则的变换应力方法，图1（a）为真实应力空间中SMP强度准则和广义Mises强度准则示意图。认为SMP强度准则在真实应力空间，变换后成为的广义Mises强度准则在变换应力空间，并认为本构模型是建立在变换应力空间中的，计算真实应力应变时只需要对屈服面求导到真实应力空间；图1（b）为π平面内上述转化示意图；图1（c）为CSUH模型屈服面在真实主应力空间和变换主应力空间中示意图；图1（d）为图1（c）的p-q截面，其中M_e为三轴伸长时的临界状态应力比。

图1 变换应力方法

利用TS方法将CSUH模型三维化时认为CSUH模型是建立在变换应力空间里的。CSUH模型在变换应力空间中塑性势面g：

式中，、分别为变换应力空间中有效平均主应力和剪应力，“~”符号表示变量在变换应力空间中，余同；p_y为塑性势面与p轴右侧交点横坐标；为特征状态应力比，几何意义为椭圆的上象限点处应力比。

CSUH模型在变换应力空间中的屈服面f：

式中：临界状态应力比M为p-q空间内常规三轴压缩临界状态线（CSL）斜率；临界状态参数χ可以改变屈服面形状，当χ=0时f退化为椭圆形屈服面，如图1（d）所示；p_s为压硬性参数；为塑性压缩系数；参数λ和κ分别为等向压缩试验中压缩指数和回弹指数；p_x0为初始条件下屈服面与p轴右侧交点横坐标；为硬化参量，

式中，为变换应力空间中当前应力比，为潜在破坏应力比，为塑性体积应变增量。剪切过程中，特征状态应力比和潜在破坏应力比都随着表示当前密实度大小的状态参量的改变而改变，如式、所示。由式、可见也是状态相关的。

式中，剪胀性参数m可以缩放的影响进而影响。状态参量如下：

式中，p_s是压硬性参数：

CSUH本构模型有M，ν，κ，λ，N，Z，χ，m共8个基本参数，不同路径和不同初始状态（初始围压p0和初始孔隙比e0）的同种土均使用同一套参数。

3.2 模型3阶柔度矩阵

根据式可以给出3阶柔度矩阵形式的增量关系式：

对于各向同性材料：

式中，E为弹性模量、ν为泊松比。

式中3阶塑性柔度矩阵：

式中，模型的屈服面对真实应力和硬化参量求偏导，塑性势面对变换应力求偏导，硬化参量对真实应变求偏导数，以上求偏导过程要用到链式求导法则。屈服面、塑性势面一般写为和的函数，而和又是应力的函数，根据基于SMP强度准则的TS方法：

式中，I₁、I₂和I₃分别是第一、第二和第三主应力不变量；其中q_c为基于SMP强度准则的变换广义剪应力。文献[11]给出了CSUH模型的式求导结果，在此不在赘述。式也可以写为：

硬化参量对塑性应变的偏导也可以写为对塑性体应变的偏导，因此对于CSUH模型，式也可以写为：

加卸载准则如下：

如果式中通过弹性柔度矩阵试算出的说明加载，需要使用式的弹性柔度矩阵加塑性柔度矩阵即弹塑性柔度矩阵重新计算式中的应力应变增量，如果通过弹性柔度矩阵试算出的则说明当前过程为卸载，试算正确，即仅使用式中弹性柔度矩阵计算正确，可以进行下一步的弹性试算。

4 动态迭代方法

4.1 常见迭代方法

利用优化算法反演本构模型参数需要大量计算，其中类似式的应力应变增量迭代式最为耗时。优化反演方法确定本构模型参数的过程就是使用优化算法向本构模型提供不同的本构模型参数，本构模型根据参数和试验提供的初始条件计算相应试验的理论应力应变曲线，然后比较计算的理论应力应变曲线与试验应力应变曲线间的误差并反馈给优化算法，优化算法根据反馈的误差重新调整并尝试向本构模型提供更好的模型参数。因此，优化算法反演本构模型参数需要大量使用本构模型应力应变增量式计算应力应变曲线。

前面介绍了不同本构模型可以使用类似式统一的应力应变增量式迭代计算应力应变曲线，迭代时可使用前向欧拉法、四阶龙格库塔法、隐式迭代等多种方法。四阶龙格库塔法虽然具有较高的计算精度和较好的稳定性，但其复杂的实现过程和较大的计算量使其在实际应用中面临挑战。隐式迭代方法（如后向欧拉法、牛顿-拉普森法等）虽然允许采用较大的增量步，但需要迭代求解非线性方程，且可能存在收敛性问题。半隐式返回映射法在处理弹塑性问题时表现出色，但其程序实现复杂，计算效率相对较低。

前向欧拉法因其简单直观、易于实现的特点而被广泛应用。该方法基于一阶泰勒展开，通过当前步的切线模量直接计算下一步的应力状态，其核心思想是将连续的应力应变关系离散化为增量形式。按前向欧拉法原理可将式及之后各应力路径下应力应变增量简化为得到：

式中，表示第n+1步的应力和应变，表示第n步的应力和应变；C_n表示根据第n步应力状态计算的柔度矩阵；表示最大主应变增量。

四阶龙格库达法（RK4）积分精度较高，前向欧拉法积分精度较低，但发现将总的轴向应变均匀分为3000份时，两种迭代方法在3000步均匀步长所得的积分曲线近似重合，如图2所示，由此可知均匀步长迭代3000步时能得到较为精确的应力应变曲线。由于RK4需要求解4次不同增量子步长下的柔度矩阵，而前向欧拉法仅需要计算1次，前向欧拉法计算量更小。因此，当迭代步数达到3000步时，前向欧拉法即能节省计算量，又能获得令人满意的结果，因此选用前向欧拉法计算本构模型应力应变曲线。相比于优化算法自身的计算，本构模型应力应变增量式的大量迭代计算才是整个反演过程中最为耗时的部分，前向欧拉法迭代3000步仍存在耗时太长的问题。

图2 不同积分方法以均匀步长迭代3000步得到的CU应力路径应力应变曲线

4.2 提出的动态迭代方法

土的应力应变曲线的斜率并不是均匀的，可以通过分段调整模型迭代增量步长减少总的迭代步，从而加快反演参数的速度。

曲线斜率较大时较小的轴向应变增量Δε₁就能引起较大的应力增量，而曲线斜率较小时较大轴向应变增量Δε₁只能引起较小的应力增量。如果不同曲线斜率时都用相同的迭代增量步长Δε₁，为了保证斜率较大时迭代的精确性，迭代增量步长Δε₁必须设置的足够小，所需迭代步数就很多，计算速度就很慢。为了提升本构模型数值计算速度，提出了分段调整模型迭代增量步长的动态迭代方法以减少总的迭代步数。在应力应变曲线斜率较大时采用较小的迭代步长，在应力应变曲线斜率较小时采用较大的迭代步长，这种做法可以在保证精度的同时减少迭代步数，加快了计算速度。呈等差数列分布的轴向应变一般占总轴向应变的30%，呈均匀分布的轴向应变一般占总轴向应变的70%，分段位置随总的轴向应变的大小变化而动态变化。

图3 分段调整本构模型迭代增量步长Δε₁

如图3所示，将总的轴向应变ε_1max分为两部分，第一部分为等差数列划分部分，此部分轴向应变S1_n与总的轴向应变ε_1max的占比为p_t，这部分轴向应变增量Δε₁呈现等差数列分布；第二部分为均匀划分部分，此部分轴向应变S2_n与总的轴向应变ε_1max的占比为(1-p_t)，这部分轴向应变增量Δε₁呈现均匀分布。如果已知总的迭代步数为N_s，均匀划分部分迭代步数N_s2=p_tN_s，则均匀划分部分轴向应变增量Δε₁：

设等差数列划分部分的初始迭代步长a₁为均匀划分部分迭代步长的1/10：

由总迭代步数N_s和均匀划分部分迭代步数N_s2可知等差数列划分部分迭代步数N_s1=(1-p_t)N_s，则等差数列划分部分的等差数列差值d和迭代增量步长a_n：

取合适的比例p_t和合适的总迭代步数N_s即可用比较少的迭代步数计算出精度较高的应力应变曲线，如图4所示，图中取p_t=0.3，N_s=50。

图4 均匀迭代增量步长和等差数列分布的迭代增量步长计算结果对比

表1 CSUH模型参数敏感性分析中基准参数^[11]

图4为CSUH模型采用表1中丰浦砂参数计算围压为500kPa和孔隙比为0.84时CU试验的应力应变曲线。根据式可知，理论上数值计算迭代增量步长dε₁应该无限小，但在实际计算中发现一般3000步已足够精确，如图2所示，此时，由于此时步长不是无限小，故将dε₁记作Δε₁，其余主应变增量和主应力增量也是如此。可以看出，如果采用50步均匀步长迭代计算，则与3000步均匀步长迭代计算的应力应变曲线相差甚远，偏差较大。如果采用分段调整模型迭代增量步长的动态迭代方法，同样是迭代50步，但其计算的应力应变曲线与均匀迭代3000步计算的应力应变曲线很接近。

可见，通过分段调整模型迭代增量步长的动态迭代方法可在基本保证精度的情况下将迭代步数降为原来的1/60。以反演长河坝堆石料的CSUH模型参数为例，反演过程中使用了1个等向压缩和4个CD试验，数据点共414个，当迭代步数为50步时用时7分钟，当迭代步数为3000步时用时332分钟，用时降为原先的1/47，如图5所示。

图5 不同迭代步数反演计算时长

可见，分段调整模型迭代增量步长的动态迭代方法在保证预测精度的前提下极大地减少了模型迭代次数，大幅提升了反演参数的速度。

5 CSUH模型迭代反演参数

5.1 Fujinomori黏土

为了验证CSUH模型描述不同类土应力应变特性的能力，采用动态迭代方法预测了Fujinomori黏土的CSUH模型参数，其初始条件表2所示。

表2 Fujinomori 黏土剪切初始状态

对于OCR=8.0的试样，p₀=98 kPa，过程为先固结到前期固结压力σ_c=784 kPa，然后卸载到p₀=98 kPa，最后在p₀=98 kPa下进行等p剪切试验，其余超固结度值的粘土经历的加卸载过程类似。CSUH模型参数如表3所示。

表3 Fujinomori 黏土的 CSUH 模型参数

（1）等p压缩试验

如图6所示，图中点代表试验结果，线代表CSUH模型的预测结果。其中，超固结度（OCR）最小为1.0，此时对应为正常固结土，密实度最小；超固结度（OCR）最大为8，此时对应为重度超固结土，密实度最大。由图6可以看出，正常固结粘土表现为应变硬化和剪缩。随着超固结度的增大，峰值应力比逐渐增大，体积变形开始由剪缩转变为先剪缩后剪胀。当OCR=8时，峰值应力比最大，体积应变的剪胀量最大。同时从图中可以看出CSUH模型计算结果与不同超固结度粘土的试验结果比较吻合，这说明CSUH模型能够合理地定量描述不同密实度土在剪胀性方面的区别。

图6 Fujinormori粘土等p压缩试验和CSUH模型预测结果^[19]

Fig. 6 Comparisons between the constant p compression test results of the Fujinormori clay and the CSUH predictions

（2）真三轴试验

Fujinomori粘土真三轴试验基本过程是：先等向固结到p=196 kPa，然后以1%/天的轴应变率进行等p（或称等）剪切试验。从图7中可以看出CSUH模型的预测结果与试验结果符合较好，也说明采用变换应力方法后CSUH模型能够定量描述不同应力洛德角θ时土在摩擦性方面的区别。

(a) 中主应力比b=0，等同于应力洛德角θ=0°

(b) 中主应力比b=0.268

(c) 中主应力比b=0.5

(d) 中主应力比b=0.732

(e) 中主应力比b=1

图7 Fujinomori粘土不同中主应力系数b对应的等p路径真三轴试验测量结果和CSUH模型预测结果^[20]（点：测量；线：预测）

Fig. 7 Comparisons between the true triaxial test measured results under constant p stress path with different medium principal stress coefficient b of the Fujinormori clay and the CSUH predictions

总之，对饱和Fujinomori粘土的常规三轴和真三轴试验进行预测结果表明了CSUH模型不仅可以定量区分不同密实度粘土剪胀性方面的差异，还能定量区分不同应力洛德角时粘土在摩擦性方面的差异。

5.2 钙质砂

翁贻令^[21]对南沙群岛吹填钙质砂（或称珊瑚礁砂）做了CD和CU试验。钙质砂地基以含砾砂土或碎石土为主，为未胶结的松散体，在此用于验证模型的为粒径小于0.075mm的钙质砂试验数据。CSUH模型参数如表4所示，对比效果见图8和图9。

表4 钙质砂的CSUH模型参数

Table CSUH model parameters for calcareous silt

(a) 轴向应变与偏应力

(b) 轴向应变与体积应变

图8 钙质砂常规三轴固结排水（CD）试验测量结果与CSUH模型预测结果^[21]

Fig. 8 Comparisons between the CD test measured results of the calcareous silt and the CSUH predictions

由图8可知，不同围压下钙质砂在CD试验中均发生了应变硬化和剪缩。

(b) 轴向应变和超静孔压

图9 钙质砂常规三轴固结不排水（CU）试验测量结果与CSUH模型预测结果^[21]

Fig. 9 Comparisons between the CU test measured results of the calcareous silt and the CSUH predictions

由图9可知，钙质砂在较小的应变时就达到了峰值强度，随后应变软化到强度几乎为零，孔压几乎接近轴向应力，土样接近液化状态。因此，若整理有效应力比可以发现测量的有效应力比出现了异常，这点文献[21]也有提到。

4.3 预测长河坝堆石料

刘恩龙等^[22]对长河坝的主堆石料（硬质闪长石构成）进行了等向压缩、CD和CU试验。长河坝堆石料的CSUH模型参数如表5所示。

表5 长河坝堆石料的CSUH模型参数

Table 8 CSUH model parameters for the rockfill material from the Changhe dam

图10 长河坝堆石料等向压缩试验测量结果与CSUH模型预测结果^[22]

Fig. 10 Comparisons between the isotropic compression test measured results of the rockfill from the Changhe dam and the CSUH predictions

图10是长河坝堆石料的等向压缩试验结果和CSUH模型的预测结果。从预测结果与试验结果的比较可以看出CSUH模型能够描述堆石料这种大颗粒粒状土的压硬性，压硬性即土在压缩过程中所表现出的模量随密度增加而增大的特性，即“越压越硬”^[23]。

(a) 轴向应变与偏应力

(b) 轴向应变与体积应变

图11 长河坝堆石料常规三轴固结排水（CD）试验测量结果与CSUH模型预测结果^[22]

Fig. 11 Comparisons between the consolidated drained triaxial compression (CD) test measured results of the rockfill from the Changhe dam and the CSUH predictions

图11对比的是长河坝堆石料在不同围压下的CD试验结果和CSUH模型预测结果。从试验结果可以看出在围压较低时土体会有轻微的应变软化，体积会发生剪胀，在围压较高时土体表现为应变硬化，体积始终保持为剪缩。这种现象主要是由堆石料的颗粒破碎造成的。由于堆石料与砂土都具有破碎特性，因此堆石料与砂土的试验结果定性上基本一致，但是堆石料比砂土更容易破碎，因此会造成二者的试验结果在定量上有所差别。例如，图8的钙质砂在初始状态为（400 kPa，0.982）和图11的堆石料在初始状态为（3500 kPa，0.367）在剪切后得到的最终体变都在6%左右。可以推论，如果初始状态相同，以上钙质砂和堆石料的最终体变量必不相同。

(a) 轴向应变与偏应力

(b) 轴向应变与超静孔压

图12 长河坝堆石料常规三轴固结不排水（CU）试验测量结果与CSUH模型预测结果^[22]

Fig. 12 Comparisons between the consolidated undrained triaxial compression (CU) test measured results of the rockfill from the Changhe dam and the CSUH predictions

图12对比的是长河坝堆石料在不同围压下的CU试验结果和CSUH模型预测结果。从图12（a）中可以看出不同围压下土体均表现为应变硬化。从图12（b）中可以看出随着轴应变的增加孔压快速增大后逐渐减小，即堆石料在不排水过程中有先剪缩后剪胀的趋势。另外，围压越大产生的孔压越大，孔压抵消一部分围压，所以相同围压p₀下不排水强度q_f小于排水强度。对比图11（b）和图12（b）可知，初始围压为3500 kPa的堆石料在排水试验中一直剪缩，而在不排水试验中还出现了剪胀。较好的预测结果表明CSUH模型也可以定量描述不同排水条件下同种材料（堆石料）剪胀性的区别。

6 结论

本文提出了一种分段调整迭代增量步长的动态迭代方法，用于提高CSUH模型的参数反演效率。通过理论分析和试验验证，得到以下主要结论：

(1)给出了弹塑性本构模型的弹塑性柔度矩阵和不同边界条件下应力应变增量关系式，在计算模型理论曲线值时不必再求解方程组，加快了计算速度。

(2)基于应力-应变曲线斜率变化特征，将总轴向应变的30%区间采用等差数列分布的迭代步长，70%区间采用均匀分布的迭代步长，可以有效平衡计算精度和效率。研究表明，该方法可将迭代步数从传统的3000步减少到50步，计算时间降低约47倍，同时保持了与四阶龙格库塔法相当的计算精度。

(3) CSUH模型能准确描述不同类型土体（如黏土、砂土及堆石料）在不同应力路径（等向压缩、不同中主应力比的等p排水、CD和CU）下的应力-应变特性。

本研究提出的动态迭代方法显著提高了反演弹塑性模型参数时迭代效率，为该模型在工程实践中的应用创造了有利条件。后续研究建议进一步探讨该方法在其他本构模型中的应用。

参考文献(References):

[1] 姚仰平, 罗汀, 侯伟. 土的本构关系[M]. 北京: 人民交通出版社股份有限公司, 2018.

[2] Zhang, Y., Chen, Y. A constitutive relationship for gravelly soil considering fine particle suffusion[J]. Materials, 2017, 10(10): 1217.

[3] Cong, S., Ling, X., Li, X., et al. Elastoplastic Model Framework for Saturated Soils Subjected to a Freeze–Thaw Cycle Based on Generalized Plasticity Theory[J]. Materials, 2021, 14(21): 6485.

[4] Dong, L., Tian, S., Yao, C., et al. A Nonlinear Constitutive Model for Remoulded Fine-Grained Materials Used under the Qinghai–Tibet Railway Line[J]. Materials, 2022, 15(15): 5119.

[5] Roscoe, K. H., Schofield, A. N., Thurairajah, A. Yielding of clays in states wetter than critical[J]. Géotechnique, 1963, 13(3): 211-240.

[6] Roscoe, K. H., Burland, J. B. On the generalized stress–strain behaviour of ‘wet clay’. In: Engineering plasticity[M]. Cambridge, 1968.

[7] 侯伟, 姚仰平. 剑桥模型与统一硬化模型对超固结土特性描述的对比分析[J]. 工业建筑, 2011, 41(09): 18-23+140.

[8] Yao, Y. P., Gao, Z. W., Zhao, J. D., et al. Modified UH model: constitutive modeling of overconsolidated clays based on a parabolic Hvorslev envelope[J]. Journal of Geotechnical Geoenvironmental Engineering, 2012, 138(7): 860-868.

[9] Yao, Y. P., Hou, W., Zhou, A. N. Constitutive model for overconsolidated clays[J]. Science in China Series E: Technological Sciences, 2008, 51(2): 179-191.

[10] Yao, Y. P., Hou, W., Zhou, A. N. UH model: three-dimensional unified hardening model for overconsolidated clays[J]. Géotechnique, 2009, 59(5): 451-469.

[11] Yao, Y. P., Liu, L., Luo, T., et al. Unified hardening (UH) model for clays and sands[J]. Computers and Geotechnics, 2019, 110: 326-343.

[12] Zhu, B. L., Chen, Z. Y. Calibrating and validating a soil constitutive model through conventional triaxial tests: an in-depth study on CSUH model[J]. Acta Geotechnica, 2022, 17(8): 3407-3420.

[13] Wu, W. A unified numerical integration formula for the perfectly plastic von Mises model[J]. International Journal for Numerical Methods in Engineering, 1990, 30(3): 491-504.

[14] Potts, D., Gens, A. A critical assessment of methods of correcting for drift from the yield surface in elasto‐plastic finite element analysis[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1985, 9(2): 149-159.

[15] Sloan, S. W., Abbo, A. J., Sheng, D. Refined explicit integration of elastoplastic models with automatic error control[J]. Engineering Computations, 2001, 18(1/2): 121-194.

[16] Yin, Z. Y., Jin, Y. F. Development of Geotechnical Optimization Platform EROSOPT[M]. Practice of Optimisation Theory in Geotechnical Engineering. Springer. 2019: 243-292.

[17] Zhu, B. L., Su, X. X., Cao, Y. Q., et al. Determining parameters of the CSUH constitutive model by genetic algorithm[J]. Japanese Geotechnical Society Special Publication, 2020, 8(6): 188-193.

[18] Yao, Y. P., Wang, N. D. Transformed stress method for generalizing soil constitutive models[J]. Journal of Engineering Mechanics, 2014, 140(3): 614-629.

[19] Nakai, T., Hinokio, M. A simple elastoplastic model for normally and over consolidated soils with unified material parameters[J]. Soils and Foundations, 2004, 44(2): 53-70.

[20] Nakai, T., Matsuoka, H., Okuno, N., et al. True triaxial tests on normally consolidated clay and analysis of the observed shear behavior using elastoplastic constitutive models[J]. Soils and Foundations, 1986, 26(4): 67-78.

[21] 翁贻令. 钙质土的抗剪强度及其影响机制研究[D]. 广西大学, 2017.

[22] Liu, E. L., Chen, S. S., Li, G. Y., et al. Critical state of rockfill materials and a constitutive model considering grain crushing[J]. Rock Soil Mechanics, 2011, 32(S2): 148-154.

[23] 姚仰平, 张丙印, 朱俊高. 土的基本特性、本构关系及数值模拟研究综述[J]. 土木工程学报, 2012, 45(03): 127-150.