弥散裂缝模型在盾构管片非线性数值模拟中的适用性分析

谢家冲，王金昌，林志军，李俊，陈页开

（1.浙江大学 建筑工程学院，浙江 杭州 310058；2.中国电建集团华东勘测设计研究院有限公司，浙江 杭州 311122；3.华南理工大学 土木与交通学院，广东 广州 510640）

0 引言

盾构隧道是由预制混凝土管片以螺栓与防水垫层连接构成的复杂地下结构体，若将各个部件抽离开则可将各个管片视为独立工作的扁平曲梁。相比普通混凝土梁结构，盾构管片具有独特的结构形式，同时也被赋予更高要求的工作性能。针对单片管片的承载力与结构非线性响应，国内外众多学者基于实际工程采用的各类管片展开足尺试验研究，取得了一系列试验成果。Meng等[1]针对钢纤维混凝土管片开展承载力试验研究，结果表明，钢纤维的掺入能够延迟开裂的影响并优化应力分布，有效提高开裂荷载和管片韧性。Cheng等[2]开展了四点循环加载的管片模型试验，并基于试验结果验证了解析模型的有效性，据此能够有效预测管片受弯力学特性。Abbas等[3]分别针对普通盾构管片和钢纤维混凝土管片开孔后的力学特性展开足尺试验研究，结果显示，普通盾构管片拥有更高的极限荷载与开孔力学特性，而钢纤维混凝土管片则具有更高的开裂荷载和更稳定的后开裂阶段。齐明山和柳献[4]以上海地铁通用管片为例展开力学性能试验研究，结果表明，纤维混凝土混合配筋管片裂缝宽度达到0.20 mm时对应弯矩、屈服弯矩、极限弯矩均明显高于普通钢筋混凝土管片的结果。周海鹰[5]针对不同支承条件下和不同配筋的管片展开了试验研究，并证明了一般力学方法计算结果的一致性。

综上所述，已有针对单片管片的力学特性研究主要通过足尺模型试验方法，并对管片材料、约束和加固条件展开参数分析，获得了管片在特定条件下精确的非线性行为。试验结果的可贵之处在于为其他计算手段提供了验证基础，据此可以有效预测不同环境与荷载条件下的盾构管片行为方式。已有部分研究基于管片试验结果，建立了管片数值计算模型并开展参数分析：Zhou等[6]建立了三维有限元模型并赋予混凝土单元随动强化弹塑性增量本构，模拟结果与试验结果一致性较好。杨雨冰和谢雄耀[7]基于断裂力学的有限元方法对单个管片、纵向接头以及整环衬砌结构3个层次讨论管片的破损机制，结果表明，单个管片承载力性能与约束条件关系较大。综上所述，数值模型是预测结构行为的有效手段，然而在针对管片非线性的研究中大多集中于管片承载力与挠度特点，忽略了裂缝扩展的特点与屈服后的特定行为方式，因此适用于管片破坏全过程的非线性分析数值手段与裂缝模型亟待研究。不同于钢筋较固定的行为方式，混凝土的开裂非线性行为更复杂，因此采用何种混凝土模型来准确地描述管片的开裂行为是值得关注的问题。

考虑开裂的混凝土模型大致可分为离散裂缝模型与弥散裂缝模型2种，最早分别由Ngo等[8]和Rashid[9]提出。其中弥散裂缝模型能够在连续介质体系下模拟结构内部各个方向的开裂，被广泛应用于混凝土结构非线性分析。本文基于Diana10.3有限元程序重点探讨3种常用弥散裂缝模型在盾构管片数值模拟中的适用性，分别为旋转裂缝模型、固定裂缝模型和多向固定裂缝模型。另外，针对断裂能展开参数分析，研究结果能够为不同尺度的盾构隧道结构分析提供借鉴。

1 弥散裂缝模型介绍

1.1 总应变裂缝模型

总应变裂缝模型通过应变的函数来描述应力变化，同时能够模拟受压、受拉和受剪条件下混凝土的非线性力学行为。该模型基本思想是将弹性应变εe与开裂应变εcr合二为一考虑，并基于裂缝方向计算评估应力大小。

根据裂缝方向是否随着主应力方向旋转可将总应变裂缝模型分为旋转裂缝模型和固定裂缝模型，图1给出2种模型裂缝应变和方向与主应力σ1和σ2之间的关系。

图1 总应变裂缝模型

从结构开裂的物理力学本质来讲，固定裂缝模型思想更接近实际。而从弥散裂缝模型的使用过程来看，旋转裂缝模型往往有更好的模拟效果而被广泛采用[10]。旋转裂缝模型中应变转换矩阵T由单元主应力方向决定，即：

在固定裂缝模型中应变转换矩阵T与开裂坐标系nst均由初期裂缝给出，并不随主应力的变化而发生改变。随后根据开裂应力-应变关系，并在开裂坐标系nst和单元坐标系xyz中得到应力结果：

其中裂缝的应力-应变关系可通过刚度矩阵D表达，即：

割线刚度矩阵D由材料刚度矩阵Dsecant得到：

1.2 多向固定裂缝模型

多向固定裂缝模型基于decomposed crack model（裂缝分解模型）思想，将应变分解为弹性应变εe与开裂应变εcr，如式（6）所示：

多向固定裂缝模型允许在一个积分点产生多条裂缝，每个积分点的裂缝都有与之相对应的应变分量eicr和应力分量sicr[10]。裂缝的数量则通过临界角α控制，新裂缝萌生条件必须满足主拉应力超过混凝土抗拉强度且当前主拉应力方向与既有裂缝法向的夹角超过α。模型局部坐标与临界角α示意如图2所示。

图2 多向固定裂缝模型

当裂缝的应力向量与应变向量存在式（7）关系时：

则每条裂缝的应力-应变关系可以解耦地表达如下：

因此，在多向固定裂缝模型中，每条裂缝的开裂应力仅与对应的开裂应变相关，简化了不同裂缝间复杂的耦合影响。与总应变裂缝模型相同，多向固定裂缝模型也遵守弥散开裂的断裂能法则，考虑了拉伸截断、受拉软化和剪切保留的非线性行为。其裂缝旋转方向与固定裂缝模型的一致，即开裂坐标系nst均由初期裂缝决定，不随主应力变化而改变。

2 基于试验的模型有效性验证

2.1 试验工况

本研究针对单片管片的数值模拟基于周海鹰[5]的试验成果进行。标准环尺寸和配筋如图3所示。管片为中心角为67.5°的标准环，环宽1.2 m、内径2.7 m、外径3.0 m。混凝土强度等级为C50，内外侧分别布置一排纵筋，保护层厚度分别为59、43 mm。

图3 标准环尺寸及配筋示意

试验采用三点加载，加载点位置如图4所示，加载方向指向隧道圆心，同时在加载点位置安置橡胶垫片以防止混凝土局部压碎。试验设置2种支承条件Pzw1和Pzw2，分别为在管片两端设置简支梁支承和轴力支承的方式。基于该试验获得的跨中挠度曲线、材料应变曲线和破坏模式，对混凝土裂缝模型的可靠性进行研究。

图4 加载和支承方式

2.2 数值模型

单个管片数值模型基于Diana10.3程序建立。模型为二维数值模型，混凝土单元和橡胶垫片均为CQ16M八节点四边形实体单元。针对混凝土非线性分析的精度要求，初步设置单元尺寸为0.02 m，即在梁高度方向划分15个单元，数值模型和网格划分如图5所示。试验工况Pzw1为管片两端底部设置铰支座，Pzw2支承条件为约束管片两端截面切向位移并设置轴力支承。非线性分析荷载步大小为10 kN，荷载步上限为200步。

图5 数值模型和网格划分

本节中混凝土采用旋转裂缝模型，混凝土材料参数如表1所示。钢筋模型采用嵌入式钢筋单元模拟，该钢筋单元不参与混凝土母体单元的网格划分，通过与母体单元结合的方式贡献刚度。在二维模型中，同一平面内的纵筋和箍筋均被简化为连续均质的钢筋网片，截面积之和作为网片单元几何的输入值。根据图3的配筋计算得到上、下侧纵筋和每环箍筋横截面面积之和分别为2035、1608.5、226.2 mm2。钢筋材料采用Von Mises塑性硬化模型，钢筋参数如表2所示。

表1 混凝土材料参数

表2 钢筋参数

2.3 模型验证

试验记录了开裂荷载、屈服荷载和破坏荷载，为了与试验结果对应，在数值模型中这3个关键荷载的判定标准分别为：出现宏观裂缝（最大裂缝宽度超过0.05 mm）、线性段结束点、荷载值最大点。管片的跨中挠度曲线如图6所示。

图6 管片的跨中挠度曲线

由图6可知，支承条件对管片的承载力影响较大，Pzw2和Pzw1工况中管片分别在荷载为960 kN和160 kN达到屈服状态，在破坏前均有明显的水平直线段。Pzw2工况的屈服点更明显，计算与试验结果均在跨中挠度为20 mm左右达到屈服，但计算结果给出的屈服荷载稍小。

表3针对3个关键荷载给出了试验值与计算值对比。

表3 关键荷载试验值与计算值对比

由表3可见，数值计算得到的各项特征荷载与实测结果吻合良好，其中开裂荷载计算值与试验值有一定差异，分析原因在于裂缝萌生初期难以观测，第1条观测得到的宏观裂缝具有一定随机性。综上所述，基于管片承载能力角度，2个工况的计算值都能够反映各自支承条件的管片位移响应特性。

图7以下排钢筋跨中位置的应变结果为例说明钢筋单元的可靠性。

图7 下排钢筋应变曲线

由图7可知，在2个工况中，数值计算结果都能与试验曲线吻合，屈服应变拐点的出现位置一致。Von Mises塑性硬化模型中将钢筋按多段直线简化，而在屈服段，Pzw2工况出现了明显的非线性段，计算应变结果偏大，但在曲线特征点能够与试验值吻合良好，钢筋单元到达屈服时的荷载P与实际一致。因此在二维数值模型中，嵌入式钢筋网单元能够以简化的形式有效模拟实际钢筋工作行为。

图8给出了2个工况管片各自达到极限荷载值时裂缝宽度云图与试验结果的对比。

图8 管片破坏特征对比

图8结果表明，支承条件同时极大程度地影响管片裂缝的开展模式。从跨中裂缝形态来看，Pzw2工况中跨中两侧的裂缝呈拱形分布，而Pzw1工况的裂缝相互平行生长，同时裂缝间距更大。从裂缝分布的范围来看，到达极限荷载时，在Pzw2工况中主要宏观裂缝集中于跨中区域，同时顶部混凝土被压溃，损伤区域扩展至拱顶外侧；而在Pzw1简支工况中，沿曲梁的长度方向均有宏观裂缝萌生，同时损伤深度更大，裂缝发展至外侧钢筋位置。Pzw2工况和Pzw1工况到达极限荷载时的最大裂缝宽度也差异较大，分别为5.77、9.42 mm。

根据与实测裂缝分布模式的对比可知，2个工况的计算结果与试验结果基本一致，能够较好地模拟各自到达极限荷载时的裂缝分布与损伤形态。Pzw2工况的破坏以跨中拱顶被压溃为标准，损伤开裂较严重区域集中于跨中位置，整体损伤程度要小于Pzw1工况的结果，Pzw2工况中轴力支承条件也更接近于实际管片工作状态。

3 计算结果与讨论

本节基于工况Pzw2试验结果进行参数分析，3种弥散裂缝模型参数如表4所示。

表4 弥散裂缝模型参数

表4中β为剪切保留系数，当β为常数时，剪切模量折减公式如式（9）所示；当剪切模量基于损伤计算时，即考虑开裂状态下泊松比的折减[11]，其公式如式（10）所示。

式中：Gcr——折减的剪切模量，MPa；

G——未考虑折减的剪切模量，MPa；

v'——考虑泊松比效应折减后的泊松比。

3.1 弥散裂缝模型对位移响应的影响

图9和表5分别给出采用不同裂缝模型下管片跨中的挠度曲线以及破坏特征值。

图9 不同弥散裂缝模型对管片跨中挠度曲线的影响

表5 破坏特征值试验与计算结果对比

由图9和表5可知，在加载前期，各种裂缝模型下的挠度曲线基本重合，当挠度到达20 mm时，各条曲线走向产生分歧。旋转裂缝模型的极限荷载最小；3条固定裂缝模型的跨中挠度曲线能达到更大的极限荷载，分别为1048、1050、1066 kN，更接近于试验值，但其达到极限荷载时对应的跨中挠度分别为23.8、24.8、24.3 mm，大于试验值的21.0 mm；采用多向固定裂缝模型模拟的管片在计算过程中未发生明显屈服，跨中挠度曲线基本呈线性增长，分析原因是其忽略了混凝土的开裂软化性能，跨中顶部未达到压溃状态。

值得注意的是固定裂缝模型中对剪切保留系数β的参数分析，当β越大，即开裂后的剪切模量越大，当时挠度曲线达到屈服荷载后表现为有较长的水平屈服段，与试验曲线接近。当β系数折减基于损伤即与泊松比效应关联后，曲线屈服段消失，在达到极限荷载后承载力迅速降低。

3.2 弥散裂缝模型对裂缝形态的影响

图10为应用不同弥散裂缝模型对裂缝宽度云图的影响。

图10 不同裂缝模型对开裂云图的影响

由图10可见，与旋转裂缝模型跨中呈“拱形”的裂缝相比，固定裂缝模型计算结果云图中裂缝扩展方向呈径向，同时在跨中部位呈现密集的裂缝。当剪力保留为基于损伤计算时，跨中出现一条主裂缝；当β=0.1时，跨中的裂缝密集且无规则；当β增加至0.5时，跨中的裂缝明显减少且呈清晰的竖向裂缝，同时最大裂缝宽度由11.5 mm降至10.6 mm；多向固定裂缝模型结果的裂缝形态则与固定裂缝模型结果类似，但跨中顶部未出现开裂。

由表5中宏观裂缝间距结果可知，多向固定裂缝模型计算获得的裂缝间距最小，仅为60.0 mm；旋转裂缝模型的裂缝间距最大，为130 mm，与试验结果较一致；固定裂缝模型的结果介于两者之间，且β的改变基本不影响裂缝间距的大小。因此综合位移响应和裂缝形态两个角度，旋转裂缝模型的计算能够更好低与试验结果吻合。

3.3 断裂能的影响

在弥散裂缝模型中断裂能是决定混凝土开裂行为最重要的因素之一，被定义为一条裂缝开裂全过程中单位面积所耗散的能量，即混凝土单轴拉应力σ与裂缝宽度wcr的积分。断裂能可分为拉伸断裂能GfI和压缩断裂能Gc，分别在拉伸软化与受压软化中控制应力-应变曲线。作为混凝土材料的一种固有属性，一般通过试验获得参考的取值范围，但取值具有较大波动性。同时各规范针对断裂能的计算方法有所差异。以欧洲规范Fib Model Code 2010[12]为例，其拉伸断裂能与压缩断裂能经验公式如式（11）、（12）所示。

式中：fcm——混凝土平均抗压强度，fcm=fck+8，N/mm2；

fck——混凝土圆柱体抗压强度，N/mm2。

立方体与圆柱体之间换算系数取0.84[13]。

另外，针对Gc和GfI的倍数关系，Feenstra[11]进行了大量试验，结果给出该倍数关系在50~100波动。综上可知，在计算中断裂能的取值存在较大主观性，因此本节对管片混凝土的两种断裂能展开参数敏感性分析。

图11为GfI和Gc对管片极限荷载的影响。

图11 断裂能对管片极限荷载的影响

由图11可知，当Gc=22 500 N/m时，随着GfI的增加，极限荷载基本不变，GfI从140 N/m提升至160 N/m时，极限荷载仅提高了10.6 kN。随着Gc的增大，GfI对极限承载力的影响显著提升，当Gc=37 500 N/m时，GfI从140 N/m提升至160 N/m，极限荷载增大了60.3 kN。若以GfI作为不变量，也可得到相似的结论，因此以管片承载力指标分析，GfI和Gc对管片承载力计算结果都具有较高敏感性同时其敏感性又与彼此的大小呈正相关。

基于欧洲规范Fib Model Code 2010给出的2种断裂能的计算公式与试件尺寸的换算系数，混凝土强度为C50时，GfI=147.6 N/m，对应的Gc为36 900 N/m，通过内插法计算得到此时管片的极限承载力约为1033.8 kN，与试验获得的1071 kN基本吻合。因此基于式（11）、式（12）的断裂能计算方法，获得的参数能够有效应用于钢筋混凝土管片的非线性分析。

4 结论

（1）分别对2种支承条件下的管片承载力试验进行数值模拟分析，通过跨中挠度曲线、钢筋应变曲线和裂缝形态等方面比较分析计算与试验结果，验证了该模型的可靠性，同时说明约束条件对管片承载力特征值、破坏模式有较大影响。

（2）在工况Pzw2数值模型基础上，研究不同混凝土弥散裂缝模型带来的影响。计算结果表明：3种裂缝模型下挠度曲线在各自达到屈服前基本重合；旋转裂缝模型模拟效果较优，跨中位移响应、裂缝间距和跨中“拱形”裂缝形态均与实际吻合良好；固定裂缝模型结果裂缝形态以径向裂缝为主，且在跨中部位更密集，剪力保留系数β对屈服后挠度曲线走向、跨中裂缝分布影响较大，但基本不影响极限承载力大小；多向固定裂缝模型在加载过程中未发生明显屈服，挠度曲线基本呈线性发展。当分析考虑屈服后行为和裂缝形态，则建议采用旋转裂缝模型，若计算仅要求结构屈服前行为并不要求裂缝分布的准确性，则3种模型均可采用。

（3）通过双因子参数分析可知，GfI和Gc对管片承载力计算结果都具有较高敏感性，同时其敏感性又与彼此的大小呈正相关，计算结果同时验证了欧洲规范给出的2种断裂能计算方法的适用性。