考虑颗粒破碎的钙质砂修正邓肯-张E-B模型*

曾凯锋 刘华北

(华中科技大学土木工程与力学学院， 武汉 430074， 中国)

0 引 言

钙质砂通常是指富含碳酸钙或其他难溶碳酸盐类物质的特殊介质，是海洋沉积物的一种，广泛分布于我国南海诸岛附近，其特殊的形成环境导致其颗粒疏松多孔、形状不规则、易发生破碎(刘崇权等， 1995； 汪稔等， 2019)。而由于颗粒破碎的存在，钙质砂的力学性质与普通陆源砂有很大差异。

现今，在颗粒破碎对岩土材料力学性质影响方面，国内外学者进行了大量的试验研究。吴京平等(1997)对钙质砂进行了三轴剪切试验，发现颗粒破碎将使得钙质砂剪胀性减小，体积收缩应变增大，峰值强度降低。Jia et al. (2017)发现颗粒破碎改变了堆石料的结构，并且对其剪胀性、摩擦角、强度、渗透性以及在地震荷载下的变形有较大的影响。陈火东等(2018)发现颗粒破碎使钙质砂体应变从剪胀逐渐发展到剪缩，且使得钙质砂峰值摩擦角减小，抗剪强度降低。由于颗粒破碎对材料的强度和变形有很大的影响。徐永福(2018)采用离散单元软件PFC2D模拟了考虑颗粒破碎影响的粗粒土的直剪试验，发现颗粒破碎对剪切强度的破碎准则有影响，颗粒破碎试样的强度包络线是幂函数关系。因此使用本构模型来模拟易碎材料的行为时，颗粒破碎是一个不可忽略的因素。

目前，将颗粒破碎考虑进本构模型的方法主要有以下几种：(1)通过考虑颗粒破碎来修正塑性硬化和剪胀关系。如：Sun et al. (2007)基于剑桥模型，考虑了颗粒破碎对等向压缩变形、应力-剪胀方程和剪切强度特性的影响，提出了一种用于颗粒材料的弹塑性本构模型。(2)基于损伤力学引入损伤因子。如：汪稔等(2002)采用边界面模型，在弹性应力-应变关系上引入颗粒破碎损伤因子，提出了一个损伤-滑移耦合作用力学模型。(3)考虑颗粒破碎对临界状态线的影响。如：Liu et al. (2013)基于临界状态土力学，通过采用临界状态线漂移的方法来描述颗粒破碎及其影响，并通过双曲函数使颗粒破碎与消耗的塑性功相关，提出一个考虑颗粒破碎的弹塑性本构模型。

现有的考虑颗粒破碎的本构模型大多采用隐式方法，通过将颗粒破碎与流动规则，屈服函数，损伤指数，临界状态线的形状或模型参数相结合来建立本构模型。只有很少的研究者在他们的模型中直接引入颗粒破碎，将模型行为与破碎指数，限制围压或输入能量联系起来。此外，用于钙质砂的考虑颗粒破碎的本构模型也较少。

在众多本构模型中，邓肯-张模型具有参数确定方便，使用简单等优点，故本文基于邓肯-张模型进行颗粒破碎修正。本文首先对某种钙质砂(钙质砂1)进行了一系列不同围压下的常规三轴排水试验，然后通过试验数据分析了相对破碎Br对邓肯-张模型参数的影响规律，将颗粒破碎直接引入邓肯-张模型中，建立了一个考虑颗粒破碎的钙质砂修正邓肯-张模型。为验证修正后模型的准确性，最后本文还使用该模型对文献中其他3种钙质砂(钙质砂2、钙质砂3以及钙质砂4，具体参数见表1)的常规三轴排水行为进行了模拟，结果表明拟合较好，并在颗粒破碎较大的情况下明显优于传统邓肯-张模型。

表1 钙质砂基本性质Table 1 The basic properties of calcareous sands

图1 钙质砂1 ε1/(σ1-σ3)-ε1曲线Fig. 1 The ε1/(σ1-σ3)-ε1 relations of calcareous sand No.1

1 考虑颗粒破碎的邓肯-张E-B模型

当邓肯-张E-B模型直接应用于钙质砂等易碎材料时，由于颗粒破碎的存在，模型将存在诸多问题。例如，邓肯-张模型假设应力-应变关系为双曲线，因此试验结果整理成形式ε1/(σ1-σ3)-ε1将呈线性关系，并以此来确定切线模量Et所需的参数。但当颗粒破碎存在时，ε1/(σ1-σ3)-ε1将不为线性，图1为钙质砂1在不同围压下的ε1/(σ1-σ3)-ε1曲线。而且当颗粒破碎较大时，Ei、B与围压的关系不再为Knpa(σ3/pa)n的形式，若依旧使用该形式进行拟合，不仅拟合度很低而且n将小于0。因此，要将邓肯-张模型应用于钙质砂，必须要对其进行颗粒破碎方面的修正。

1.1 颗粒破碎的度量

为了准确度量颗粒破碎的大小，大量学者进行了研究并提出了各种颗粒破碎指标。Lee et al. (1967)使用B15=D15i/D15f(D15为级配曲线上重量百分比为15%时所对应的颗粒粒径大小，其中下标i、f对应试验前后的级配曲线)来描述颗粒破碎的大小。Marsal(1967)使用Bm来表示描述颗粒破碎大小(Bm=∑|Bmai-Bmbi|，其中，Bmbi、Bmai分别表示试验前后某一粒组含量)。Lade et al. (1996)基于有效粒径D10提出了破碎指标B10的概念(B10=1-D10f/D10i，D10为重量百分比为10%时的颗粒粒径大小，其中下标i、f对应试验前后的级配曲线)。Hardin(1985)采用相对破碎Br=Bt/Bp来度量颗粒破碎的大小。Bp为破碎潜能(即0.074imm的粒径线与初始级配曲线与之间的面积)；Bt为总破碎量(即试验前后两级配曲线之间的面积)。

其中Hardin提出的相对破碎Br这一指标很好地克服了单独考虑某一粒径或含量的缺点，对于不同材料在破碎过程中抓住了整体级配的变化。因此，在本文使用Br来度量颗粒破碎的大小。

1.2 E-B模型

邓肯-张E-B模型建立在广义虎克定律的弹性理论基础上，因为其可以反映土的非线性，加之其参数不多且易确定，故其成为土的最为普及的本构模型之一。

常规三轴排水试验中，模型描述的应力-应变关系为(Duncan et al.，1970)：

(1)

式中，σ1、σ3为最大、小主应力；a、b为试验常数；ε1为轴向应变。

在常规三轴压缩实验中，切线模量Et表达式为：

(2)

对于粗粒土，普遍认为黏聚力c为0，内摩擦角φ随围压的增加而减小(Duncan et al.，1980)，表达式为：

φ=φ0-Δφlg(σ3/pa)

(3)

式中，φ0与Δφ为内摩擦角相关的材料常数。

E-B模型中使用体积模量B来近似描述土的体变特性。体积模量B的定义为：

(4)

式中， (σ1-σ3)70%与(εv)70%为(σ1-σ3)达到70%qf时的偏应力和相应的体应变的试验值。体积模量B与围压的关系为：

(5)

式中，Km、m为与体积模量相关的材料常数。

整个E-B模型共有7个参数，分别为：Kn、n、Rf、Km、m以及φ、c(或φ0、Δφ)。

1.3 修正E-B模型

基于邓肯张E-B模型，采用相对破碎Br度量颗粒破碎大小，综合考虑颗粒破碎对其相关参数的影响，本文提出了一个考虑颗粒破碎的钙质砂修正E-B模型。具体如下：

由于E-B模型中Ei由ε1/(σ1-σ3)-ε1直线的截距确定，但颗粒破碎的存在会导致其不为直线，故本文模型中采用偏应力为50%qf时的割线模量E50代替Ei进行分析。E50与Ei的关系为：

(6)

张家铭等(2008)研究发现颗粒破碎将导致钙质砂的强度降低，压缩性增加。Yu(2017)对珊瑚砂进行了一系列三轴排水剪切试验后发现颗粒破碎导致峰值状态摩擦角减小。因此，假定相同围压下E50、φ与B随Br的增加而减小。本文通过对钙质砂1进行一系列常规三轴排水试验并对E50、φ与B进行分析后发现，E50、φ、B与σ3及Br的关系为：

(7)

φ=φ0-Δφ1lg(σ3/pa)-Δφ2Br

(8)

(9)

钙质砂1采用式(7)～式(9)计算的E50、φ与B与通过试验结果计算的E50、φ与B对比如图2所示。可以发现，拟合较好。考察公式(7)～公式(9)可以发现，当无颗粒破碎时(Br=0)，模型退化为初始的E-B模型。

对于破坏比Rf，由于颗粒破碎对其影响不大。为简化模型，故不考虑颗粒破碎对破坏比Rf的影响。

同时，由于无法直接确定各应力状态下的相对破碎Br，因此Br必须使用另一状态量来描述。有研究表明颗粒破碎与输入土体的能量密切相关(Lade et al.，1996)，而输入能量可以通过应力及应变计算。故假定相对破碎Br是输入能量E函数。本文通过对钙质砂1进行一系列不同试验条件的三轴试验和侧限试验后发现，钙质砂1相对破碎Br与输入能量E的关系如图3所示。从图3中可以发现，当输入能量E较小时Br与E的关系为双曲线关系，较大时为指数关系，且这种关系受实验条件影响不大。由于双曲线关系只在输入能量很小时存在，为简化模型，假定相对破碎Br与输入能量E呈指数关系并认为钙质砂初始存在一定的破碎。故钙质砂相对破碎Br与输入能量E的关系表达式为：

Br=abE

(10)

式中，a、b为描述土体颗粒破碎的模型参数。

在常规三轴排水试验中每个计算点的输入能量E与应力及应变的关系为：

(11)

图2 钙质砂1 E50、摩擦角、B与围压及Br关系图Fig. 2 The E50、friction angle and B-confining pressure and Br relationships of calcareous sand No.1

图3 钙质砂1相对破碎Br与输入能量E关系曲线Fig. 3 The Br-E relationship curve of calcareous sand No.1

式中， dε1， dεv为轴向应变增量与体应变增量；SOT为试验开始点；POT为当前计算点。

最终，切线模量Et按式(12)计算，体积模量B按式(9)计算。

(12)

2 模型参数的确定方法

本文提出的考虑颗粒破碎的钙质砂修正邓肯-张模型共有12个参数，分别为：与E50相关的3个常数Kn1、Kn2、n； 与内摩擦角相关的3个常数φ0、Δφ1、Δφ2； 与体积模量相关的3个常数Km1、Km2、m； 用于描述颗粒破碎的两个常数a、b； 破坏比Rf。模型的常数大多具有明确的物理意义，可以通过传统的实验室试验进行校准。为保证Br-E拟合的准确性，需要进行至少一组围压较低的常规三轴排水试验。但围压较低时，钙质砂容易发生剪胀，剪胀将导致体积模量B偏大。故为保证B拟合的准确性还需进行3组围压较高不存在剪胀的常规三轴排水试验。故本模型参数确定总共需要进行至少4组常规三轴排水试验。后文将以钙质砂4(翁贻令， 2017)为例详细介绍模型参数的确定方法。钙质砂4进行的一组围压较低的常规三轴排水试验的围压为50ikPa，进行的3组围压较高的常规三轴排水试验的围压分别为100ikPa、200ikPa以及400ikPa。

2.1 颗粒破碎相关常数a、b

首先通过式(11)与试验前后级配曲线计算出试验结束时的输入能量E以及相对破碎Br。不同围压下钙质砂4试验结束时的(E，Br)分别为(44.￣85， 0.￣0372)、(74.￣47， 0.￣0536)、(124.￣47， 0.￣0876)及(190.￣62， 0.￣1129)。然后按式(10)对Br及E进行拟合即可得到参数a、b。对于钙质砂4，函数关系式为Br=0.￣033i65×1.￣006i55E，a、b分别为0.￣033i65、1.￣006i55。

2.2 内摩擦角相关常数φ0、Δφ1、Δφ2

首先按应力-应变曲线有峰值时取峰值，无峰值时取轴向应变15%时的应力点的原则，寻找破坏时的偏应力qf，再通过qf计算出内摩擦角。不同围压下钙质砂4内摩擦角分别为46.97°、43.52°、41.19°以及35.39°。然后按式(11)计算破坏时输入能量E，并通过E与Br的关系计算出破坏时的相对破碎Br。对于钙质砂4，Br分别为0.041i182、0.050i718、0.075i828以及0.116i782。最后将摩擦角φ作为应变量，破坏时的相对破碎Br与lg(σ3/pa)作为自变量按式(8)形式进行双变量线性拟合即可得到参数φ0、Δφ1、Δφ2。对于钙质砂4，φ0、Δφ1、Δφ2分别为48.02、6.04以及75.1。

2.3 E50相关常数Kn1、Kn2、n

首先确定各围压下50%qf对应的割线模量E50。对于钙质砂4，其分别为136ikPa、221ikPa、386ikPa以及550ikPa。然后按式(11)计算50%qf时的输入能量E，并通过E与Br的关系计算出50%qf时的相对破碎Br。对于钙质砂4，其分别为0.034i12、0.034i61、0.038i16以及0.045i57。最后按式(7)形式对E50、Br及σ3进行双变量非线性拟合即可得到参数Kn1、Kn2、n。对于钙质砂4，Kn1、Kn2、n分别为91.98、0.233以及264。

2.4 体积模量相关的3个常数Km1、Km2、m

首先通过试验数据找到偏应力为70%qf对应的偏应力以及体应变，按式(3)计算出各围压下的体积模量B。对于钙质砂4，由于围压50ikPa时发生剪胀，故只使用围压100ikPa、200ikPa以及400ikPa的数据来确定参数Km1、Km2、m。各围压下体积模量B分别为2991ikPa、3104ikPa以及2792ikPa。然后按式(11)计算70%qf时的输入能量E，并通过E与Br的关系计算出70%qf时的相对破碎Br。对于钙质砂4，Br分别为0.￣037i24、0.￣044i05以及0.￣059i26。最后按式(9)形式对B、Br及σ3进行双变量非线性拟合即可得到参数Km1、Km2、m。对于钙质砂4，Km1、Km2、m分别为29.￣98、0.￣079以及0.￣2211。

2.5 破坏比Rf

破坏比Rf的确定方法与邓肯-张模型相同。首先通过式(13)计算出极限偏应力qult。对于钙质砂4，其分别为311ikPa、546ikPa、1074ikPa以及1891ikPa。然后通过式(14)计算出每组实验的破坏比。对于钙质砂4，破坏比Rf分别为0.875、0.809、0.718以及0.582。最后对不同围压下的破坏比求平均得到破坏比Rf，钙质砂4的破坏比为0.746。

(13)

(14)

式中， 70%、95%为偏应力为70%qf及95%qf时的相关试验数据。

3 试验验证

本节中使用本文提出的修正邓肯-张模型对4种钙质砂的常规三轴排水行为进行了模拟。表1列出了4种钙质砂的基本性质以及数据来源。表2列出了4种钙质砂的模型参数。模型参数的确定基于上一节中阐述的方法进行。

表2 模型参数Table 2 Model constants

由于邓肯-张模型无法考虑应力软化现象，因此模拟只模拟到应力峰值点。另外，模拟结果还与原始邓肯-张E-B模型模拟结果进行了对比。由于颗粒破碎导致4种钙质砂的原始邓肯-张E-B模型的体积模量B不再为Kmpa(σ3/pa)m的形式，若依旧使用该形式进行拟合，不仅拟合度很低而且m将小于0，因此本文未给出修正前后模型的体应变模拟的对比。

3.1 钙质砂1

本文对钙质砂1进行了300～1000ikPa不同围压的常规三轴排水试验。 对于钙质砂1， 由于其应力没有峰值， 因此取15%轴向应变对应的数据点作为破坏点， 图4为钙质砂1的模拟结果。 从图中可以看到， 修正模型对钙质砂1的常规三轴排水行为模拟较好。 由于钙质砂1粒径较大、 围压较大， 颗粒破碎明显， 使得Ei与围压的关系不再为Knpa(σ3/pa)n的形式，若依旧使用该形式进行拟合，不仅拟合度很低而且n将小于0，故原始邓肯-张模型无法对钙质砂1的常规三轴排水行为进行模拟，因此本文也未给出修正模型与原始邓肯-张模型的对比结果。

图4 钙质砂1试验验证Fig. 4 Model verification with calcareous sand No.1

3.2 钙质砂2(张家铭， 2004)

张家铭(2004)对钙质砂2进行了100～3200ikPa不同围压的常规三轴排水试验。图5为钙质砂2的模拟结果。从图5a中可以看出，虽然原始邓肯-张模型的Ei仍为Knpa(σ3/pa)n的形式，但拟合效果随着颗粒破碎的增加逐渐下降，而使用本文提出的修正邓肯-张模型拟合度较好。确定钙质砂2的模型参数时，发现1600ikPa时摩擦角反而比800ikPa时大，因此推测1600ikPa时数据可能有误，故围压1600ikPa时应力拟合相差较大。

图5 钙质砂2试验验证Fig. 5 Model verification with calcareous sand No.2

从图5b中可以看出，在低围压下，体应变拟合效果不理想，这是由于低围压下钙质砂发生剪胀，而邓肯-张模型无法模拟剪胀导致的，而当围压大于400ikPa时，剪胀基本消失，体应变也拟合较好。

3.3 钙质砂3(王刚等， 2018)

王刚等(2018)对钙质砂3进行了200～1600ikPa不同围压的常规三轴排水试验。图6为钙质砂3的模拟结果。可以看出，使用本文提出的修正邓肯-张模型拟合度较好，但使用原始邓肯-张模型拟合的效果随着颗粒破碎的增加逐渐变差。

图6 钙质砂3试验验证Fig. 6 Model verification with calcareous sand No.3

3.4 钙质砂4(翁贻令， 2017)

翁贻令(2017)对钙质砂4进行了50～400ikPa不同围压的常规三轴排水试验。钙质砂4为单一粒径，只有2～5imm的颗粒，颗粒较大，故在试验围压不大时，破碎量依然较大。图7为钙质砂4的模拟结果。可以看到，本文提出的修正邓肯-张模型拟合度较好。由于Ei与围压的关系不为Knpa(σ3/pa)n的形式，若依旧使用该形式进行拟合，不仅拟合度很低而且n将小于0，故未给出修正模型与原始邓肯-张模型的对比结果。

图7 钙质砂4试验验证Fig. 7 Model verification with calcareous sand No.4

4 结论与讨论

本文通过考虑相对破碎Br对邓肯-张E-B模型参数的影响规律，建立了一个考虑颗粒破碎的钙质砂修正邓肯-张E-B模型。关于修正模型的相关结论如下：

(1)颗粒破碎将导致钙质砂的强度降低以及压缩性增加，反映到邓肯-张模型中为切线模量Et、内摩擦角φ和体积模量B随相对破碎Br的增加而减小。

(2)对于钙质砂，相对破碎Br与输入能量E呈指数关系，且其关系受试验条件的影响不大。输入能量E可以通过应力及应变计算得到，因此通过输入能量E可以将颗粒破碎引入邓肯-张模型中。

(3)本文使用修正邓肯-张模型对4种钙质砂的常规三轴排水行为进行了模拟，发现拟合效果较好，验证了修正模型的合理性。并且当颗粒破碎较大时，修正后模型明显优于原始邓肯-张模型。

同时，本文提出的修正模型基于颗粒破碎与输入能量呈指数关系。但实际颗粒破碎不会随输入能量无限增加，推测输入能量达到一定程度，颗粒破碎的增加趋势将趋于平稳，但由于仪器限制以及工程上一般不会超过本文研究的应力范围，故本文未对更高围压下钙质砂颗粒破碎规律做详细研究。因此，本文提出的考虑颗粒破碎的钙质砂修正邓肯-张模型可能不适用于模拟更高应力水平下钙质砂的行为，更高应力水平下钙质砂行为模拟有待进一步研究探讨。