制样方法对砂土液化力学性质影响的离散元模拟

摘要

基于PFC2D软件对湿捣法和落砂法制备的砂土试样进行不排水循环剪切试验模拟，分析了砂土颗粒在动力加载过程中细观结构的变化规律。结果表明，施加荷载后，试样颗粒长轴方向在压缩侧偏向水平，拉伸侧趋向竖直，而接触法向的趋势则与之相反，二者的各向异性系数均逐渐上升。试样整体配位数随荷载的作用逐渐下降，液化后于0~3范围内波动。湿捣法试样比落砂法试样的各向异性程度更小，液化所需的加载循环数更多，即具有更高的抗液化强度。密实度越高，两种试样的液化强度差别越大。

关键词

砂土液化; 制样方法; 数值模拟; 结构特征

砂土液化是一种破坏性很强的地质灾害。砂土液化有可能引发大面积地基失稳、建筑物倒塌、桩基滑移、地下管线上浮等工程问题，给人民的生命财产安全构成巨大威胁。为了深入理解液化灾害的成因和机理，学者们针对砂土液化问题开展了大量的室内试验研究。由于砂土不具有黏聚力，常规的取样设备较难获得原状土样，因此现有的室内试验大多采用重塑法进行试样制备。学者^［1-2］通过试验发现，在相同的应力状态和加载条件下，不同制样方法制备的砂样呈现出不同的液化行为和强度特征。这主要是由于不同制样方法产生的不同初始结构会影响砂土的力学行为特性^［3-6］。目前，已有多种方法可对砂土细观结构进行直接观测。然而，这些方法大多只能获取某一时刻或某一平面的试样细观结构特征，无法观测整个加载过程中细观结构的变化。此外，在试验中，大部分细观结构的获取过程都不可避免会对试样的细观结构产生扰动。近年来，学者采用数值模拟方法研究土体内部结构特征，可避免以上研究中的不足。离散单元法作为一种能够反映颗粒材料力学特性的方法^［7］，在模拟和分析砂土材料细观结构变化过程方面具有独特的优势。PFC是最常用的离散元法软件之一，不仅能模拟砂土的宏观力学响应，也能用于研究砂土试样内部的细观结构演化特征。已有研究^［8-11］证明离散单元法能有效模拟砂土在土工试验中的宏观力学响应，并能够分析试样的细观结构特征。然而，现有的研究大多关注土体的强度特征，针对制样方法对砂土液化特性影响的分析报道较少。

因此，本文采用PFC2D软件首先模拟了两种常用制样方法（湿捣法和落砂法）的试样制备过程。接着，对制备好的试样进行不排水循环加载模拟，分析了不同密度和循环强度条件下，不同制样方法制备的试样的宏观液化性质和细观结构演化规律。最后，利用细观结构变化特征解释试样的宏观力学响应，从而揭示制样方法对砂土液化特性的影响机理。

1 制样过程模拟

实际场地中的砂土颗粒具有不规则棱角，这将增加颗粒间的接触数量和面积，进而影响试样的结构分布和强度特征。因此，为了更符合实际情况，采用PFC2D软件内置的Clump功能，将三个圆形颗粒结合成不可破坏的近椭圆形颗粒，形状如图1所示，其中左右两个次圆的切点为中间主圆圆心。研究表明^［9-10］，峰值内摩擦角和试样发生破坏时的破坏应变随颗粒长短轴比的变化而变化。当颗粒长短轴大于1.5后，颗粒长短轴的变化对峰值强度的贡献不明显。因此，学者通常采用长短轴为1.4模拟具有亚菱角状的颗粒材料，因此本研究设置颗粒长短轴比l/h为1.4。基于已有研究^［11-13］，砂土颗粒法向刚度设为2×10⁸ N·m^-1，切向刚度设为1×10⁸ N·m^-1，颗粒间摩擦系数设为0.5。为了防止颗粒在制样及加载过程中溢出墙体，将墙体刚度进行放大，并取颗粒与墙体摩擦系数为0以减小边界效应的影响。颗粒生成时依据面积等效原则，保证不规则clump颗粒与直径0.6~1.8mm的圆形颗粒面积相当，具体颗粒级配曲线如图1所示。本次模拟颗粒间应力-应变关系采用线弹性模型，其他参数如表1所示。

图1  颗粒示例及粒径级配曲线

Fig.1  Particle shape and particle size distribution curve

表1  颗粒参数

Tab.1  Parameters of the simulation tests

属性名称	数值
颗粒密度/（kg∙m3）	2 650
颗粒法向刚度/（N∙m-1） 颗粒切向刚度/（N∙m-1）	2×108 1×108
墙体法向刚度/（N∙m-1） 墙体切向刚度/（N∙m-1）	2×109 1×109
循环加载时颗粒间摩擦系数	0.5
颗粒-墙体间摩擦系数	0
循环加载时阻尼系数	0.7

本文模拟了土工试验中两种最常见的砂土制样方法：落砂法和湿捣法。PFC软件中有墙体和颗粒两种实体，模拟制样时由颗粒模拟土体，墙体模拟容器，通过调整颗粒生成方式和墙体加载过程来模拟不同制样方法。常用的模拟制样方法有半径膨胀法、重力沉积法、分层压实法等^［12-13］。基于土工试验中不同制样方法的特点，在数值模拟中分别采用重力沉积法模拟落砂法，分层压实法模拟湿捣法。

落砂法制样步骤如下：

（1）生成左右及底部墙体，并于底部区域内生成前述近椭圆形颗粒；

⑵施加重力让颗粒自由下落；

⑶于较高区域生成新颗粒；

⑷新颗粒下落与旧颗粒均匀接触；

⑸重复以上步骤至试样达到预设高度，并于此高度生成顶部墙体，将所有位于墙体范围外颗粒删除，至此落砂法制样完毕。

湿捣法步骤如下：

⑴生成四周墙体，并于底部生成近椭圆形颗粒；

⑵其他三面墙体不动，顶部墙体下压至设定高度，压实颗粒；

⑶顶部墙体重新上升，于较高处空间生成新颗粒；

⑷墙体再次下降并压实颗粒；

⑸重复以上步骤直到试样到达指定高度，至此湿捣法制样完毕。

两种方法的实现过程示意图如图2所示。其中，落砂法依靠颗粒自身重力完成压实过程，而湿捣法则通过墙体施压来压实土样。试样尺寸为40×80 mm，共约2 100个颗粒。制样完成后，通过施加200 kPa的墙体围压进行试样固结。本文共设置4个模拟试样，分别为落砂法制备的中密砂（e=0.167）和密砂（e=0.155），湿捣法制备的中密砂（e=0.167）和密砂（e=0.153）。

图2  落砂法和湿捣法在PFC软件中制样过程

Fig.2  Simulation process of deposition method and layer method in PFC

2 循环加载试验结果

本文对落砂法和湿捣法生成的4个试样，分别进行了不排水循环加载试验模拟。其中，循环加载速率通过设置墙体的移动速度来实现，不排水条件通过保持试样面积在循环加载过程中不变来实现^［14］。模拟试验考虑了密砂和中密砂两种初始状态，以及动应力比（CSR）为0.20和0.15两种加载情况。循环加载过程采用等应力幅值控制。加载时，上下墙体向内收缩，左右墙体向外扩张，当达到预设动应力比时，所有墙体速度反向，以实现卸载过程。在同等密实度条件下，湿捣法和落砂法试样的应力状态和加载条件完全一致，其中中密砂的CSR为0.15，密砂的CSR为0.20。中密砂和密砂的应力路径和孔压时程曲线如图3和4所示。

图3  中密砂试样应力路径和孔压时程曲线

Fig.3  Stress paths and pore water pressure history of medium dense sand

图4  密砂试样应力路径和孔压时程曲线

Fig.4  Stress paths and pore water pressure history of dense sand

从图3可以看出，湿捣法试样到达液化状态（双幅轴向应变为5%）所需的循环数N_c比落砂法更大，这意味着在相同应力状态条件下，湿捣法制备的试样比落砂法试样强度更高，更能抵抗外力荷载的作用。根据Ni等^［15］研究，砂土在第一循环内的超孔压累积量影响其到达相变线的速度，从而影响抗液化强度。从湿捣法和落砂法试样在第一循环内有效应力的下降量（即超孔压的增长量）可见，落砂法试样更快到达相变线，因此更快发生液化。图4所示为密实状态下两种砂的应力路径，可以发现类似的强度规律。值得注意的是，湿捣法制备的试样到达液化的N_c比落砂法大，且两者的差距在密砂中体现得更明显，这与Ni等^［16］进行的室内动三轴试验（如图5所示）得出的结论一致。综上，在中密实及密实状态下，湿捣法制备的砂样结构性更强，能抵抗更大的外力荷载作用。

图5  中密砂和密砂的应力路径（动三轴试验结果^［16］）

Fig. 5  Stress paths of medium dense and dense sands (cyclic triaxial test results^[16])

3 细观结构分析

为了深入理解上节中观测到的宏观力学行为，本节通过一系列细观结构特征分析不同制样方法生成的试样细观结构差异性，主要包括颗粒长轴方向、颗粒接触法向方向和颗粒配位数。

利用DEM能够采集到试样的结构特征，通过统计和分析这些结构特征，进一步可得试样的孔隙率、配位数、颗粒长轴方向、颗粒接触方向等各种结构参数。Rothenburg^［17］在研究中指出，频率分布图能够描述平面细观结构的分布规律。对于颗粒长轴方向及接触法向分布，可由公式（1）进行傅里叶变换来近似。

$$E\left(\theta \right)=\frac{\mathrm{1}}{\mathrm{2}\pi }\left[\mathrm{1}+a\cdot \mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\left(\theta -{\theta }_a\right)\right]$$

（1）

式中： $E(\theta )$为分布的频率函数；${\theta }_a$为粒分布的主方向；a为颗粒分布的各向异性系数，a值越大，各向异性程度越大，满足：

$${\int }_{\mathrm{0}}^{\mathrm{2}\pi }E\left(\theta \right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta \mathrm{d}\theta =\frac{a}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}{\theta }_a$$

（2）

$${\int }_{\mathrm{0}}^{\mathrm{2}\pi }E\left(\theta \right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}\theta \mathrm{d}\theta =\frac{a}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}{\theta }_a$$

（3）

3.1　颗粒长轴方向分布特征

为了观察循环加载过程中细观结构的变化及演化规律，在图4的应力路径图中，取偏应力为压缩侧幅值对应为A系列点，拉伸侧幅值对应B系列点和加载对称轴对应C系列点，共计15组应力状态。A系列点分别标记为A1（主应力200 kPa）、A2（主应力150 kPa）、A3（主应力100 kPa）、A4（初始液化状态）和A5（液化后第3个循环）。B、C系列点与A系列对应的点处于同一加载循环内。

图6为湿捣法制备的中密砂试样在循环加载过程中颗粒长轴方向分布图。由颗粒长轴方向获取的各向异性系数a反映了试样整体颗粒分布的各向异性程度。随着循环荷载的进行，试样强度逐渐下降，但液化前（点1~3）颗粒长轴方向分布没有明显的变化，在此期间试样变形不明显。当试样到达初始液化状态（点4~5），颗粒长轴方向分布发生明显的变化，当进行到压缩侧时（A点），颗粒长轴偏向水平分布，a值较液化前逐渐增大；当进行到拉伸侧时（B点），颗粒长轴趋于竖直分布，a值较之前也逐渐增大；当加载到对称轴时（C点），颗粒长轴分布也偏向水平方向，与压缩侧（A点）趋势更接近。基于Ishihara和Okada^［18］的研究，应力状态在拉、压两侧的相变线之间变化时，土体的变形较小，结构变化不明显；当应力状态越过相变线后，试样结构分布发生明显变化，较大的变形将被观察到。图6所示结果可直观上描述应力状态在相变状态前后的颗粒分布特征。

图6  湿捣法试样循环加载过程中颗粒长轴方向分布（中密砂）

Fig.6  Distribution of long axis direction of layered sand particles during cyclic loading (medium dense sand)

图7所示为不同试样在A、B和C状态点颗粒长轴方向的各向异性系数a变化图。其中0位置代表砂土加载前的初始各向异性系数，1~5表示的位置如图4所示。从图7中可以看到，对于同种制样方法，密度越高，试样的初始各向异性系数越大，即各向异性程度越高。在相同密度条件下，湿捣法试样的a值比落砂法的小，说明湿捣法制备的试样颗粒长轴分布更趋于各向同性，这与室内试验中湿捣法的结构特征类似。整体而言，砂土到达液化状态之前（点1~3），各向异性系数相差不大。结合颗粒长轴方向分布图特征可知，此时颗粒尚未产生明显的移动和转动。当试样到达液化状态后（点4~5），湿捣法试样的各向异性系数略有所上升，压缩侧、拉伸侧以及加载对称轴上试样的a值差别不大，但由图6可知其分布方向不同。落砂法试样在压缩侧（点A）和加载对称轴（点C）的a值几乎一致，但随液化程度的加大呈下降趋势。不同的是，试样在拉伸侧（点B）的a值出现了明显的下降，此时颗粒长轴发生了明显变化，部分颗粒发生转向不再偏水平分布。此外，对于中密砂和密砂，拉伸侧的a值下降的幅度相近。由此可见，液化后砂土试样的结构发生明显的破坏和重组。对于未经受液化历史的砂土，砂土液化强度取决于土体的初始结构特征，各向异性程度越高，抵抗外力荷载的能力越差。

图7  试样在不同位置点处的颗粒长轴方向各向异性系数变化

Fig.7  Relationship between anisotropy coefficient of long axis direction and the loading locations for sand samples

3.2　颗粒接触法向方向分布特征

在PFC模拟中，颗粒接触力定义为颗粒与墙体或颗粒之间根据接触来传递相互作用力。基于式（2）~（3）和接触力法向（即接触力的方向）分布图获取的各向异性系数a，能够反映试样内部颗粒间接触力方向分布的均匀程度。图8为湿捣法制备的中密砂试样在循环加载过程中接触法向的分布图。可以看到，随着循环加载的进行，在 A、B、C应力状态点，砂颗粒间接触法向分布的各向异性系数a逐渐增加。当试样进入液化状态后（点4~5），由于外力荷载作用方向的不同，A、B点的接触法向方向分布完全不同，试样在压缩侧（点A）的接触法向集中在竖直方向分布，在拉伸侧（点B）则集中在水平方向分布。试样在加载对称轴（点C）的接触法向分布与压缩侧类似。液化后由于各向异性系数已处于较大值，呈缓慢增长趋势。

图8  湿捣法试样循环加载过程中接触法向分布（中密砂）

Fig.8  Distribution of contact direction of layered sand particles during cyclic loading (medium dense sand)

图9为不同试样在A、B和C状态点接触法向的各向异性系数a变化图。可以发现，对于相同制样方法制备的试样，密实度越高，各向异性系数a越大。而对于相同密度的砂土，落砂法试样的a值大于湿捣法，即湿捣法试样的接触法向在各个方向分布更均匀。到达液化之前（点1~3），湿捣法试样A点的a值较初始值大约增长了0.12，B点的a值亦有所增加但增量小于A点，C点变化不明显，各位置点各向异性系数整体呈上升趋势。对于落砂法试样，A状态的a值均增大了0.08，增长幅度虽小于湿捣法，但实际数值仍大于湿捣法试样，然而B状态点明显下降了0.15，这主要是受到轴向拉伸荷载的作用，C状态点稍有下降但变化不大。密实度大小则对液化前各位置点的a值变化幅度及演化规律影响不大。

图9  不同试样各位置点接触法向各向异性系数变化

Fig.9  Relationship between anisotropy coefficient of contact direction and the loading locations for sand samples

当试样发生液化后（点4~5），所有试样在A、B、C状态处的a值均显著增大，说明土体发生液化破坏了初始的颗粒接触结构，重新形成的接触法向在轴向压缩侧（点A）集中于竖直方向分布，在拉伸侧（点B）则集中在水平方向，呈明显的各向异性分布。此时，湿捣法和落砂法制备的中密砂试样在A、B、C状态点的各向异性系数a差别不大，但落砂密砂试样的a值略大于湿捣密砂试样。

由图6及图8可以发现，施加循环荷载后，砂颗粒接触法向的a值在液化前（点1~3）缓慢增长，并且在压缩侧（点A）偏向竖直，在拉伸侧（点B）偏向水平，但是在拉、压两侧的颗粒长轴方向没有明显差异，说明接触法向的变化先于颗粒长轴方向。液化后（点4~5）接触法向a值显著增大，各向异性程度明显，颗粒长轴方向a值亦开始增长，但幅值仍小于前者。其中，砂颗粒在点A状态的接触力呈竖直方向分布，而颗粒长轴方向趋向水平分布，B点的方向分布正好相反。换言之，试样受到循环荷载后，内部接触力各向异性程度增加，导致试样强度下降。当试样进入初始液化状态后，颗粒长轴受接触力影响在重排列时也表现出更大的各向异性。无论是接触法向还是颗粒长轴方向，各向同性结构具有较高的强度，但接触法向的影响较为明显。

综合图7与图9可知，相同制样方法制备的砂样，密实度越高，初始各向异性系数越大；在密实度相同条件下，落砂法试样的初始各向异性系数大于湿捣法。在循环荷载作用下，落砂法试样在点A、C状态的a值增长幅度略小于湿捣法，但因其初始各向异性系数较大，a值仍高于湿捣法。结合图3和图4的应力路径图发现，在相同的应力状态及加载条件下，颗粒接触法向分布各向异性系数a值较小的试样到达液化所需的循环数N_c比系数较大的试样更多。换言之，各向异性程度越小，砂土抗液化能力越强。

3.3　颗粒配位数特征

配位数为试样接触总数与颗粒总数的比值，是反应试样结构特性的定量参数之一。配位数越大，每个颗粒与其他颗粒的平均接触就越多，颗粒间越难产生相对位移，试样配位数的下降往往伴随试样强度的降低。

图10给出不同试样配位数随运算步数变化图，其中虚线标记了试样进入液化状态时的运算步数。图10 a为湿捣法中密砂试样的变化图。如图所示，在初始加载阶段，试样配位数与初始配位数相近，维持在3.7附近，颗粒间接触紧密，试样具有较高强度，能够抵抗外力荷载的作用。随着循环次数的增加，配位数减小至3左右，颗粒间仍接触较充分。当试样进入液化状态时，配位数急剧下降至1以下，此时颗粒已与周围颗粒无直接接触，呈“悬浮”状态，试样强度急剧丧失。接着轴向应变增大，颗粒间重新形成接触直至足够承担外部施加的循环荷载，配位数回升至3附近，此后配位数往复波动于0~3之间。同样对于中密砂，落砂法制备的试样配位数如图10c所示，试样配位数随循环荷载作用呈持续下降趋势，配位数无维持阶段。进入循环液化后，配位数在0~3之间循环波动。对于密砂，湿捣法和落砂法试样的配位数在加载过程中的变化规律与密砂类似。湿捣法密砂试样的配位数下降速度远小于落砂法密砂试样。

图10  不同试样配位数随运算步数变化

Fig.10  Process of coordination number with different samples

比较图中虚线位置可以发现，随着密实度的增大，湿捣法试样的强度急剧上升，因此中密砂与密砂之间的强度差别越来越大。而对于落砂法，中密砂与密砂之间的抗液化强度差别相对较小，这与Ni等^［17］的试验结果一致。根据Ni等^［17］的研究，这种差别主要是由于土体各向异性程度、局部孔隙分布以及孔隙大小的变化造成的。

湿捣法试样与落砂法试样的对比结果表明，即使初始配位数几乎相同，湿捣法试样比落砂法试样配位数下降更缓慢，具有更高的抗液化强度，到达液化需要更多次数的循环荷载。液化以后各试样配位数均随循环荷载反复波动，波动范围及趋势相近。

4 结论

本文利用离散元方法模拟了落砂法和湿捣法两种制样方式，对这两种不同方式制备的试样进行了不排水循环加载试验，分析了土体的应力路径，应力-应变关系及其液化强度特征。并且采集了从初始试样到加载至液化后全过程的细观结构数据，分析了颗粒长轴方向、接触法向以及试样配位数随循环荷载作用的变化规律，得到以下结论：

（1）密砂及中密砂试样在相同的应力条件下，湿捣法所得试样达到液化状态所需的加载循环数N_c大于落砂法，且密砂的循环数差距大于中密砂。无论是落砂法还是湿捣法，液化后密砂的轴向应变幅值均小于中密砂。

（2）湿捣法试样的颗粒长轴方向分布更趋向各向同性，而落砂法颗粒明显偏向水平，各向异性系数较大。湿捣法试样的接触法向分布无明显偏向，各向异性系数较小，落砂法接触法向偏竖直方向，各向异性系数更大。两种制样方法随着密度的增加，颗粒长轴方向和接触法向的各向异性系数都有所增加。

（3）颗粒接触法向的各向异性系数在加载后即开始上升，增长到一定程度后试样液化，轴向应变增大，此时颗粒长轴方向才发生变化。接触法向的变化先于颗粒长轴方向的变化，且影响了其变化趋势。

（4）对于初始配位数相近的试样，液化前湿捣法所得试样受循环荷载时配位数下降较落砂法更缓慢，需要更多的循环数才能液化，具有更强的抗液化性能。液化后各试样配位数波动范围和趋势较为一致。

作者贡献声明

叶斌：研究思路指导及论文修订。

宋思聪：数值模拟及论文撰写。

倪雪倩：研究内容制定及算例验证。