盾构技术的发展为城市地铁隧道提供了安全、高效、环保的施工方法，然而近年来部分城市可供利用的中浅层地下空间日趋饱和，加之一些地区缺乏整体的地下空间开发规划，盾构掘进在城区地层中面临着愈发复杂的施工环境和层出不穷的障碍限制. 密集的房屋、桥梁等建筑物所附属的地下桩基、连续墙等混凝土结构常常成为盾构掘进路线上的重大障碍^[1]，导致盾构机不可避免地需要直接切削破除此类障碍物，迫使一些本用于切削土体的盾构刀具承担起混凝土切削的任务，如苏州轨道交通2号线切削广济桥桩基、广州地铁3号线切削居民楼桩基、上海轨道交通10号线切削沙泾港桥桩基等^[2-4]. 除此之外，盾构机在进、出洞时切除混凝土封门的情况则更为频繁，这给刮刀等软土刀具及刀座、盾构机刀盘与轴承等部件的强度，以及所需匹配的功率提出了更高要求.

现阶段国内外工程材料切削模型的研究成果主要集中在破岩掘进或土体耕挖方面. 1985年，Rånman等^[5]针对岩石材料的切削过程提出了局部压碎模型，认为切削时刀具和岩石间会积蓄弹性能，直至切削力达到岩石破坏强度后发生岩石局部压碎而释放弹性能，从而完成切削动作. Ucgul等^[6]在离散元软件DEM中采用了滞回弹簧接触模型，探讨了旋耕刀在砂壤土中的拉力与刀具底部长度、刀具前角、切入夹角之间的关系. 叶勇等^[7]指出岩石在受刀具切削时的破坏模式是由塑形破坏模式到过渡模式，再到脆性破坏模式，并将离散元模拟结果与实际试验进行了对比分析. 尽管盾构工程中刀具直接切削混凝土的情况越来越多，现阶段对混凝土受切削破坏过程的描述与刀具受力计算方法的研究仍十分有限. 王飞^[8]基于盾构直接切削桩基的工程进行了数值模型和现场切桩试验，探究了刀具破除混凝土的机理并提出了较为完整的盾构切桩施工技术体系. 周里群等^[9]在PFC 3D中通过接触黏结模型模拟了沥青混凝土的切削过程，采用虚拟试验校核了模型参数，对切入角度等因素进行了正交试验研究. 随着计算机技术的发展，建立数值模型进行数值仿真计算成为材料切削问题最为常用而且有效的研究方法之一. 混凝土在受切削过程中会出现挤压、破碎、崩落等现象，选用机理清晰、参数具体的本构模型准确描述这一系列现象是数值计算的关键.

本文引入HJC混凝土动态本构模型对混凝土材料受盾构刀具切削破坏过程进行数值计算，阐述HJC模型参数的确定方法，并根据混凝土试块的室内切削试验结果对相关参数进行对比修正，进而为预估刀具刀盘在切削混凝土材料时的负载水平、确定盾构破除混凝土障碍施工的可行性与切削混凝土时合理选择施工参数提供一定的数值计算方法和参考依据.

Holmquist-Johnson-Cook 混凝土动态本构模型(HJC模型)^[10]是针对混凝土类材料提出的一种综合考虑应变率效应、损伤演化效应、围压效应和压碎、压实效应影响的本构模型，较好描述了混凝土类材料在大变形、高应变率和高静水压力下的力学行为，且形式简单，参数物理意义明确，已被LS-DYNA有限元分析软件引入，被广泛应用于冲击爆炸等强动载作用下混凝土类材料的动态响应分析中^[11]. 盾构刀具切削混凝土过程同样是材料受动态侵彻而发生损伤、压碎、破坏的过程，因此可以考虑将HJC模型用于盾构刀具切削混凝土类材料的数值计算中. 该模型包含强度模型、损伤演化模型与状态方程3个部分（见图1、2、3）.

如图1所示，在HJC模型中，材料的强度模型采用归一化的等效应力来描述：

(1) ${\sigma _{{\rm{eq}}}}^{\rm{*}} = [A(1 - D) + B{({p^*})^N}](1 + C\ln {\dot \varepsilon ^*}).$

式中： $\!{\sigma _{{\rm{eq}}}}^{\rm{*}}{\!\rm{=}}{\sigma _{{\rm{eq}}}}\!/\!{f_{\rm{c}}}\;$为归一化等效应力，且 ${\sigma _{{\rm{eq}}}}^{\rm{*}}\!\!\leqslant\!\!S{_{{\rm{FMAX}}}}$， ${S_{{\rm{FMAX}}}}$为材料最大归一化等效应力， ${f_{\rm{c}}}$为准静态单轴抗压强度， ${\sigma _{{\rm{eq}}}}$为实际等效应力；D为损伤变量，其值由损伤模型确定； ${p^*} = p/{f_{\rm{c}}}$为归一化静水压力，p为静水压力； ${\dot \varepsilon ^{\rm{*}}}$为材料应变率；C为应变率影响参数；A为归一化黏度强化系数；B为归一化压力硬化系数；N为压力硬化指数；A、B、C、N统称为模型参数.

在HJC模型中，材料的损伤由其塑性应变累积形成，而塑性应变包括剪切塑性应变 $\varepsilon _{\rm{p}}^{\rm{f}}$与体积塑性应变 $\;\mu _{\rm{p}}^{\rm{f}}$，如图2所示. 损伤变量表示为

(2) $D = \sum {\frac{{\Delta {\varepsilon _{\rm{p}}} + \Delta {\mu _{\rm{p}}}}}{{\varepsilon _{\rm{p}}^{\rm{f}} + \mu _{\rm{p}}^{\rm{f}}}}},$

(3) $\varepsilon _{\rm{p}}^{\rm{f}} + \mu _{\rm{p}}^{\rm{f}}{\rm{ = }}{D_1}{({p^*} + {T^*})^{{{{D}}_2}}} \geqslant {E_{\rm{FMIN}}}.$

式中： $\Delta {\varepsilon _{\rm{p}}}$、 $\Delta {\mu _{\rm{p}}}$分别为数值计算中一个积分步之内的等效塑性剪切应变增量与等效塑性体积应变增量； $\varepsilon _{\rm{p}}^{\rm{f}}{\rm{ + }}\mu _{\rm{p}}^{\rm{f}}$为材料在当前积分步下的塑性应变值； ${T^*}{\rm{ = }}T/{f_{\rm{c}}}$为材料归一化抗拉强度，T为材料抗拉强度； ${D_1}$、 ${D_2}$为损伤常数； ${E_{\rm{FMIN}}}$为材料破坏时的最小塑性应变.

HJC模型中体积应变 $\mu $与静水压力p之间的关系分3个阶段描述，分别为弹性阶段（OA）、压密阶段（AB）与密实变形阶段（BC），如图3所示.

1）弹性阶段， $\mu $与p之间呈线性关系，即

(4) $p = K\mu ;\;p \leqslant {p_{\rm{c}}},\;\mu \leqslant {\mu _{\rm{c}}}.$

式中：K为体积模量； ${p_{\rm{c}}}$为弹性极限静水压力； ${\mu _{\rm{c}}}$为弹性极限体积应变.

2）压密阶段，材料内部空隙逐渐被压缩而产生塑性体积变形，此时， $\mu $与p之间的关系为

(5) $p = {p_{\rm{c}}} + {K_l}(\mu - {\mu _{\rm{c}}}); \;{p \leqslant {p_{\rm l}},\;\mu \leqslant {\mu _{\rm l}}} .$

(6) ${K_{\rm{l}}} = ({p_{\rm{l}}} - {p_{\rm{c}}})/({\mu _{\rm{l}}} - {\mu _{\rm{c}}}).$

式中： ${K_{\rm{l}}}$为材料压密阶段静水压力与体积应变的比例系数； ${p_{\rm{l}}}$为压密静水压力， ${\mu _{\rm{l}}}$为对应的体积应变.

3）密实变形阶段，材料被完全压碎，内部无空隙，此时 $\mu $与p之间的关系为

(7) $p = {K_1}\overline \mu + {K_2}{\overline \mu ^2} + {K_3}{\overline \mu ^3}.$

式中： ${K_1}$、 ${K_2}$、 ${K_3}$为压力常数；由于材料进入该阶段时可能发生软化现象而引入修正的体积应变 $\overline \mu $，其表达式为

(8) $\overline \mu {\rm{ = }}(\mu - {\mu _{\rm{l}}})/(1 + {\mu _{\rm{l}}}).$

在LS-DYNA显示动力学分析程序中，通过定义HJC关键字文件可以将材料参数赋予有限元模型，HJC模型的关键字文件选项卡如图4所示. 图中，MID为参数所赋予的材料编号， $\rho $为材料密度（图中用RO表示），G为材料剪切模量，ESP0为材料参考应变率. $\rho $、G、T、f_c等混凝土材料基本力学性质按照试验标准^[12]测得.

在不考虑材料损伤的准静态常规三轴试验条件下（D=0， ${\dot \varepsilon ^{\rm{*}}}{\rm{ = }}{10^{{\rm{ - }}5}}\sim {10^{{\rm{ - }}4}}$），假设材料强度满足Mises屈服准则，偏平面上HJC模型屈服面方程简化为

(9) ${\sigma _{{\rm{eq}}}}^* = A + B{({p^*})^N}.$

式中： ${\sigma _{{\rm{eq}}}}^*{\rm{ = (}}{\sigma _1} - {\sigma _3})/{f_{\rm{c}}}$， ${({p^*})^N} = ({\sigma _1} + {\sigma _2} + {\sigma _3})/(3{f_{\rm{c}}})$； ${\sigma _1}$为轴压， ${\sigma _2}{\rm{ = }}{\sigma _3}$为围压，该屈服面在偏平面上的投影为Mises圆.

根据塑性力学相关理论，当罗德角 $\theta {\rm{ = }}{60^ \circ }$时，偏平面上Mises屈服面的压缩子午线与Mohr-Coulomb(M-C)准则包络线相交于纯剪切和单轴压缩2个特征点，如图5所示.

屈服准则的表达式为

(10) ${\sigma _{{\rm{eq}}}}{\rm{ = }}c + p\tan \varphi .$

式中：c为黏聚力， $\varphi $为内摩擦角，两参数可通过不同围压下的混凝土三轴试验测得. 若Mohr-Coulomb屈服准则满足线性关系，则式（10）经过特征点（0， ${\tau _0}$），即

(11) ${\tau _0}{\rm{ = }}c + 0 \cdot \tan \varphi .$

式中： ${\tau _0}$为混凝土纯剪切状态下的抗剪强度. 同理，式（9）经过特征点（0， ${\tau _0}$）时，有

(12) ${\tau _0}/{f_{\rm{c}}} = A + 0 \cdot B.$

在三轴试验中，Mohr-Coulomb准则中的黏聚力c可根据 ${\sigma _1}$与 ${\sigma _3}$的线性关系斜率k计算^[13]：

(13) $c = {f_{\rm{c}}}/(2\sqrt k ).$

联立式（11）~（13），可得归一化的黏度强化系数A的表达式为

(14) $A = 1\big/\left(2\sqrt k\right). $

因此，可通过不同围压下混凝土圆柱体试件的三轴试验得到轴压 ${\sigma _1}$与围压 ${\sigma _3}$的线性关系斜率k，从而根据式（14）计算得到A值，同时测得不同静水压力p及其所对应的轴压与围压之差 ${\sigma _{{\rm{eq}}}}$，进而归一化得到多组坐标（ $\sigma _{{\rm{eq}}}^*$， ${p^*}$），再根据式（9）对B、N的值进行线性拟合求解. 试验数据及拟合结果如图6所示.

混凝土压缩第一阶段的状态方程参数通常采用三轴试验测定，一般有 ${p_{\rm{c}}} = {f_{\rm{c}}}/3$与 ${\mu _{\rm{c}}} = {p_{\rm{c}}}/K$^[10]；对于第二、三阶段，通常采用Hugoniot试验或经验公式确定，由于实验成本较高，一般可借鉴相关文献所述的实验数据对K₁、K₂、K₃进行拟合. 本文采用王永刚等^[14-17]开展的混凝土Hugoniot试验得到的数据进行拟合，如图7所示. 根据图7，同时可确定状态方程中不同阶段分界点的静水压力 ${p_{\rm{c}}}$、 ${p_{\rm{l}}}$及其对应的体应变 ${\mu _{\rm{c}}}$、 ${\mu _{\rm{l}}}$.

混凝土材料的抗压强度随应变率的增大而增大，为了确定应变率对抗压强度的影响系数C，一般基于混凝土的SHPB实验测得归一化等效强度 $\sigma _{{\rm{eq}}}^*$随应变率对数 $\ln {\dot \varepsilon ^{\rm{*}}}$的变化关系，然后进行线性拟合. 但由于设备条件不同，试验结果较为离散，一般C值的拟合结果在0~0.02^[18]. 本文对相关文献结果^[19-23]进行整理，采用最小二乘法进行线性拟合，如图8所示.

HJC模型假设材料的损伤参数与混凝土强度无关，因此损伤参数按照文献[10]中所述方法取值，取 ${D_1} = 0.04$， ${D_2} = 1.0$， ${E_{\rm{FMIN}}} = 0.01$.

为了验证HJC模型在盾构刀具切削混凝土数值仿真计算中的适用性与准确性，本文基于北京交通大学城市地下工程教育部重点实验室的盾构刀具切削试验平台进行混凝土试块切削试验，如图9（a）所示. 平台由反力装置、动力装置、控制系统和数据采集系统4部分组成，试验所用切刀长度为116.9 mm，宽度为90 mm，侧面轮廓为38.7~70.0 mm的梯形，刀刃前、后倾角分别为15°与10°，而被切削混凝土试件模型为棱长为200 mm的立方体.

为了尽量避免由切削路径中的非均匀介质带来的切削阻力偶然性，所用混凝土为无粗骨料混凝土，即水泥砂浆. 按配合比^[24]制作正方体与圆柱体混凝土试件，养护28 d后分别开展基本力学性能指标测试（测试结果见表1）与混凝土切削试验. 根据第2章内容，可初步确定混凝土基本力学参数与图4中HJC模型关键字选项卡中各参数取值，汇总结果如表1所示.

切削试验初步设定刀具贯入度H_c=10 mm，切削基准速度v_c=13.3 mm/s，采集并记录2个液压杆中压力值，进而计算得到刀具与试块之间的接触面法向切削阻力F_n与切向切削阻力F_τ，两者与总切削阻力F的关系如图10所示.

当刀具刚接触试件时，两者发生冲击碰撞；随着刀具的继续切削，刀具与试件之间的接触关系为“接触、挤压、压碎、剥落”的不停循环，直至切削结束，整体过程表现为刀具在小幅跳动中完成切削动作；当切削将要结束时，由于试件下部有临空面，容易出现刀具将试件整块剥落的现象；试件被切削的痕迹如图9（b）所示，作为数值计算的对比数据，绘制接触面法向切削阻力随时间的变化规律如图13所示.

采用ANSYS LS-DYNA前处理工具建立盾构切刀有限元模型与被切削混凝土试件有限元模型，合理调整被切削区域与其他区域的网格比例后单元数量至8.91×10⁶. 由于盾构切刀的主体材料为高强度钢，在短暂切削过程中相对混凝土的质量损失与变形可以忽略，同时为了兼顾计算效率与精度，设置切刀模型为刚体. 整体有限元模型如图11所示.

模型中混凝土材料性质符合HJC模型，按照表1中参数设置HJC关键字文件，刀具贯入度、切削速度与试验保持一致，采用LS-DYNA求解器对模型进行求解，同时输出接触面法向切削阻力. 求解完成后在LS-PrePost中查看土体模型被切削过程中的塑性应变区动态分布，如图12所示.

输出法向切削阻力F_n随时间t的变化情况，如图13（a）所示. 将其与混凝土试块切削试验中测得对应时间的法向切削阻力数据进行对比，并进行统计分别得到2组数据法向切削阻力的均值F_n^ave和标准差 $\sigma ({F_{\rm{n}}})$：

$ F_{\rm{n}}^{{\rm{ave}}}{\rm{ = }}\frac{1}{m}\sum\limits_{i = 1}^m {{F_{\rm{n}}}}, \;\;\sigma ({F_{\rm{n}}}){\rm{ = }}\sqrt {\frac{1}{m}\sum\limits_{i = 1}^m {({F_{\rm{n}}} - } F_{\rm{n}}^{{\rm{ave}}}{)^2}}. $

式中：m为F_n的采集样本数.

根据图13（a），混凝土试件的切削试验与基于HJC模型的数值计算所得到的法向切削阻力均随时间上下波动，其中试验结果的波动范围为15.1 ~ 24.9 kN（开始与结束时刻除外，下同），数值计算结果的波动范围较大，为8.3 ~ 32.1 kN；分析表明数值计算所得的切削阻力平均值较大（比试验结果大15%左右）. 另外，室内切削试验与数值模拟结果均表明，盾构刀具在切入混凝土表面时，切削阻力的波动幅度较大；而在切削至接近试块下方自由面时会出现剩余材料整块脱落、切削力骤降为0的现象，这一过程在数值模拟中相对平缓.

尽管在大小与波动幅度上均与切削试验所测数值存在一定的偏差，考虑到理想的数值模型与实际条件下混凝土试件切削之间在刀体材质、约束条件、接触关系等方面的差异，本文认为在通过相应试验初步确定各材料参数后，基于HJC模型的数值计算结果基本可以反映盾构刀具在切削混凝土材料过程中的切削阻力水平及其变化规律.

本文对状态方程与率效应中的参数选取采用文献试验数据拟合的方法，模型参数取值不能与切削试验中的材料性质完全吻合，这在一定程度上造成数值计算结果的偏差，因此，有必要对相关参数取值对数值计算结果的影响进行分析. 同时，在盾构施工中，切削速度由刀盘转速与刀具的安装半径决定；刀具切削深度（贯入度）则由盾构机的推进速度与刀盘转速之间的相对关系决定，即

(15) ${H_{\rm{c}}} = \frac{v}{{N\cos \alpha }}.$

式中：H_c为刀具贯入度，v为盾构机推进速度，N为盾构机刀盘转速，α为切刀前角。刀盘转速与推进速度均是盾构掘进中重要的可控参数，两者的适当调节对于控制盾构刀具切削阻力、保障盾构机顺利掘进有重要意义。

首先对率效应参数C值的影响进行分析，分别就C值为0.002、0.006、0.012与0.020时的数值模型进行计算，分别绘制计算所得法向切削阻力F_n随时间变化曲线，如图14所示. 对图14中数据进行统计分析，计算不同率效应参数C值时法向切削阻力的均值和标准差. 结果表明：率效应参数C主要影响切削阻力的标准差，即波动幅度，而对其平均值影响不大.

其次，由于损伤参数D₂对数值结果的影响极小^[17]，只对D₁对法向切削阻力的影响进行分析，取D₁值分别为0.01、0.02、0.04、0.08时的数值模型进行计算法向切削阻力，如图15所示. 同样统计不同D₁值时的F_n^ave与 $\sigma ({F_{\rm{n}}})$，结果表明：随着D₁的增大，切削阻力的波动幅度与平均值均有一定程度的增大.

根据室内混凝土切削试验数据与图14、15中数据的统计规律指标对表1中的率效应参数C与损伤参数D₁进行对比修正，修正后的C=0.006，D₁=0.02. 采用修正后的HJC模型参数重新进行数值计算，计算得到的法向切削阻力与室内切削试验结果对比如图13（b）所示. 经修正，基于HJC模型的法向切削阻力数值计算结果的平均值、波动幅度均与室内切削试验结果更加吻合.

保持刀具贯入度H_c=10 mm，采用参数修正后的HJC模型取切削基准速度v_c分别为13.3、19.5、26.60与33.3 mm/s计算法向切削阻力，如图16所示. 当切削速度增大时，刀具切削阻力明显增大，这是由于HJC模型中混凝土的压碎破坏与材料应变率有关. 计算法向切削阻力平均值F_n^ave随切削速度的变化关系，如图17所示.

根据HJC模型数值计算结果，盾构刀具切削阻力平均值随切削速度呈指数形式增加，因此在切削混凝土材料时，减小刀盘转速可以明显降低单把刀具的切削阻力及其带来的整体刀盘的扭矩水平，而且切削阻力会随着刀具的安装半径呈指数形式增大.

保持切削速度v_c=13.3 mm/s，分别取切削深度为5、10、15与20 mm计算法向切削阻力，如图18所示.

当切削深度增加时，不仅刀具切削阻力增大，而且其波动幅度也急剧增大，这是由于HJC模型可以反映混凝土材料破坏时切削阻力波动时峰值与谷值的叠加效应，即较大深度的切削过程相当于若干较小深度切削过程的叠加，与实际情况符合. 计算法向切削阻力平均值F_n^ave随切削深度的变化关系，如图19所示.该图可进一步说明法向切削阻力随切削深度线性增大.

根据式（15），当切削混凝土的深度随盾构机的推进速度增大时，不仅单把刀具的切削阻力及其对应的刀盘扭矩平均水平会线性增大，而且由此带来的扭矩峰值会成倍增加.

（1）通过混凝土性能力学试验与相关文献分析初步确定了模型参数，基于HJC模型的数值计算结果可以反映盾构刀具在切削混凝土材料过程中的切削阻力水平及其变化规律.

（2）室内试验与数值模拟结果均表明，刀具在切入混凝土表面时，切削阻力的波动幅度较大；在切削接近试块下方自由面时会出现剩余材料整块脱落、切削力骤降为0的现象，但这一过程在数值模拟中相对平缓；

（3）HJC模型中率效应参数C主要影响切削阻力的波动幅度；随着损伤参数D₁的增大，切削阻力的波动幅度与平均值均有一定程度的增大；可以根据混凝土切削试验对参数C与D₁进行对比修正；

（4）数值计算结果显示，盾构刀具切削阻力平均值随切削速度呈指数形式增加，随切削深度线性增加；HJC模型可以反映混凝土压碎破坏与材料应变率的关系及切削阻力随深度的线性叠加效应.