砂土作为土木工程、交通运输和水利工程等领域中常用的建筑材料，其力学特性是工程建筑基础研究的重点，也是土力学研究的热点。砂土液化是指在外荷载作用下，饱和砂土由于孔隙水压力升高而引起剪切强度丧失和有效应力降低，最终导致砂土表现出类似于流体特征的一种现象。近年来，由降雨入渗和地震引发的静动力液化而导致的滑坡、泥石流等地质灾害和长期交通荷载作用下的路基累积变形等问题日益受到人们的关注^[1-4]。对松砂的静动力液化行为也开展了大量室内试验工作，包括三轴剪切试验^[5-6]和循环三轴试验^[7-8]，在砂土液化的宏观力学机制方面等取得了可喜的成果。但是砂土是离散的介质，松砂的静动力液化行为非常复杂，基于室内试验获得的宏观力学响应通常是唯象的，不能从本质上揭示砂土液化的发生机制。

随着计算机技术的发展，由Cundall和Strack^[9]提出并逐步发展起来的离散元(DEM)数值方法，可以在获得颗粒材料宏观响应的同时追踪颗粒的微观组构演化。PFC(Particle Flow Code)作为不断完善发展的离散元程序，广泛应用于砂土等颗粒材料的宏细观力学行为研究。周健等^[10]通过开展PFC2D常体积循环双轴试验，对循环荷载作用下砂土的液化现象进行了模拟。Yang和Dai^[11]利用DEM方法开展了一系列双轴剪切试验的数值模拟，指出砂土的不稳定行为是土体颗粒在等体积剪切过程中重新排布的结果。Guo和Zhao^[12]研究了砂土在剪切过程中的各向异性演化规律及临界状态力学行为。史旦达等^[13]研究了初始组构对砂土液化势的细观机理，指出砂土的液化过程是一个组构各向异性的演化过程。目前，对砂土液化的细观机理研究大多是基于二维离散元数值试验，亟需在三维条件下开展进一步数值研究。另一方面，为进一步揭示液化现象的发生机制，探讨砂土静动力液化的内在关联，有必要对砂土静动力液化现象进行数值模拟。

基于此，利用三维离散元方法，分别在静动力加载条件下，对松散砂土数值试样进行三轴不排水剪切试验模拟，研究松散砂土静动力液化的宏微观力学特性。数值模拟试验表明松砂试样在静动力荷载作用下均会发生液化现象，是砂土颗粒微观量演化的结果。数值试验结果为砂土液化的细观机理研究提供了强有力的数值试验依据，丰富了震陷、泥石流等液化失稳灾变的分析方法。

离散元方法的力学原理是基于牛顿第二定律和颗粒间的接触法则，通过显示的方法来求解颗粒的位移，实现对颗粒材料的力学模拟^[14]。采用能够描述岩土工程材料的力学特性简化的Hertz非线性接触本构模型和能消除边界对数值模拟结果影响的周期边界。其材料的细观参数如表 1所示。颗粒几何特性和颗粒级配对砂土的力学特性具有重要的影响。为了克服计算数值上的困难，试样采用球形颗粒，试样级配曲线的粒径范围为0.1~0.4 mm之间，可以代表被学者广泛采用的试验砂枫丹白露砂的级配曲线。首先在一个正方体放入小球，通过设置摩擦系数来得到不同密实度的试样。试样固结采用各向等向固结。图 1是固结结束之后的松砂试样。数值试验中通过控制试样的体积不变来实现不排水条件。循环加载采用等应力幅的加载方式进行。单调荷载作用下和循环荷载作用下的数值模拟试验方案如表 2和3所示。

图 1 数值试样 Fig. 1 Numerical sample

在常体积条件下单调荷载和循环荷载作用过程中侧向有效围压的变化反映了试样内超孔隙水压力的变化情况。在三维状态分析中，平均正应力p定义为p=(σ_xx+σ_yy+σ_zz)/3、剪应力q定义为q= 和孔隙水压力定义为u=σ_zz0-σ_zz，其中σ_xx, σ_yy, σ_zz分别为x, y, z方向的正应力。

不排水单调加载条件下松砂宏观响应如图 2所示，包括75，100，120 kPa和150 kPa等4种围压。从图 2(a)可以看出，在加载初期，偏应力迅速上升直至峰值，此时试样变形极小(ε_a < 0.5%)。而后，随着轴向应变不断累积，偏应力逐渐减小至接近于0，试样完全丧失强度。值得注意的是，砂土作为一种压硬性材料，其强度将随围压的增大而增大；数值试验结果也揭示了这一性质，P150-M试样峰值强度最大，约为72 kPa。另一方面，峰值强度所对应的应变也与围压有关；一般地，围压越大，峰值强度所对应的轴向应变越大，也即应变软化发生得越迟。图 2(b)所示为试样有效应力路径，由图可知，在不排水加载过程中，120 kPa和150 kPa围压下的两组试样，由于孔隙水压力的不断累积，有效平均正应力不断减小，试样发生剪缩；最终有效应力和偏应力均趋近于0，这表明土体不具有任何抗剪切的能力，即试样发生初始静态液化^[15]。而对于75, 100 kPa围压下的试样，有效正平均正应力先增大后减少，试样发生先剪胀后剪缩现象，这表明松砂在低围压下会出现密砂的剪胀特性。室内试验表明，静态液化在松砂中较为常见^[5]，本数值试验结果较好地模拟了室内试验中松砂的液化现象和松砂在低围压下的剪胀性。

图 2 单调荷载作用下试样宏观响应 Fig. 2 Macro response of samples under monotonic loading

两组典型不排水动力试验结果(P100-C2和P100-C3)如图 3、4所示。循环荷载作用下孔隙水压力累积速度逐渐加快(图 3(a)、4(a))，最后急剧上升至初始围压；同时有效应力逐渐减小至0，试样最终液化(图 3(b)、4(b))。由图 3(c)、4(c)可知，试样初始变形极小，轴向应变没有明显的循环效应和累积效应；而当临近液化时，应变急剧发展，试样从受拉侧瞬间破坏，完全丧失强度。进一步地，图 3所示P100-C3试样在第6个加载周期发生液化，而图 4所示P100-C2试样在第77个振次才发生液化，这说明当循环荷载幅值保持不变时，初始围压越大，砂土抗液化能力越强。P100-C4和P100-C25组试验也呈现出类似规律，限于篇幅，试验结果不在文中描述。Sze等^[8]利用不同制样方法获得了不同密实度和围压下的砂土试样，开展了一系列循环荷载试验。试验结果表明试样初始围压越大，试样越不容易发生液化现象。本试验的数值结果与Sze等^[8]的室内试验结果一致，三维离散元数值模拟可以模拟砂土的动力液化特性。

图 3 循环荷载作用下试样P100-C3宏观响应 Fig. 3 Macro response of sample P100-C3 under cyclic loading

图 4 循环荷载作用下试样P100-C2宏观响应 Fig. 4 Macro response of sample P100-C2 under cyclic loading

图 5给出了试样液化时的破坏振次N_f与循环应力比CSR之间的关系，其中CSR为动剪应力(q/2)与初始围压之比。可以看到，随着循环应力比的减小，试样破坏所需振次N_f将逐渐增大，这与曹久亭^[16]等人进行的常围压室内单元体试验结果相一致。需要指出的是，本次试验中砂土抗液化能力不仅与循环应力比有关，还需考虑试样围压的影响。根据临界状态土力学理论^[17]，砂土的“松散”程度不只由试样相对密实度决定，还受围压水平的影响。一般而言，围压越大，试样相对“越松”，因此当循环应力比相同时其动强度越低，更易发生液化。

图 5 试样破坏振次和循环应力比关系 Fig. 5 Relationship between number of vibrations and cyclic stress ratio

砂土颗粒的微观力学特性通常用接触法向、接触力、颗粒方向、配位数等微观量进行描述。通常采用组构张量对接触法向进行分析，其可以反映材料的各向异性程度。组构张量G可表示为：

(1)

其偏量部分F表示为：

(2)

其中w(m^k)为权重因子。组构偏张量F通过方向n_F和模量F两部分描述，组构张量的特性可以表述为：

(3)

而组构张量的方向nF可用洛德角θ_F来描述，两者间的关系为：

(4)

其中0°≤θ_F≤60°, 三轴压缩试验中θ_F=0°，而三轴拉伸试验中θ_F=60°。

配位数是指单个颗粒所含有接触的数目。通常采用平均配位数Z来表征其接触形态，即：

(5)

式中N_C和N_B为总的接触数和颗粒数。

不排水单调加载条件下松砂微观响应如图 6所示。从图 6(a)可以看出，所有的试样液化之后，其组构的模量都达到稳定值，约为0.62，与围压无关。Guo和Zhao^[12]也得到相同的结果。由图 6(b)可知，组构张量的方向θ_F在剪切开始阶段随轴应变迅速减小，在ε_a>1%后，趋于稳定值0°，即组构张量方向与加载方向一致。从图 6(c)可以看出，试样在剪切过程中，配位数不断下降，直至液化之后稳定在3.2。试样的配位数与颗粒材料内部结构的稳定性直接相关。由摩擦球形颗粒构成的三维试样，大于等于4的配位数是保证系统处于稳定状态的必要条件^[18]。当偏应力达到峰值之时，配位数就减少到4以下，试样处于不稳定状态。

图 6 单调荷载作用下试样微观响应 Fig. 6 Micro response of samples under monotonic loading

两组典型不排水动力试验微观响应(P100-C2和P100-C3)如图 7、图 8所示。由图 7(a)、图 8(a)可知，试样的初始组构模量很小，各向异性很弱，而临近液化时，组构的模量演化加快，直至达到稳定值0.59，比静态液化时的稳定值小。由图 7(b)、图 8(b)可知，组构方向θ_F随着加载方向的变化而不断变化，出现循环效应，当临近液化时，组构的方向θ_F趋近于60°, 也表明试样在受拉侧破坏。周建等^[10]研究结果表明组构的方向变化要滞后于应力方向的变化。循环荷载作用下配位数(图 7(c)、8(c))，最后急剧下降至3.4，比静态液化时的值大。

图 7 循环荷载作用下试样P100-C3微观响应 Fig. 7 Micro response of sample P100-C2 under cyclic loading

图 8 循环荷载作用下试样P100-C2微观响应 Fig. 8 Micro response of sample P100-C2 under cyclic loading

通过利用三维离散元软件，对松散砂土数值试样进行三轴不排水剪切试验模拟，在得到砂土静动力液化特性的同时，对加载过程中试样内部组构和配位数的变化进行了分析，并初步探讨了砂土液化的细观力学机制。得到的主要结论有：

(1) 三维离散元数值模拟可以很好的模拟松砂的静动力液化现象，室内试验中松砂的“初始液化”和低围压下松砂也可能发生剪胀等现象均可以在数值模拟中实现。

(2) 砂土颗粒的组构张量和配位数的演化规律在宏观上的体现就是砂土的宏观力学特性。

(3) 在动力加载条件下，砂土是在受拉侧发生破坏，破坏时其微观参量组构张量模量比静态液化时的稳定值小，而配位数比静态液化时的值大。

本研究仅对松砂的静动力液化现象进行了三维离散元模拟和分析，可以增加不同密实度的砂土进行数值模拟分析，进一步分析砂土动力循环液化的细观机理。