不同胶结程度MICP固化珊瑚砂的无侧限压缩离散元分析

2022-05-13章懿涛方祥位胡丰慧姚志华吴焕然申春妮

土木与环境工程学报 2022年4期

章懿涛，方祥位，胡丰慧，姚志华，吴焕然，申春妮

(1.重庆大学土木工程学院，重庆 400045；2.空军工程大学机场建筑工程系，西安 710038；3.重庆科技学院建筑工程学院，重庆 401331)

虽然目前针对微生物加固技术已开展了较多研究，但对在受力条件下固化试样微观结构的演化和细观介质的损伤、破坏机理认识很不充分，这些微观机理的研究对分析微生物固化土的宏观响应具有重要意义。一些微观结构分析试验方法，如SEM扫描电镜、GIS图形分析、MIP压汞试验及X射线衍射试验被用于微生物固化试样分析。然而，这些测试项目较难开展，耗时且费用昂贵，无法获取MICP固化材料在剪切作用下的胶结破坏机理。由Cundall等[6]提出的离散单元法(Discrete Element Method,DEM)是一种能够在颗粒尺度上建立物理过程并进行微、细观分析的数值方法。颗粒离散单元法采用非连续介质分析方法，通过颗粒介质构建计算模型，以力-位移方程和牛顿第二定律为理论基础，基于中心差分法对整个时间范围内颗粒的运动方程进行迭代求解，能够自动反映微观介质的连续非线性应力-应变关系，同时捕捉到介质破裂和微裂纹形成、扩展及相互作用的动态过程[7]。

离散元已被广泛用于模拟无黏性材料(例如砂土)和胶结材料(例如混凝土、胶结砂)，研究材料的物理力学性能及破坏机制[8-9]。Li等[10]采用三维离散元模拟了强、中、弱3种胶结强度的水泥胶结砂在三轴压缩下的力学行为，在黏结强度较强和中等的胶结试样中观察到了应变软化响应，基于离散元分析发现，在剪切开始时断键速率较为温和，随后逐渐增大，最终屈服形成剪切带，而弱胶结强度试样表现出应变硬化响应。徐东升等[11]采用离散元软件对不同级别的珊瑚砂水泥胶结体的破坏行为进行分析，发现随着砂颗粒级配范围及平均粒径的增大，微裂隙的发展会朝某个方向集中，试样裂缝数量减少，破坏程度减轻。

目前，已有一些学者针对MICP技术固化砂土进行离散元模拟研究。Yang等[12]基于三维离散单元法，模拟了不排水三轴压缩条件下微生物诱导碳酸钙沉淀加固砂土的宏观力学响应，揭示了其微观机理。Feng等[13]采用离散单元法模拟研究了MICP胶结砂土的特性，探讨了数值试样胶结键的断裂模式，并分析了剪胀与应力-应变行为之间的关系。Khoubani等[14]提出了一种新的胶结模型，即两个颗粒由几个较小的碳酸钙颗粒组成的环连在一起，能够模拟两个颗粒之间胶结接触的渐进破坏。以往MICP胶结模拟采用的方式是在砂颗粒之间直接施加胶结模型[8,12-13]或在颗粒接触处生成一圈规则的碳酸钙颗粒[14]，而实际情况是在砂颗粒表面及接触处均生成碳酸钙颗粒。

笔者采用离散元软件PFC3D对室内无侧限压缩试验进行模拟分析，通过在珊瑚砂颗粒表面及接触处生成微小碳酸钙颗粒来模拟MICP胶结，考虑了珊瑚砂-珊瑚砂颗粒、珊瑚砂-碳酸钙颗粒及碳酸钙-碳酸钙颗粒的不同接触模型。对建立的数值试样开展无侧限压缩试验，基于室内试验结果标定模型细观参数。在此基础上分析不同胶结程度微生物珊瑚砂固化体的颗粒位移、微裂纹发展及微裂纹分布等细观特征。

1 室内试验

1.1 珊瑚砂

试验所采用的珊瑚砂如图1所示，图2左侧曲线为珊瑚砂的粒径分布，其颗粒粒径分布较均匀，粒径均小于0.6 mm，且主要集中分布于0.18～0.425 mm之间。图3给出了珊瑚砂放大50倍的微观结构图片，从图中可以看出，试验所用珊瑚砂粒径较均匀，表面较为光滑平整。

图1 珊瑚砂[15]

图2 颗粒级配曲线Fig.2 Particle gradation

图3 珊瑚砂微观结构照片Fig.3 Microstructure photograph of coral

1.2 细菌和胶结液

试验所用的细菌为在中国普通微生物菌种保藏管理中心购买的巴斯德氏芽孢杆菌(Bacilluspasteurii)。将冷冻菌粉在NaCl溶液中复活2 h，随后加入到配制好的液体培养基中。将培养基放入恒温摇床中，在150 r/min转速下培养至细菌活性达到2 mmol/(min·L)时取出并进行二次扩培，然后将扩培后得到的菌液放入冰柜冷藏备用[16]。

胶结液用于提供脲酶水解反应及沉积过程中所用到的化学组分，包括钙源和尿素。采用的胶结液由2 mol/L氯化钙溶液和2 mol/L尿素溶液等体积混合而成[15]。细菌与尿素溶液体积比选1∶1是为了确保尿素能完全水解。菌液的浓度、活性、胶结液浓度的确定依据来自已有的MICP固化砂土的研究[17]。

1.3 MICP固化过程

试验采用的固化模具为5 cm的透明有机玻璃管(对半模具)，首先在玻璃管底部放入一层生化膜，然后放入珊瑚砂样，砂样放置完毕后，在其顶端放置一层生化膜[15]。有机玻璃管砂样体积约为180 cm3，砂样干密度为1.67 g/cm3。

试样装填完毕后，将培养好的菌液(0.5倍珊瑚砂体积)从有机玻璃管顶部注入珊瑚砂中，底部流出；为保证后续生化反应充分进行，重复2～3次。取相当于珊瑚砂体积的混合液，由有机玻璃管顶部注入珊瑚砂中，从底部流出，重复2～3遍，注入过程中注意控制流速稳定。将菌液和培养液按上述步骤每天注入1次，进行3～6 d，得到不同胶结程度的珊瑚砂柱[15]。

1.4 室内试验结果分析

对固化后的珊瑚砂柱表面进行磨平，随后置于50 ℃的烘箱48 h烘干处理，得到1#～4#珊瑚砂柱，其干密度依次为1.762、1.806、1.851、1.903 g/cm3，其中1#试样胶结程度较低，4#试样胶结程度较高。1#～4#试样碳酸钙生成量占未胶结试样质量的百分比分别为5.5%、8.1%、10.8%、14.0%。随后开展无侧限抗压强度试验，1#～4#固化试样均为脆性破坏，不同胶结程度珊瑚砂柱的应力-应变曲线如图4所示，试验结果如表1所示。从图4可以看出，无侧限压缩应力-应变曲线可大致分为3个阶段。第1阶段，应力随应变缓慢增长，呈压密状态；第2阶段，应力随应变快速增长；第3阶段，应力发生突降，试样产生压裂脆性破坏。

表1 珊瑚砂柱试样抗压强度Table 1 Compressive strengths of coral sand column

图4 不同胶结程度珊瑚砂柱的应力-应变曲线Fig.4 Stress-strain curve of coral sand column with

2 离散元模拟

采用离散元软件PFC3D(Particle Flow Code3D)对室内试验进行模拟分析。为了使离散元模拟结果与室内试验结果具有可比性，按照试验样品中所涉及的相似构型制备了数值试样，具体细节如下：

2.1 数值试样制备

为了模拟室内无侧限压缩试验，制备了如图5所示的试样，数值试样尺寸与室内试样大小一致。试样由两种类型的颗粒组成：珊瑚砂颗粒和碳酸钙颗粒。

图5 生成的数值试样Fig.5 Generated numerical

数值试样中的砂颗粒由球形颗粒表示，并按照室内试验珊瑚砂级配放大生成。以往的离散元数值试验结果表明，对颗粒的粒径进行相应放大对材料的物理力学行为影响较小[10,18]。为了使计算机能够负担得起运算，珊瑚砂颗粒的半径放大了12倍(具体见图2)。生成的砂颗粒总数是7 476个，颗粒数量足以产生稳定的试验结果[10,19]。

离散元模拟中常通过对颗粒与颗粒之间施加黏结模型来模拟胶结状态。过去模拟MICP胶结是通过在砂颗粒之间直接施加胶结模型[8,12-13]或在颗粒接触处生成一圈规则的碳酸钙颗粒[14]，而实际试样在砂颗粒表面及接触处均生成碳酸钙颗粒。为了能够较好地模拟MICP胶结，在砂颗粒表面及接触处生成微小的碳酸钙颗粒，并施加平行黏结模型将其附着在珊瑚砂颗粒上，生成的数值试样如图5所示。数值模型考虑了珊瑚砂-珊瑚砂颗粒、珊瑚砂-碳酸钙颗粒及碳酸钙-碳酸钙颗粒之间的接触，能够很好地模拟两个砂颗粒胶结的渐进破坏。这种采用微小颗粒填充孔隙的方法在Wang等[20-21]的水泥固化砂二维离散元数值试验中被首次提出，随后，Brugada等[22]在三维中用这种方法来模拟具有填孔水合物分布的甲烷水合物土，两者都取得了良好的模拟效果。这种方法生成的数值试样不仅能够较真实地模拟胶结情况，而且能够较好地控制胶结物的含量[10,20-21]。在本文中，通过室内试验中砂颗粒胶结前后的干密度可以计算出1#～4#试样的碳酸钙含量，相对应的碳酸钙含量也被分配给4个数值试样。离散元模拟中珊瑚砂颗粒的密度取2 700 kg/m3，碳酸钙颗粒的密度取3 150 kg/m3，通过4个数值试样不同的碳酸钙含量计算出对应的颗粒数量。最终得到数值试样的碳酸钙颗粒个数分别为:1#试样16 760个、2#试样24 775个、3#试样32 972个、4#试样42 445个。通过对试样上下两面刚性墙施加一定的速度，加载速率为5.4×10-3m/s，以达到加载数值试样的目的，并在压缩过程中对试样的应力-应变及微裂纹等进行监测。

2.2 接触力学模型及细观参数标定

离散元试样接触模型如图6所示，珊瑚砂-珊瑚砂颗粒采用线性接触模型、珊瑚砂-碳酸钙颗粒及碳酸钙-碳酸钙颗粒采用平行黏结模型来模拟真实试样中的胶结状态。平行黏结模型用以描述有限距离内颗粒间有夹层材料或胶结材料的本构特性，在颗粒的界面处可以是矩形或者圆形，颗粒与颗粒之间的平行黏结可以看作是一组均匀分布且具有恒定法向刚度和切向刚度的弹簧分布于接触面内[23]。平行黏结模型利用胶结键(bond)把接触间隙(gap)内的颗粒胶结在一起，胶结的两个颗粒发生相对运动时，胶结键会产生作用力与力矩，同时，力与弯矩再反作用于颗粒。当两个颗粒之间的法向力或切向力任意一者超过黏结强度值时，胶结发生破坏并产生相应的微裂纹。其中，由于法向应力过大导致胶结破坏产生拉裂纹，切向应力过大导致胶结破坏产生剪裂纹。

图6 离散元模型接触Fig.6 Contact of discrete element

与有限元能够直接利用材料的宏观参数不同，离散元模拟需要输入颗粒的细观参数。珊瑚砂颗粒之间采用线性接触模型，颗粒之间的有效模量E*通过式(1)～式(3)(PFC 5.0 手册)计算。

(1)

其中：

L=R(1)+R(2),ball-ball

(2)

A=πr2

(3)

式中：kn为法向刚度，通过珊瑚砂的三轴数值试验参数标定[24]，确定为2×106N/m；R(1)、R(2)为ball-ball接触中两个球的半径；r为ball-ball接触中两个球的半径较小值。

确定珊瑚砂颗粒之间的线性接触参数后，对数值试样开展无侧限压缩试验，试样的上下两端有两块方形加载板(沿z轴分配数值相等和方向相反的速度)对数值试样进行加载。根据室内试验2#试样的应力-应变曲线标定数值试样的平行黏结细观参数，具体数值模拟参数如表2、表3所示。

表2 材料细观参数表Table 2 Table of material meso-parameters

表3 接触细观参数表Table 3 Table of contact meso-parameters

3 数值模拟结果分析与讨论

离散元模拟选取2#试样的无侧限压缩试验应力-应变曲线进行标定，图7为2#离散元模拟和室内试验应力-应变曲线的对比。需要特别说明的是，PFC软件模拟中应力-应变曲线不能体现试样压缩的初始压密阶段[25-26]。从二者的曲线可以看出，数值试样的无侧限抗压强度与室内试验极为接近，数值试样的应变及峰后应变的走势与试验试样曲线对应得很好。

图7 2#试样数值试验和室内试验无侧限压缩曲线对比图Fig.7 Comparison of unconfined compression curves between numerical test and laboratory test for sample

为了进一步验证参数选取的正确性，在同一套参数下对1#、3#、4#试样开展数值试验。各试样的无侧限抗压强度与试验对比结果如图8所示，由于模拟中珊瑚砂颗粒分布的随机性，试验结果与模拟结果存在一些误差，但总体上，数值试样抗压强度与真实试样接近，说明在1#～4#试样胶结程度范围内参数选取合理，故可以进一步对数值试样无侧限压缩的模拟结果进行分析。

图8 无侧限压缩试验的模拟结果与试验结果对比图Fig.8 Comparison of simulation results and test results of uniaxial compression

3.1 数值试样破坏分析

在数值试验中，对胶结破坏进行全过程监测，并基于离散裂缝网络(Discrete Fracture Network，DFN)获取胶结键破坏的位置、大小及方向，根据上述信息创建微裂纹。通过PFC软件对数值试样无侧限压缩过程微裂纹及破坏碎片进行跟踪，进而对不同胶结程度的数值试样展开压缩破坏分析。离散元数值试样宏观破坏如图9所示。从图中可以看出，4个数值试样的破坏主要发生在端部，并伴随有其他位置的碎片脱落。其中，胶结程度低的试样发生破坏时以大块胶结体一起脱落破坏为主，而胶结程度高的试样压缩破坏时表现为小块或零散颗粒脱落破坏。这可能是因为胶结程度低的试样胶结并不均匀，一些地方胶结很弱，受到外力作用时首先发生破坏并相互贯通，最终呈现出大块破坏的宏观表现，而胶结程度高的试样胶结比较均匀，不会产生大块脱落破坏的情况，以小块或零散颗粒脱落为主。随着试样胶结程度的提高，试样从大块脱落破坏向小块或零散颗粒脱落破坏转变。

图9 数值试样破坏图Fig.9 Damage diagram of numerical

图10为1#试样和4#试样的室内无侧限压缩破坏形态图。从图中可以看出，胶结程度较低的1#试样以大块脱落破坏为主，而胶结程度较高的4#试样未出现大块脱落破坏，表现为小块或零散颗粒脱落破坏。数值模拟结果与室内试验结果较吻合。

图10 室内无侧限压缩破坏形态图Fig.10 Failure pattern diagram of unconfined compression

3.2 数值试样颗粒位移分析

图11为不同胶结程度试样压缩破坏后的颗粒位移情况。图中箭头方向表示颗粒移动的方向，不同颜色的箭头分别对应着不同的颗粒位移大小。从图中可以看出，胶结程度较低的试样在压缩破坏后，其右上角有一近似三角区域颗粒位移与其他部分颗粒位移有所脱离，该区域颗粒位移整体向右下方，且位移较大，在1 mm以上。胶结程度高的试样未发生较大区域与其他部分脱离的现象，可以看到零散箭头较多，以小块或零散颗粒脱离破坏为主，且底部黄色区域明显大于胶结程度低的试样。

图11 试样颗粒位移图Fig.11 Particle displacement diagram of

从整体上来看，4个数值试样在无侧限压缩后的位移云图可大致分为3部分：上下两端及中间部分。试样上下两端位移较大，颗粒位移在1 mm以上；中间红色区域位移较小，颗粒位移在0～1 mm之间，说明在竖向压缩的情况下，主要是试样两端压缩，而中间部分相对压缩较少。

3.3 数值试样微裂纹发展分析

胶结试样的宏观破坏是由微裂纹的萌生和扩展引起的，这些信息可以通过监测微裂纹的生成和记录微观拉裂纹和剪裂纹的数目变化来获取。不同胶结程度数值试样微裂纹发展如图12所示。从图中可以看出，数值试样的离散元试样破坏过程可分为以下3个不同的阶段：

图12 数值试样微裂纹发展图Fig.12 Microcrack development diagram of

第1阶段：压密阶段。这一阶段无微观裂纹产生，试样压缩变形主要来自微生物胶结试样缺陷的闭合和内部孔隙的压密。

第2阶段：裂纹扩展阶段。首先产生的是拉裂纹，经过一段时间后，试样开始萌生剪裂纹，随后，拉裂纹和剪裂纹快速扩展，随着应力的增大，颗粒间的胶结接触发生大量破坏，拉裂纹和剪裂纹数目分别以较为恒定的速度增加，且拉裂纹的增长速度大于剪裂纹，应力增长到最大值标志着这一阶段的结束。

第3阶段：裂纹急剧增长阶段。峰值应力的出现标志着这一阶段的开始，峰值过后，侧向应变和体积应变继续增大，应力发生突降，颗粒之间的胶结接触破坏更加剧烈，数值试样的拉裂纹和剪裂纹数目也急剧增加，最终导致试样宏观破坏。

3.4 数值试样微裂纹分布分析

1#～4#数值试样微裂纹最终的生成情况如图13所示，蓝色表示由于剪切导致的黏结破坏，其剪切强度已经超过黏结强度值；红色表示由于拉伸导致的黏结破坏，其拉伸强度已经超过颗粒之间的拉伸强度值。

图13 数值试样微裂纹生成图Fig.13 Microcrack formation diagram of

4个试样黏结接触均发生了拉伸和剪切破坏，产生了微观拉裂纹和剪裂纹。其中，1#数值试样胶结程度最低，在加载破坏后，拉裂纹数目为4 964个，剪裂纹数目为2 546个，拉裂纹和剪裂纹数目之比为1.95；在加载破坏后,胶结程度较低的2#试样拉裂纹和剪裂纹数目分别为6 094、3 411个，其拉裂纹和剪裂纹数目之比为1.78；在加载破坏后，胶结程度较低的3#数值试样拉裂纹数目为9 462个，剪裂纹数目为5 605个，拉裂纹和剪裂纹数目之比为1.69；无侧限压缩破坏后，微生物胶结程度较高的4#试样拉裂纹和剪裂纹数目分别为19 964、15 969个，拉裂纹和剪裂纹数目之比为1.25。在外力压缩作用下，胶结程度低的试样主要是胶结弱的地方首先发生破坏并相互贯通，而胶结程度高的试样压缩破坏时表现为小块或零散的颗粒脱落，并且随着胶结程度的提高，胶结键数量也随之增加，故微裂纹的数目会随着胶结程度的提高而增多。

从4个不同胶结程度的试样在无侧限压缩破坏后的拉裂纹和剪裂纹数目可以看出，生成的拉裂纹数量多于剪裂纹数量，以拉破坏为主。在压缩破坏的情况下，随着微生物胶结程度的提高，拉裂纹和剪裂纹数目比值逐渐减小。说明拉、剪裂纹数目比值越小，试样的无侧限抗压强度越高。

为了表征及量化不同胶结程度的微生物固化试样在无侧限压缩后内部微裂纹在空间上的分布，通过fish语言遍历所有微裂纹，得到其在xoy平面上的位置分布信息，统计每一个方向上的裂纹数目，并绘成如图14所示的玫瑰图。玫瑰图的极角θ表示微裂纹的位置向量与x轴的夹角，玫瑰图上的点距离圆心越远，表示在该方向上微裂纹数目越多。从图中可以看出，胶结程度最低的1#试样各个方向微裂纹数目差别较大，在100°、155°、210°及310°方向附近微裂纹数目较多，其余方向微裂纹分布较少，试样在受到压缩时，微裂纹在各个方向扩展并不均匀，裂纹数目较多的方向首先发展成宏观裂缝，破坏时的强度较低；而胶结程度最高的4#试样微裂纹各个方向分布比较均匀，在受到压缩时，各个方向裂纹扩展比较均匀，当裂纹扩展达到一定程度时，试样发生破坏，所需的压缩强度较高。2#试样和3#试样的微裂纹分布玫瑰图介于1#试样和4#试样的过渡状态。