不同应力路径下旋转硬化方程的统一表示及验证

张浩，吴志鸿，张峰，戴自航

(1.中建海峡建设发展有限公司，福州 350003；2. 福州大学 土木工程学院, 福州 350116 )

原状土体是在一定应力条件下形成的，土体颗粒结构产生定向排列，在宏观结构上形成一定的各向异性。同时，土体的变形与加载应力路径密切相关，采用不同的应力加载路径会得到不同的应变，引起初始各向异性的变化，形成土体颗粒之间的新的排列，这种新的各向异性称为诱发性各向异性。在屈服面模型中，这种各向异性通过屈服面的旋转来表示，这就是模型中采用的旋转硬化表达式。

研究人员对旋转硬化规律采用了不同的表达形式。Wheeler等[1]提出的旋转硬化表达式同时考虑了塑性体积应变和塑性剪切应变的影响，但旋转硬化参数不能直接确定，需要通过反分析，并且没有考虑洛德角的影响。Whittle等[2]提出的旋转硬化表达式采用临界状态线作为旋转的边界线，只考虑了塑性剪切应变对旋转硬化的影响，这与试验现象不吻合。Taiebat等[3]没有考虑塑性剪切应变的影响，同时，参数值不能直接确定，需要通过反分析,屈服面的旋转存在边界线。Zhang等[4]依据交变加载试验结果提出了一个旋转硬化表达式，该表达式考虑了塑性剪切应变的作用，但是，参数需要通过反分析才能确定，给计算带来不便。Hashiguchi[5]认为屈服面旋转过程中塑性体积应变没有影响，采用K0线作为屈服面的倾斜线，这与已有的试验结果不符。Hueckel等[6]同时考虑了塑性体积应变和塑性剪切应变对旋转速度的影响，但是，旋转硬化参数的确定缺乏依据。Newson等[7]采用应力表示屈服面的旋转速度，屈服面的旋转只考虑了每一步的应力增量，没有考虑剪切应力的作用，同时，由于该本构模型采用非相关联流动法则，使得该旋转硬化表达式非常复杂。王立忠等[8]提出的表达式同时考虑了塑性体积应变和塑性剪切应变的作用，认为屈服面的旋转存在一定的范围，超过该范围时，对应屈服面的倾斜角度平衡值为0。对土体旋转硬化规律表达形式方面的研究成果还有很多[9-14]，目前提出的旋转硬化表达式均是针对某一种具体的应力路径所提出的，当应力路径发生改变时，旋转硬化表达式的适用性有待验证。

1 旋转硬化方程的统一表示

1.1 旋转硬化方程的分类

研究人员针对旋转硬化规律采用了不同的表达形式，不同表达形式之间，塑性剪切应变和塑性体积应变对应的平衡值也不相同。依据加载应力比η和相对应屈服面的倾斜角度平衡值α之间的关系，这些旋转硬化表达式可划分为3种类别。

第1种：η与α之间为线性关系，不存在旋转极限线，方程式为

α=xη

(1)

式中：x为η与α之间的线性比例关系参数，其取值一般为1和0，这种关系式比较典型的旋转硬化规律包括Whittle等[2]、Taiebat等[3]及Zhang等[4]提出的旋转硬化表达式。

第2种：η与α之间为线性关系，但是存在旋转极限线，超出旋转范围极限线之后，屈服面对应的倾斜角度平衡值α为0，表达式为

(2)

这种关系式比较典型的旋转硬化规律包括王立忠等[8]提出的旋转硬化表达式。

第3种：η与α之间为曲线关系，比较典型的包括Wheeler等[1]提出的旋转硬化规律表达式。

3(3η-4α)(M2-η2)=±8β(3α-η)(η-α)

(3)

为了分析这3种表达方式的合理性，收集了14种数据的应力比η/M和与应力比相对应的平衡值α/M之间的试验结果([1], [15-17])，M为考虑了洛德角影响的临界状态线的斜率，试验结果与采用这3种表达式进行模拟的结果对比如图1所示。

图1 归一化的应力比η/M与平衡值α/MFig.1 Normalization of the stress ratio

从图1可以看出，在归一化应力比η/M较低的情况下，除了第1种表达式中x取0的情况外，依据这3种表达式得到的模拟结果与试验结果之间都吻合得比较好，当η/M的值达到一定值之后，塑性剪切应变的影响程度开始增加，随着归一化应力比η/M的增加，α/M反而开始下降。采用第2种关系式得到的曲线虽然也能反映出这种趋势，但是模拟曲线与试验结果之间存在明显差距。

1.2 旋转硬化方程的统一表示及参数确定

1.2.1 旋转硬化表达式 从图1可以看出，第3种表达式能较好地反映α/M随η/M的增加而下降的趋势，但试验结果与模拟结果之间存在一定误差，尤其在拉伸状态下，这种误差比较明显。同时，通过试验数据也可以看出，在应力比η较低的情况下，塑性体积应变占主导位置，η与α之间相等，当应力比η趋向于M时，α接近于0。基于此，提出旋转硬化规律的统一表达式为

(4)

式中总共包含4个参数，分别是μ、β、a和b，当这4个参数取不同的值时，就可以得到不同的旋转硬化规律表达式。同时，从图1可以看出，在应力比η/M较低时，η/M近似等于α/M，当达到一定程度之后，α/M随着η/M的增大逐渐减小，最后趋向于0。为了反映这种趋势，本文提出参数a和b分别取1和0，得到的旋转硬化规律表达式为

(5)

为验证该旋转硬化方程的合理性，将该表达式嵌入修正剑桥模型，有关修正剑桥模型的具体介绍可参考文献[18]。建立了一个能描述各向异性的弹塑性本构模型，模型的屈服面方程表示为

f=(q-αp)-(M2-α2)(pc-p)p=0

(6)

1.2.2 旋转硬化表达式参数的确定 旋转硬化表达式包含两个参数，分别是μ和β，这里介绍对这两个参数的推导过程。若继续一维压缩过程，屈服面不发生旋转，即

(7)

一维压缩过程中

(8)

式中：dε1和dε3分别表示压缩过程中的单元体竖向应变和水平应变，联合式(7)和式(8)，得到

(9)

结合方程式(6)，一维压缩中塑性体积应变和塑性剪切应变之间的关系表示为

(10)

得到

(11)

将式(11)代入式(9)，可得到

(12)

如果将土体式样各向同性固结到先期压力的2～3倍左右，则土体的各向异性将基本消失[19]，此时,发生的塑性体积应变可表示为

(13)

各向同性固结时，应力比η=0, 塑性剪切应变和塑性体积应变之间的比例关系可表示为

(14)

将式(14)代入式(5),得

(15)

对式(15)左边积分，得

(16)

假定屈服面累计旋转到初始值α0的0.1倍时，土体显示各向同性性质，这样联合式(13)、式(15)和式(16)就可求出参数μ的表达式。

(17)

联合式(5)、式(12)和式(17)，可建立α/M和η/M之间的关系曲线，如图2所示。

图2 应力比η/M与平衡值α/MFig.2 Stress ratio η/M and equilibrium value

从图2可以看出，采用本文提出的旋转硬化表达式得到的曲线与试验结果之间吻合得较好，能很好地反映屈服面的倾斜角度先增后减，最后趋向于0的过程。将提出的旋转硬化表达式导入修正剑桥模型，模拟屈服面的旋转以及不同应力路径下的土体试样的应力-应变行为，模型需要输入的模型与修正剑桥模型需要输入的参数相同，总共包含6个，分别为e0、λ、κ、M、υ和pm0, 参数的物理意义可参考文献[18]。

2 不同应力路径试验的模拟

为了分析旋转硬化对不同加载路径下土体应力-应变关系的影响，选择3种不同的应力路径进行分析，分别是温州软土三轴排水试验[20]、BBC(Boston Blue Clay)土的三轴不排水试验[21]、Otaniemi土的循环加载试验[22]，这3种土需要输入模型的参数见表1。

表1 输入模型的参数Table 1 Input parameter of the model

2.1 三轴排水试验

沈恺伦[20]对温州软土进行了三轴排水试验，试验过程先依照原状土的应力状态进行固结，然后再分别进行排水状态下的三轴压缩试验。图3和图4表示的是应力路径为BCD50和BED33时，三轴压缩试验结果和采用不同旋转硬化表达式得到的模拟结果之间的对比。

图3 应力路径BCD50的试验结果和模拟结果Fig.3 Simulated and experiment results of

图4 应力路径BED33的试验结果和模拟结果Fig.4 Simulated and experiment results of

图3是应力路径为BCD50时试验结果和采用不同旋转硬化表达式时得到的模拟结果之间的对比，从图中可以看出，总体上模拟结果与试验结果之间吻合得比较好，采用不同旋转硬化表达式得到的模拟曲线之间非常接近，在模拟p-εv时，随着体积应变的发展，模拟曲线之间开始出现差距。当体积应变达到3%时，模拟曲线之间开始出现一定的差距，与试验结果之间的差距也慢慢增大，在模拟q-εs时，采用不同旋转硬化表达式得到的模拟结果与试验结果之间非常接近，都能较准确地模拟q-εs的发展趋势。

图4是三轴拉伸试验BED-33时试验结果与模拟结果之间的对比，从图中可以看出，在模拟p-εv曲线时，采用不同旋转硬化表达式对模拟结果的影响不是很明显，采用4种旋转硬化表达式得到的模拟曲线之间非常接近，都能较准确地模拟p-εv之间的发展趋势，但在模拟q-εs曲线时，模拟结果与试验结果存在一定的误差，当剪切应变较小时，模拟曲线之间比较接近，随着剪切应变的发展，模拟曲线之间开始出现一定程度的误差，采用本文提出的旋转硬化表达式得到的模拟结果与试验结果之间吻合较好。

2.2 三轴不排水试验

Ladd等[21]对BBC(Bonston Blue Clay)土进行了三轴不排水压缩试验，所取土试样首先各向同性固结到150 kPa, 然后根据不同的超固结比值Roc，分别卸载到相应的应力值，这里选取的Roc分别为1和4，土样在整个加载过程中保持不排水条件，需要输入模型的土体参数见表1。Yu等[23]也对该试验进行了数值研究。

图5 不同超固结比的试验结果与模拟结果Fig.5 Simulated and experiment results of samples with different overconsolidation

图5为模拟结果与试验结果之间的对比，从图中可以看出，在对应力路径p-q进行模拟时，采用不同旋转硬化表达式得到的模拟曲线与试验曲线之间比较接近，采用这4种旋转硬化表达式都能较准确地反映出p-q之间的变化趋势。在模拟q-εs时，当OCR为1，即当试样为正常固结时，模拟曲线位于试验曲线的上部，在出现相同的剪切应变εs时，模拟得到的剪切应力q要大于试验过程中的剪切应力，同时，采用不同旋转硬化表达式得到的模拟曲线之间也不是非常接近，相比较而言，采用本文提出的旋转硬化表达式得到的模拟曲线与试验曲线之间最为接近。当OCR为4时，在模拟应力路径p-q时，不同旋转硬化表达式得到的模拟结果与试验结果之间非常接近，都能准确地反映出应力路径p-q的变化趋势，在模拟q-εs时，模拟结果与试验结果之间存在一定的误差，但这种误差要小于OCR为1时试验结果与模拟结果之间的误差，总体而言，这4个旋转硬化表达式都能较准确地反映出变形趋势，采用本文提出的旋转硬化表达式得到的模拟结果与试验结果之间最为接近。

2.3 循环加载试验

循环加载试验选取Otaniemi重塑土进行了研究，这种土的相关参数及试验过程见Karstunen等[22]的研究成果，在三轴压缩之前，先将试样按照一定的应力路径η0进行固结，随后将试样进行卸载到一定值, 然后对试样进行加载卸载循环试验。本文选取了两组试样，分别是试样CAE3516R和试样CAE3519R，试样CAE3516R模拟结果如图6所示。

图6 CAE3516R试验结果与分析结果对比Fig.6 Simulated and experiment results of

从图6可以看出，对于CAE3516R试样，在模拟应力-应变曲线p-εv时，采用不同旋转硬化表达式对模拟结果有较大影响。初始加载时，采用Whittle等和王立忠等提出的旋转硬化表达式得到的模拟曲线位于试验曲线的上部，采用本文提出的旋转硬化表达式与Hueckel等提出的旋转硬化表达式得到的模拟曲线之间比较接近。卸载阶段，不同旋转硬化表达式得到的模拟结果之间的差距比加载阶段不同旋转硬化表达式得到的模拟结果之间的差距要明显一些，相比较而言，采用本文提出的旋转硬化表达式得到的模拟结果与试验结果之间更为接近。在模拟应力-应变关系q-εs时, 采用4种旋转硬化表达式都能反映出应力-应变关系的变化趋势，但是，相比较应力-应变关系p-εv而言，不同旋转硬化表达式得到的模拟曲线之间的差距更明显。

图7为CAE3519R的试验结果与模拟结果之间的对比，从图中可以看出，在模拟应力-应变曲线p-εv时，加载阶段，不同旋转硬化表达式得到的模拟结果之间存在一定差距，卸载阶段，不同模拟结果之间则比较吻合。在模拟应力-应变关系q-εs时，不同旋转硬化表达式得到的模拟曲线之间差距比较明显，采用本文提出的旋转硬化表达式得到的模拟曲线与实测数据之间最为接近。

图7 CAE3519R试验结果与分析结果对比Fig.7 Simulated and experiment results of

3 结论

对提出的旋转硬化规律表达式进行总结，依据加载应力比η和相对应屈服面的倾斜角度平衡值α之间的关系，将提出的旋转硬化表达式划分为3种类型，分别是线性关系同时不存在旋转极限线、线性关系同时存在旋转极限线以及曲线关系。建立了这3种关系式的函数曲线，函数曲线与试验结果之间的对比结果显示，采用曲线表达形式与试验结果之间最为吻合。

在采用曲线型表达方程的基础上，提出了一个适用于不同应力路径下的旋转硬化表达式，同时考虑了塑性体积应变和塑性剪切应变对旋转速度的影响。表达式共包含两个参数，通过理论分析，得到这两个参数的确定方法。将所提出的旋转硬化表达式导入修正剑桥模型，建立一个能描述各向异性的本构模型，采用此模型模拟了不同应力路径下的应力-应变关系，模拟过程中，同时采用了其他3种旋转硬化表达式。将试验结果与模拟结果之间以及采用不同旋转硬化表达式得到的模拟结果之间进行对比分析，分析结果验证了提出的旋转硬化表达式的合理性。