地震引起的场地液化及其导致的地基失效，如：不均匀沉降、流滑、侧向扩展、强度下降等，是地震灾害中建筑、桥梁、码头等损坏的重要原因。饱和或接近饱和的松散至中密的砂土与粉土等粒状土在地震烈度7度及以上的震动条件下容易发生液化，饱和的粗砾土、堆石料及较密实砂土在震动强度较高且排水不畅时也可能发生液化[1]。碎石桩加固是场地抗液化地基处理的有效措施之一[2－3]。因施工简单且成本较低，强震地区的工程中常采用碎石桩加固处理地下水位较高的砂土或粉土等可液化场地以提高其抗液化能力。Seed和Booker[2]最先对饱和砂土中碎石桩抗液化效果进行研究并给出可供工程应用的评估图表。其后研究者对碎石桩抗液化效果在地震现场调查、理论分析、模型试验，及数值模拟等方面[4―13]进行了广泛研究。碎石桩加固有利于增加土体密实度，其良好的渗透性可减少周边土体中超静孔压累积并加快孔压消散，且碎石桩体比周围土体刚度高可减少桩周土体中的应力分配并对场地变形起到一定的约束效应，因此，碎石桩可改善可液化场地的抗液化能力，其对周围土体在地震作用下的效应是多重因素发挥影响的综合效应。
土体液化复杂力学行为本身需要精细和合理的描述，是岩土领域研究的热点之一[1]，而碎石桩为散体材料桩，其力学特性和变形又会受到周边土体围压与变形的影响，因此，可液化场地中碎石桩复合地基的动力响应研究就更为复杂。碎石桩复合地基的变形、应力、孔压变化等动力特性还有很多待研究之处，如：复合地基变形及孔压空间变化特征，周边未加固土体发生震动液化对碎石桩加固区边缘与中间部位的抗液化效果的影响，及桩土应力分配变化等。
Brennan[5]对碎石桩加固饱和砂土抗液化效果进行了系列离心动力试验，研究了排水及群桩作用的抗液化效果。Adalier等[6]通过离心模型试验研究了碎石桩刚度对饱和粉土的抗液化影响，试验用不透水薄膜将碎石桩隔离而消除了其排水效应作用。Badanagki等[7]对碎石桩单桩在水平场地和碎石桩群在微倾斜场地的抗液化效果进行了离心动力试验研究，场地为包含可液化砂土和粉土的多层土场地。这些试验分别揭示了碎石桩的刚度、排水或对土体加密等某些方面效应及相关认知，可作为数值模拟分析的标定基础。Papadimitriou等[8]和Li等[9]分别采用不同的二维砂土液化模型对Brennan[5]的试验AJB-8a和Badanagki等[7]的试验1、试验2进行了模拟，将碎石桩按刚度相同等效为碎石墙进行平面应变条件下分析，而实际工况的响应具有明显三维空间效应[13]。Tang等[10]对Adalier等[3]试验中的碎石桩及周边土体按照对称性选取单桩单元区采用三维砂土液化模型进行了响应分析，碎石桩的群桩效应及加固区与非加固区的响应差异在分析中未能反应出来。牛琪瑛等[12]以碎石桩加固饱和砂土震动台试验为基础进行了相应的系列试验和数值模拟分析，因限于1g震动条件下的模型比尺关系而仅能反应碎石桩加固在地表浅层的抗液化效果。
本文研究为渗透系数相对较高(约10-5 m/s～10-3 m/s)的可液化砂土、粉质砂土、或砂质粉土场地。采用砂土液化大变形模型模拟可液化饱和砂土，采用等效非线性增量模型模拟碎石桩，应用FLAC3D有限差分软件，本文对Brennan完成的碎石桩加固约19 m厚饱和砂土场地的离心动力模型试验[5]进行三维数值模拟与抗液化效果分析，在此基础上，对碎石桩的加密效应及不同震动强度下加固可液化场地的动力响应进行模拟研究，重点关注在地震作用下可液化场地经碎石桩加固复合地基抗液化效果，超静孔压的累积与消散，初始小变形到液化大变形的变形发展，应力与应变变化，及加固区内部与外部的响应差异等。对渗透系数相对较低(约10-8m/s～10-5 m/s)的可液化粉土或砂质粉土场地需另考虑细粒含量与渗透系数对孔压消散和应力-应变响应的影响。
1 砂土液化模型与碎石材料模型
张建民等[14－17]针对地震液化大变形问题，通过对震害调查和室内试验现象的认识，从基本规律、物理机制、本构描述，到实用分析，提出了一个饱和砂土震动液化大变形的本构理论和预测方法。张建民[15]基于排水循环扭剪试验现象将砂土的体应变εv分解为有效球应力变化引起的体应变εvc和剪切引起的体应变εvd并进一步分解为可逆性剪切体应变εvd,re与不可逆性剪切体应变εvd,ir，揭示了循环加载条件下砂土液化后剪应变与三个体应变分量之间的内在联系。基于试验现象研究，张建民和王刚[16－17]发展了定量描述砂土液化后零有效应力状态的剪应变与体应变关系，在此基础上发展建立了可以描述饱和砂土液化后大变形的弹塑性本构模型。王睿和张建民等[18]在对砂土液化大变形机理充分认识的基础上，引入临界状态变量，符合三维应力空间中边界面和剪胀面上的应力映射规则，建立了三维砂土液化大变形本构模型。该模型能够实现对不同状态砂土单调和循环加载以及液化前后力学行为的统一描述。
本文利用该三维砂土液化大变形本构模型模拟可液化土体。邹佑学等[19]将该模型在FLAC3D中二次开发，对不排水循环扭剪、不排水/排水三轴压缩、循环三轴、及三维离心动力试验模拟分析，模拟与试验结果符合度均较高，验证了模型模拟砂土从液化前小变形到液化后大变形方面的能力和程序在处理三维计算问题上的有效性。该模型在OpenSees平台上的可液化地基单桩基础和群桩基础动力分析中也得到较好应用[20―24]。模型公式完整描述请参见文献[18―20]。
碎石材料模型考虑在动荷载作用下材料模量随围压与剪应变变化而变动，并考虑模型适用性，采用等效非线性模型并根据计算程序的需要表达为其增量形式[13]。根据Seed等[25]和Rollins等[26]对碎石剪切模量的建议，碎石材料初始最大等效剪切模量Gmax表示为：
式中：
/kPa为平均有效应力；(K2)max为与相对密度或剪切波速相关的模量系数，约为80～180[25]；系数α=21.7；pa为标准大气压。
动剪切模量比[25―26]和以Sigmoidal函数表示的剪切模量衰减曲线[27]表示为：
式中：G为动剪切模量(割线模量)；γ为动剪应变；Ms为剪切模量比；a、b、x0和y0为拟合参数。
将应力-应变关系骨干曲线写成增量形式为：
式中：Gt为切线动剪切模量；τ为平均剪应力；dτ为剪应力增量；dγ为动剪应变增量。
体积模量表示为：
式中，ν为碎石材料小应变泊松比。
将上述碎石材料等效非线性模型增量形式，并采用摩尔-库仑强度准则描述屈服，进行二次开发植入FLAC3D程序[27]用于模拟碎石桩。
2 离心动力试验描述
为验证所采用模型与方法分析碎石桩加固可液化场地地震动力响应及抗液化效果的有效性，本文对 Brennan[5]在剑桥大学进行的离心动力模型试验 AJB-8a进行模拟对照分析。试验的离心加速度为50 g，下文涉及到的相关尺寸等参数均用原型比尺表述。试验原型为深约 19 m的均匀饱和砂土与碎石排水桩群，碎石桩埋深至砂层底部，试验布置和测点如图1所示。碎石桩共17根，桩径1.2 m，呈三角形布置，桩间距 3.0 m。砂土为颗粒级配 E级的硅砂(silica sand)，相对密度为37.4%(e=0.864)，碎石排水桩材料为颗粒级配B级硅砂，平均粒径约比E级大6.4倍，渗透系数分别为1.0×10-4 m/s和0.7×10-2 m/s，材料参数如表1所示。试验模型建在叠环型模型箱内，能够对在地震荷载下自由场的响应进行模拟。为匹配离心试验的动力与渗流时间比尺，采用粘滞系数为50cSt的硅油作为孔隙饱和液体。从模型底部水平输入的强震动持续时间约25 s，峰值加速度为 0.163 g，震动频率约 0.82 Hz，记录的模型箱底部振动时程如图2所示。模型共布置了7个加速度传感器和16个孔压传感器。
由于试验中所采用硅砂的材料参数不足够用于标定拟采用的砂土液化大变形本构模型的模型参数，参照文献[8]采用Toyoura砂模型参数模拟试验 ABJ-8a的硅砂在分析砂土液化时取得较好的效果，本文在模拟离心试验的级配E级硅砂时，本构模型参数也采用经过材料试验和单元试验验证过的Toyoura砂(丰浦砂)模型参数。Toyoura砂和硅砂的颗粒特征与摩擦角接近，硅砂摩擦角φcr约32o，Toyoura砂在不同文献中φcr约31.2～33.5o，材料参数对比如表1所示。碎石排水材料采用等效非线性增量模型模拟，考虑碎石排水材料在强地震作用下的非线性力学行为，及作为桩周砂土所产生超静孔隙水压的消散通道，不考虑其产生超静孔隙水压。
图1 离心模型试验AJB-8a布置[5,8]
Fig.1 Layout of centrifuge test AJB-8a[5,8]
表1 硅砂和Toyoura砂的材料参数[5,8,17]
Table 1 Material properties of Silica sand and Toyoura sand[5,8,17]
图2 试验AJB-8a记录和滤波后的输入加速度时程[5]
Fig.2 Recorded and filtered base motion during test ABJ-8a[5]
3 分析模型
3.1 计算模型和材料参数
模拟ABJ-8a试验的平面布置和网格划分如图3所示，长 33.6 m，宽 12.6 m，高 19 m，共 23712个单元(zones)。试验中的碎石桩截面为圆形，桩径d为1.2 m，桩间距S为3.0 m。计算中按照桩截面面积等效为正方形桩分析，边长为1.07 m，桩中心之间距离 3.0 m。因试验中饱和溶液采用了 50 cSt
的硅油，计算中砂土和碎石桩的渗透系数采用原型值1.0×10-4 m/s和0.7×10-2 m/s进行分析。砂土模型计算参数见表2[18―20]；碎石材料模型参数见表3，与Rollins等[26]统计分析了多项研究中980个碎石材料试验数据的碎石模量衰减曲线相比较，
曲线拟合如图4所示，在动剪应变较大(γ＞0.5%)时趋于平缓，以提高计算稳定性和避免局部区域因剪切模量衰减过多引起变形异常。
图3 模型网格布置及三维视图
Fig.3 Schematics of model layout and 3D isometric view
图4 碎石材料G/Gmax～γ曲线(图片来源：Rollins等[26])
Fig.4 G/Gmax～γ curve for the gravel material (sourced form Rollins et al [26])
模型的力学和流体边界条件为：1)底部固定，地震动从底部沿x轴方向水平输入；2)考虑叠环型模型箱模拟自由场条件，边界左侧和右侧对应高度节点的自由度进行捆绑约束；3)侧向边界在y向固定约束；土层表面自由；4)底部和周边为不排水边界，顶面为自由排水边界，设置固定孔压为0。
在模拟分析离心试验 AJB-8a的基础上，本文继续开展了不同砂土相对密度与地震动输入强度等6个模型参数分析，见表4。
表2 饱和砂土材料模型参数[18―20]
Table 2 Model parameters for saturated sand[18―20]
表3 碎石材料模型参数
Table 3 Model parameters for stone columns
表4 计算模型及参数
Table 4 Simulation Model and parameters
注：αmax为输入震动的最大加速度值。
3.2 地震波输入
模拟 ABJ-8a试验时，采用模型箱底记录的输入波，并经SeismoSignal程序进行滤波和基线校正后的波形，如图2所示。
模型参数分析时，地震动输入选用Loma Prieta地震CDMG STATION 47381台站记录的0 s～26 s波形，其加速度时程经滤波和基线校正后的输入波形(加速度、速度和位移时程)如图5所示，最大加速度αmax约 0.37 g，出现在 5.53 s附近；模型SC-4～SC-6的输入波形按αmax值参数比例调整。计算时叠加阻尼比为0.01的瑞利阻尼，以调整低应变幅振动时阻尼偏小和降低高频噪声影响。
3.3 计算步骤
先进行弹性、塑性和渗流计算生成初始应力场与孔压场，将位移和速度置为零，再启动动力计算，输入震动进行三维弹塑性流固耦合动力响应及震后固结分析。
图5 输入震动波形
Fig.5 Input excitations after baseline adjustment
4 模拟结果与分析
4.1 离心动力试验模拟结果
图6所示为模拟计算的加固区内部和外部不同位置、不同深度(约 5 m、10 m、18 m)的砂土中超静孔压比ru时程，与离心试验记录的对比。ru为超静孔压Δu与初始竖向有效应力
之比。试验和模拟显示在 10 m以上相应位置都达到了初始液化(ru=1.0)，而在PPT6673位置(约深18 m)超静孔压比达到约 ru ≈ 0.84，未达到初始液化。模拟值在浅层部位(约5 m以上)与试验值超静孔压产生和液化基本同步，液化从浅层逐步向深层发展，与试验相符，而试验向深部液化发展时速度稍快。
图6 超静孔压比的试验与数值模拟对比
Fig.6 Experimental and computed the excess pore pressure ratio at various depths and locations of test AJB-8a
图7为模拟计算的加固区内部和外部不同位置对应测点的加速度响应时程与试验记录的对比，二者总体上符合较好。模拟计算和试验都反应出了土体在液化后的加速度响应衰减，和深层加速度向上的有效传递。对液化后浅层土体，试验的加速度保持了小幅度的震动，模拟计算加速度震动幅度较小，衰减更明显，保持有小幅的液化砂土往返活动性。
图7 加速度响应试验与数值模拟对比
Fig.7 Experimental and computed the ground acceleration at various depths and locations of test AJB-8a
4.2 碎石桩加密效应及不同震动强度的动力响应
在离心动力试验验证所用模型和方法模拟可液化场地碎石桩复合地基有效性的基础上，对砂土相对密度Dr为50%/60%/70%的工况SC-1、SC-12、SC-13进行了分析，研究碎石桩加密效应的抗液化效果，图8所示为不同位置和深度超静孔压比响应对比。从图可见，加固区p1点和加固区外p9点(距桩边约2.2 m)，在5 m深度处SC-1和SC-2均达到液化，而SC-3则在p9点达到液化状态，p1点未达到液化。SC-1的10 m深度处p1和p9点都达到液化，在震动后期p1点孔压即开始消散而p9点则尚未减少，而SC-2和SC-3在两个点都未达到液化但p1点孔压比明显低于p9点。加固处内距桩边较近的p5点，在10 m深处SC-1、SC-2、SC-3都没有达到液化，其孔压相应低于 p1点孔压。加固区外边缘的p8点(距桩边约0.5m)，在15 m深度三种工况均未达到液化。分析可见：1)加固区内抗液化效果与加固区外相比，得到改善；2)加固区内土体在强震的作用下浅层可能会发生液化，碎石桩对加固区内距桩较近部位抗液化效果改善明显优于距桩远的部位，非加固区为密实砂土(Dr=70%)情况下，在强震作用时浅层也可能会发生液化；3)碎石桩的加密效应可显著改善可液化场地的抗液化效果；4)在同样深度发生液化的情况下，加固区内的孔压消散明显快于加固区外，可更快地离开液化状态。
图8 碎石桩加固不同相对密度砂土的超静孔压响应对比
Fig.8 Comparison of the excess pore pressure ratio for different relative density sands improved by SC
图9所示为SC-2工况下，h=5 m深度和y=6.3 m剖面在t=15 s及30 s时的超静孔隙水压力等高线示意图。从图9(a)与图9(c)可见，加固区外约2 m～3 m范围内孔压梯度较陡而等高线分布密集，加固区内则孔压较低且梯度分布较缓，碎石桩加固明显降低了超静孔压累积，有益于改善加固区内土体抗液化效果。在t=15 s震动较强时桩周附近孔压及其梯度大于在t=30 s时，碎石桩周边的土体孔压在强震过后震动持续时就在快速消散。从图9(b)与图9(d)可见：1)碎石桩排水的抗液化有效影响范围至加固区外约2.5倍～3倍桩径，在浅层小而随深度增加；2)在深层部位孔压有较大幅度消散时，浅层部位消散较少，深层部位孔压消散快于浅层部位；3)中间桩区域的孔压消散明显快于外围桩区域的孔压消散；加固区外的超静孔压向碎石桩消散，外围桩发挥着加固区排水屏障和非加固区排水通道的作用。由此，在碎石桩加固处理可液化场地时，增强围护桩的排水效果，如：加大桩径或加密布桩及保持其排水通道畅通，有利于提高加固区的整体抗液化效果。
图9 碎石桩加固场地 t= 15s/30s 的孔压云图(模型SC-2)
Fig.9 Contour of the pore water pressure for model SC-2 at t=15s/30s
在震动强度amax分别为 0.37 g、0.30 g、0.25 g和0.20 g的SC-2、SC-4、SC-5、SC-6工况下，加固区内部和外部不同位置、不同深度的超静孔压比时程曲线如图10所示。从图10可见，在5 m深度处 p9点 SC-2、SC-4、SC-5工况均达到液化，且SC-2、SC-4工况液化状态持续至震动结束，而SC-5工况下孔压在震动持续时即缓慢降低，液化向浅层消退。在5 m深度处p1点SC-2、SC-4工况达到液化并在震动结束前孔压消散离开液化状态，其他工况则未达到初始液化。在10 m深度，各工况在p9、p1、p5点均未达到液化，且明显加固区内的孔压比低于加固区外的孔压比，加固区内距离碎石桩近位置的孔压比低于距桩远位置的孔压比。在 SC-6工况下，5 m及以下均未达到初始液化；从计算结果，SC-6工况加固区外的p9和p3点液化深度约为4 m，加固区内除桩间土浅层1 m～2 m的少量部位外基本未出现土体液化。分析可见，地震动强度对碎石桩加固可液化场地的超静孔压影响显著，震动强度越大，各位置的超静孔压比越大，土层液化越深。
图10 不同震动强度下超静孔压响应对比
Fig.10 Comparison of the excess pore pressure ratio under shaking with different vibration intensity
场地地表p1和p2点在不同工况下的加速度响应时程如图11所示。非加固的自由场p2点地表位置在SC-1～SC-6六个工况都达到了液化，在初始液化后加速度响应衰减明显，即土体抗剪切能力衰减明显。加固区内p1点的地表位置在SC-1～SC-6六个工况也都达到液化或短暂接近液化，但加速度保持着一定的震动幅度，土体相对密度越大加速度衰减越少，震动强度较低时加速度的衰减较少。碎石桩加固虽然未能阻止 p1点地表位置的液化，但因提高加固区的整体刚度、改善下层土体受力性能，提高了其液化势或抗液化能力，因此 p1点地表相对 p2点在液化后保持有更强的往返活动性与震动加速度响应，土体抗剪切能力得到加强。
图11 地表加速度响应时程对比
Fig.11 Comparison of the surface acceleration responses
加固区和非加固区的不同位置5 m深度处的剪应力-剪应变(τxz～γxz)响应对比如图12～图14所示。从图可见，土体的大剪应变发生在初始液化后的零剪应力状态与往返活动性时的低剪应力状态；未固区p3和p9点的剪应变响应大于加固区内p1点的剪应变响应，表明碎石桩加固有助于降低可液化土体的剪应变响应，且对其边缘的非加固区土体的剪应变响应降低也有所帮助；砂土相对密度较大工况的剪应变响应小于相对密度较小工况的剪应变响应，而剪应力则相对密度越大产生的剪应力响应越大。SC-6工况在5 m深度未出现初始液化，处于小应变阶段，与图10超静孔压比所示相符合。
图12 剪应力-剪应变响应(p1：h=5 m)
Fig.12 Calculated stress-strain relations (p1：h=5 m)
图13 剪应力-剪应变响应(p9：h=5 m)
Fig.13 Calculated stress-strain relations (p9：h=5 m)
图14 剪应力-剪应变响应(p3：h=5 m)
Fig.14 Calculated stress-strain relations (p3：h=5 m)
在SC-1、SC-3、SC-2、SC-6等工况不同高度、不同位置的剪应力-有效球应力
路径的模拟结果对比如图15～图17所示。未加固区自由场p3点在5 m深处SC-1、SC-3、SC-2工况均达到零剪应力的液化状态，土体相对密度较高的工况经过更多次的震动循环才达到初始液化，而 SC-6工况因震动强度较低未达到初始液化。加固区 p1点在5 m深处SC-1、SC-2达到液化，且与p3点相比，经过多几次震动循环才达到初始液化状态，而SC-3/6则未达到液化。在加固区p5点的10 m深处4个工况均未达到液化。从图15～图17可见：1)土体相对密度高的工况中剪应力响应峰值较大，土体的剪应力比峰值较大，即剪应力与有效球应力之比的峰值较大；2)SC-2与SC-6工况对比显示在土体相对密度相同时，震动强度较高时土体中剪应力响应峰值较大；3)在相同部位，土体相对密度较大、震动强度较小的工况中有效球应力的衰减较少。
图15 应力路径(剪应力-有效球应力p3：h=5 m)
Fig.15 Calculated effective stress paths (p3：h=5 m)
图16 应力路径(剪应力-有效球应力p1：h=5 m)
Fig.16 Calculated effective stress paths (p1：h=5 m)
图17 应力路径(剪应力-有效球应力p5：h=10 m)
Fig.17 Calculated effective stress paths (p5：h=10 m)
加固区和非加固区的地面在不同位置x向位移时程及震动输入位移时程的比较如图18所示。从图18可见，在震动前期阶段，地面不同点基本上随基底输入震动而同步震动，侧向位移的振幅则有所放大，砂土相对密度较低时振幅放大较多，如p1点位置在 SC-1和 SC-3工况的侧向位移分别为24.15 cm和 22.51 cm，放大系数分别约 23.6%和15.2%。在震动过程进行且场地逐步接近液化或出现初始液化后，地表点的震动与‘’基底震动逐步出现偏离，不同位置震动偏离出现有先有后，偏离基底震动的分布则不同工况响应不统一。砂土较密或震动强度较低时，地表震动保持与基底震动相似的震动形态程度越高，此点与图11所示加速度响应相符合。
图18 侧向变形时程曲线
Fig.18 Experimental and computed the ground acceleration at various depths and locations of test AJB-8a
5 结论
采用砂土液化大变形本构模型模拟饱和砂土和等效非线性增量模型模拟碎石桩，应用 FLAC3D软件在模拟离心动力模型试验的基础上，进行碎石桩加固可液化场地复合地基三维液固耦合地震动力响应分析，结果和分析表明：
(1)所用模型与方法可合理地描述碎石桩加固可液化场地复合地基在地震作用下场地的动力响应特性和抗液化效果。模拟碎石排水桩加固可液化场地离心动力试验符合较好，验证了模型与方法的有效性及模拟能力。
(2)碎石桩加固区外约2 m～3 m范围内孔压梯度较陡、等高线分布密集，加固区内孔压较低且梯度分布较缓，碎石桩加固可明显改善可液化土体的抗液化能力，有效影响范围至加固区外约2.5倍～3倍桩径，在浅层小而随深度增加。
(3)碎石桩加固区的中间桩区域的孔压消散明显快于外围桩区域的孔压消散，加固区外的超静孔压向碎石桩消散，外围桩发挥着加固区排水屏障和非加固区排水通道的作用。在碎石桩加固处理可液化场地时，增强围护桩的排水效果，如：加密围护桩布桩或加大其桩径及保持排水通道畅通，有利于提高加固区的整体抗液化效果。
(4)碎石桩加固使可液化土体相对密度增大时，可提高土体的剪应力比峰值，增加土体抗剪切的峰值剪应力，降低剪应变响应，土体中的超静孔压累积趋缓而孔压消散速度趋快，场地达到初始液化的深度变浅。
(5)在震动过程中场地接近液化或出现初始液化后，地表点的震动与基底震动会逐步偏离，砂土较密或震动强度较低时，地表保持与基底相似的震动形态程度越高。在所分析的 0.2 g震动强度工况(SC-6)下，碎石桩加固区除桩间土浅层1 m～2 m少量部位外整体上未出现土体液化，基本达到消除场地液化的加固目的。
[1]Seed R B,Cetin,K O,Moss R,et al.Recent advances in soil liquefaction engineering and seismic site response evaluation [C]// Shamsher Prakash,Proceedings of 4th International Conference on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics,San Diego：University of Missouri-Rolla,2001(SPL-2)：1－45.
[2]Seed H B,Booker J R.Stabilization of potentially liquefiable sand deposits using gravel drains [J].ASCE Journal of Geotechnical Engineering Division,1977,107(7)：757－768.
[3]Adalier K,Elgamal A.Mitigation of liquefaction and associated ground deformation by stone columns [J].Engineering Geology,2004,72(2004)：275－291.
[4]Hausler E A.Influence of ground improvement on settlement and liquefaction：a study based on field case history evidence and dynamic geotechnical centrifuge tests [D].Berkeley：University of California,2002.
[5]Brennan A J.Vertical drains as a countermeasure to earthquake-induced soil liquefaction [D].UK：University of Cambridge,2004.
[6]Adalier K.Elgamal A,Meneses J,Baez J I.Stone columns as liquefaction countermeasure in non-plastic silty soils [J].Soil Dynamics and Earthquake Engineering,2003,23(7)：571－584
[7]Badanagki M,Dashti S,Kirkwood P.Influence of dense granular columns on the performance of level and gently sloping liquefiable sites [J].Journal of Geotechnical and Geoenvironmental Engineering,2018,144(9)：04018065-1―04018065-14.
[8]Papadimitriou A,Moutsopoulou M E,Bouckovalas G,Brennan A.Numerical investigation of liquefacton mitigation using gravel drain [C]// Kyriazis Pitilakis,Proceedings of 4th International Conference on Earthquake Geotechnical Engineering,Thessaloniki,Greece：International Society of Soil Mechanics and Geotechnical Engineering,2007(paper No.1548)：1―12.
[9]Li P,Dahti S,Badanagki M,Kirkwood P.Evaluating 2D numerical simulations of granular columns in level and gently sloping liquefiable sites using centrifuge experiments [J].Soil Dynamics and Earthquake Engineering,2018,110(2018)：232－243.
[10]Tang L,Zhang XY,Liang XZ.Numerical Simulation of centrifuge experiments on liquefaction mitigation of silty soils using stone columns [J].KSCE Journal of Civil Engineering,2016,20(2)：631－638.
[11]Rayamajhi D,Nguyen T V,Scott A A,et al.Numerical study of shear stress distribution for discrete columns in liquefiable soils [J].ASCE Journal of Geotechnical and Geoenvironmental Engineering,2014,140(3)：1－9.
[12]牛琪瑛,郭英,吴永娟,等.碎石桩加固不同密实度液化土的孔隙水压力变化规律探讨[J].工程力学,2010,27(增刊 I)：159－163.Niu Qiying,Guo Ying,Wu Yongjuan,et al.Study on pore water pressure variation of different density liquefiable sand soil reinforced by gravel pile [J].Enigineering Mechanics,2010,27(Suppl I)：159－163.(in Chinese)
[13]邹佑学,王睿,张建民.可液化场地碎石桩复合地基地震动力响应分析 [J].岩土力学,2019,40(6)：2043－2055.Zou Youxue,Wang Rui,Zhang Jianmin.Analysis for the seismic response of stone columns composite foundation in liquefiable soils [J].Rock and Soil Mechanics,2019,40(6)：2043－2055.(in Chinese)
[14]张建民.砂土动力学若干基本理论探究[J].岩土工程学报,2012,34(1)：1－50.Zhang Jianmin.New advances in basic theories of sand dynamics [J].Chinese Journal of Geotechnical Engineering,2012,34(1)：1－50.(in Chinese)
[15]张建民.砂土的可逆性和不可逆性剪胀规律[J].岩土工程学报,2000,22(1)：12－17.Zhang Jianmin.Reversible and irreversible dilatancy of sand [J].Chinese Journal of Geotechnical Engineering,2000,22(1)：12－17.(in Chinese)
[16]张建民,王刚.砂土液化大变形的机理[J].岩土工程学报,2006,28(7)：835－840.Zhang Jianmin,Wang Gang.Mechanism of large post-liquefaction deformation in saturated sand [J].Chinese Journal of Geotechnical Engineering,2006,28(7)：835－840.(in Chinese)
[17]Zhang J M,Wang G.Large post-liquefaction deformation of sand,part I：physical mechanism,constitutive description and numerical algorithm [J].Acta Geotechnica,2012,7(2)：69－113.
[18]Wang Rui,Min Jian Zhang,Wang Gang.A unified plasticity model for large post-liquefaction shear deformation of sand [J].Computers and Geotechnics,2014,59(2014)：54－66.
[19]邹佑学,王睿,张建民.砂土液化大变形模型在FLAC3D中的开发与应用[J].岩土力学,2018,39(4)：1525－1534.Zou Youxue,Wang Rui,Zhang Jianmin.Implementing a plasticity model for large post-liquefaction deformation of sand into the FLAC3D program [J].Rock and Soil Mechanics,2018,39(4)：1525－1534.(in Chinese)
[20]Wang R.Single Piles in Liquefiable Ground：Seismic response and numerical analysis methods [M].Berlin：Springer,2016.
[21]Wang R,Fu P,Zhang JM.Finite Element model for piles in liquefiable ground [J].Computers and Geotechnics,2016,72：1－14.
[22]王睿,张建民.可液化地基中单桩基础的三维数值分析方法及应用[J].岩土工程学报,2015,37(11)：1979－1985.Wang Rui,Zhang Jianmin.Three-dimensional elastic-plastic analysis method for piles in liquefiable ground [J].Chinese Journal of Geotechnical Engineering,2015,37(11)：1979－1985.(in Chinese)
[23]刘星,王睿,张建民.液化地基中群桩基础地震相应分析[J].岩土工程学报,2015,37(12)：2326－2331.Liu Xing,Wang Rui,Zhang Jianmin.Seismic response analysis of pile groups in liquefiable foundations [J].Chinese Journal of Geotechnical Engineering,2015,37(12)：2326－2331.
[24]Wang R,Liu X,Zhang J M.Numerical analysis of the seismic inertial and kinematic effects on pile bending moment in liquefiable soils [J].Acta Geotechnica.2017,12 (4)：773－791.
[25]Seed H B,Robert T W,Idriss I M,et al.Moduli and damping factors for dynamic analyses of cohesionless doils [J].Journal of Geotechnical Engineering(ASCE),1986,112 (11)：1016－1032.
[26]Rollins K M,Evans M D,Diehl N B,et al.Shear modulus and damping relationships for gravels [J].Journal of Geotechnical Engineering(ASCE),1998,124(5)：396－405.
[27]Itasca Consulting Group Inc..Fast Language Analysis of continua in 3 dimensions,version 5.0,user's manual.Itasca Consulting Group Inc.,2012.