饱和砂土自由场地地震液化的离散元模拟
史旦达+王飞
摘要：为研究饱和砂土自由场地的地震液化过程，基于颗粒离散元与计算流体动力学耦合方法，在PFC3D程序中通过二次开发实现对自由场地基底不规则地震波的输入，研究在实测Kobe地震波激励下该液化过程的宏细观力学响应，并与已有的福建中细砂离心机模型试验结果进行对比。数值模拟结果表明：试样发生液化的时刻对应于基底地震加速度峰值时刻，浅部土层先于深部土层发生液化，超静孔隙水压力的消散是从深部土层逐渐向浅部土层发展的；伴随液化过程，各深度土层的有效应力路径均从初始应力状态向应力空间的原点移动；数值试验呈现出的上述规律均与实际砂土离心机模型试验呈现出的规律保持一致。数值试样宏观的液化特征在微观机理上可由流体对颗粒的归一化拖曳力、局部孔隙比、平均配位数、颗粒位移的变化规律来相应表征。
关键词：
饱和砂土； 自由场地； 地震荷载； 液化； 离散元模拟
中图分类号： TU441.5
文献标志码： A
Discrete element simulation on seismic liquefaction
of saturated freefield sand deposits
SHI Danda1， WANG Fei2
（1. College of Ocean Science and Engineering， Shanghai Maritime University， Shanghai 201306， China；
2. East China Architectural Design & Research Institute， Shanghai 200070， China）
Abstract：
In order to research the seismic liquefaction of saturated freefield sand deposits， using the coupled method of the particle discrete element method and computational fluid dynamics， both macro and micro mechanical responses of the liquefaction are explored subjected to the real Kobe earthquake wave excitation， where the input of irregular seismic waves on freefield base is realized by selfdeveloped procedures in PFC3D. The numerical results are compared with the centrifugal model test results on Fujian sand， which is a uniformly graded medium to fine sand. It is found that the exact instant for sample liquefaction corresponds to the peak time of the base seismic acceleration. Liquefaction happens in shallow sand layers prior to deep sand layers， whilst the excess pore water pressure dissipates earlier in deep layers than in shallow layers. During liquefaction， the effective stress path of each sand layer moves from the initial stress condition to the coordinate origin in the stress space. The above laws obtained by numerical tests are consistent with the laws obtained by the real sand centrifugal model test. The macroscopic liquefaction behaviors are correspondingly characterized by the normalized fluid drag force on particles， the porosity， the average coordination number and the particle displacement in the micro mechanism.
Key words：
saturated sand； free field； earthquake loading； liquefaction； discrete element simulation
0引言
在地震災害中，饱和砂土地基的液化是造成大量建筑物震害的重要原因。关于砂土液化宏观机理的研究，自SEED等[1]于1966年提出“初始液化”概念以来已近半个世纪，积累的成果较多[23]。然而，关于砂土液化细观机理的研究，由于试验条件和研究手段等方面的制约，目前成果相对较少。离散元方法（Discrete Element Method， DEM）由于制样便捷、细观组构信息丰富且容易获取等优点，在粒状土静、动力特性及其机理研究方面发挥了重要作用。
目前采用DEM模拟砂土液化通常有两种方法。一种称为常体积法，SITHARAM[4]、GONG等[5]、史旦达等[6]、SUZUKI等[7]均采用该方法模拟了砂土的振动液化特性，该方法的优点是计算效率高，缺陷是试样本身未设置流体单元，不能考虑固液二相细观相互作用。另一种为流固耦合方法，即颗粒离散元与计算流体动力学（Computational Fluid Dynamics， CFD）耦合，HAKUNO等[8]最早提出可用于饱和砂土液化分析的流固耦合方法，但其提出的在孔隙尺度上的流体效应模拟计算效率很低。为提高计算效率，NAKASA等[9]对HAKUNO等提出的方法进行了优化，流体采用比孔隙尺度稍大的流体网格来模拟，每个流体网格内能够容纳若干个颗粒，流体流速与压力梯度之间满足Darcy定律。ZEGHAL等[10]在前人研究基础上提出了一种更为实用的耦合方法，流体运动满足均一化NavierStokes方程，流体对颗粒的作用力采用已有的半经验公式表示，目前应用较广的离散元程序PFC3D提供的流固耦合分析模块就采用了这种方法。刘洋等[11]基于PFC3D流固耦合模块，研究了饱和砂土的液化响应，但局限于规则的简谐荷载输入，并未考虑实际地震荷载，数值模拟结果也未与物理试验结果展开比较。
本文基于PFC3D流固耦合模块，通过二次开发实现基底不規则地震波的输入，研究实测Kobe地震波激励下饱和砂土自由场地的液化宏细观力学响应，并与典型的离心机试验结果进行对比。除分析超静孔隙水压力、有效应力路径等宏观响应外，重点探讨流体对颗粒的拖曳力、孔隙比、平均配位数和颗粒位移等细观力学响应，用来揭示液化过程中宏观液化规律与细观液化机理之间的宏细观关联。
1流固耦合原理简介
1.1液相流体
假定液相流体不可压缩且密度不变，其运动满足均一化NavierStokes连续方程和动量方程[12]，
nt+·（nvf）=0 （1）
ρf（nvf）t+·（nvfvf）=·（nτf）-
fi+nρfg（2）
式中：n为孔隙率；ρf为流体质量密度；vf和τf分别为流体平均速度矢量和平均应力张量；fi为单位体积内平均流
体颗粒拖曳力矢量；g为重力加速度；
为哈密顿算子。
1.2流体颗粒拖曳力
流体颗粒相互作用拖曳力可采用Ergun半经验公式[13]计算，
fi=150μf（1-n）2n·dp2（vf-vp）+
1.75（1-n）ρfvf-vpdp（vf-vp） （3）
式中：μf为流体黏滞系数；dp为流体单元内的颗粒平均粒径；vp为颗粒平均速度矢量。
1.3固相颗粒
固相颗粒采用离散元模拟，其运动满足牛顿第二定律，且考虑孔隙流体的作用，单个颗粒的平动和转动方程分别为
mpv·p=mp·g+cfc+fd
（4）
Ipω·p=crc×fc （5）
式中：mp和Ip分别为颗粒的质量和惯性矩；vp和ωp分别为颗粒的平动和转动速度矢量；fc为接触点处粒间作用力矢量；rc为颗粒中心至接触点的方向矢量；fd为流体施加于颗粒的作用力，包括浮力项和流体颗粒拖曳力项。
fd=-pf+fi1-nVp（6）
式中：Vp为颗粒体积；pf为孔隙流体平均压力梯度。
2试样制备及地震荷载输入
2.1数值试样制备
在PFC3D中可以通过对数值试样施加一定倍数的重力加速度来还原现场原型的场地条件，数值模型尺寸为140 mm（长）×60 mm（宽）×250 mm（高），对数值试样施加30倍重力加速度，还原至现场原型的场地尺寸为4.2 m（长）×1.8 m（宽）×7.5 m（高）。本文涉及的土层深度和分析得到的物理力学量值均为换算至场地原型的数值，换算方法满足离心试验相似性原理[14]。试样制备采用重力沉积法，颗粒为三维圆球颗粒，颗粒粒径范围5～10 mm，颗粒平均粒径为7.65 mm，初始孔隙比为0.727，参考福建标准砂的最大孔隙比emax=0.848和最小孔隙比emin=0.519，计算得到数值试样的相对密度约为37%。生成的颗粒数量为5 198个，数值试样如图1a）所示，试样顶部为自由面，底部为刚性墙，侧面采用周期边界，地震荷载通过底部（基底）墙体输入。
试样内布置固定位置的流体网格，如图1b）所示。流体网格在X，Y，Z三个方向上的尺寸分别为28 mm，20 mm和21 mm，每个流体网格内约有颗粒32个。试样内部设置4个量测圈，用于监测孔隙比、平均配位数等细观组构量变化。量测圈的设置见图1a），量测圈布置在试样竖向中心轴位置，第1，2，3，4号量测圈的中心位置距试样顶面的高度分别为0.13h，0.29h，0.54h和0.79h，h为试样高度。第1，2，3，4号量测圈分别对应于第1，2，3，4层土层，且分别对应于图1b）中Z方向第11，9，6，3号流体网格。
数值试样中颗粒和流体的细观参数设置见表1。颗粒之间的接触采用线性接触模型，流体黏滞系数为5.02 Pa·s。DEM采用中心差分格式，计算得到的平均时步为2.5×10-6 s。为保证计算效率，每执行100步离散元计算后执行1步CFD计算，即CFD计算时步为DEM计算时步的100倍，为2.5×10-4 s。
2.2地震荷载输入
地震荷载通过基底墙体输入，采用实测Kobe地震波（见图2a）），加速度峰值为0.32g，持续时间为50 s。数值模拟中在基底墙体不能直接输入加速度，需要通过对加速度时程曲线积分获得速度时程曲线（见图2b）），再通过给基底墙体赋予速度来实现。对模型施加的加速度为30g，根据相似性原理，模型基底激振持续时间的相似比为1/30；原型激振
持续时间为50 s，换算至模型激振持续时间即为1.67 s。在给基底墙体输入速度的过程中，在速度时程曲线输入文件中每隔0.000 4 s读取一个数据点（选择较小的时间间隔能确保输入速度曲线的精度），通过PFC3D内嵌的Fish语言编程将速度值赋给基底墙体，完成基底地震荷载输入。
3超静孔隙水压力及归一化拖曳力分析
3.1超静孔隙水压力
超静孔隙水压力比ru，即为激振引起的试样超静孔隙水压力u与初始竖向有效应力之比，其数值大小及随时间变化的规律已被广泛用于试样液化分析中。本文数值模拟中，激振引发的u可以通过监测流体网格中心节点处的流体压力得到。4个深度土层的监测结果见表2，4个深度土层的ru随时间t的变化曲线见图3a）。
a）数值模拟结果
b）离心机模型试验结果[17]
图3数值模拟与离心机模型试验超静孔压比曲线对比
分析图3a）并结合表2可知：（1）整个试样发生液化的时刻约为10 s时，液化发生的时刻与输入的地震加速度峰值时刻（见图2a））保持一致；浅部土层ru峰值最大，随着深度的增加ru值略有减少，但数值上仍超过0.9，可以认为整个试样均发生了液化；从液化发生的时刻看，浅部土层先发生液化，深部土层后发生液化；已有的室内试验和液化现场实例调查均表明[1516]，对于饱和砂层的液化分布，随着埋深的增加，土层所受的初始竖向有效应力逐渐增加，因此深部土层较难达到初始液化。（2）第10～18 s，各深度土层的ru均出现迅速下降，u消散较快。（3）大概25 s后，试样的u消散过程趋于稳定，残余ru随土层深度的增加而减小，第1层土层的残余ru约为0.48，到第4层土层时，其残余ru减小为0.16左右。
梁孟根等[17]进行了饱和砂土自由场地液化响应的离心机振动台试验，试验砂料为福建中细砂，模型箱内的试样尺寸为730 mm（长）×330 mm（宽）×400 mm（高），加速度为50g，基底输入德阳波250 s，加速度峰值为0.392g，加速度峰值出现的时刻为
20 s。试样内沿高度方向自上而下共布置5个孔压计（P1～P5），为了与数值试样4个土层深度位置尽量保持对应，选取P1，P2，P3和P5的量测结果进行对比分析（P1，P2，P3和P5孔壓计距试样顶面的深度依次为0.16h，0.33h，0.50h和0.85h，h为试样高度）。该离心机模型试验输入的为德阳波，本文数值模拟采用Kobe地震波，在加速度峰值等特性上存在一定差异，因此仅从定性规律角度对比分析离心机模型试验和数值模拟的u变化规律。
图4为福建中细砂与本文数值试样的级配对比曲线。分析图4可得：福建中细砂平均粒径d50=0.16 mm，不均匀因数Cu=1.6，曲率因数Cc=0.96；数值试样d50=8.1 mm，Cu=1.422和Cc=1.051。与福建中细砂相比，数值试样除
适当放大颗粒平均粒径外（减少一定的颗粒数量，保证数值模拟可以实施），Cu和Cc值均
与福建中细砂的相近，可以认为数值试样的级配分布基本与福建中细砂的一致。
图3b）给出了文献[17]离心机模型试验实测的ru随时间变化的曲线。对比分析图3a）与图3b）可得：（1）从各深度ru曲线的整体发展形态看，数值试样和离心试样均经历了u增长、试样发生液化和液化后u消散的过程。（2）数值试样和离心试样发生液化的时刻均对应于输入加速度峰值的时刻（离心机模型试验加速度峰值时刻为20 s），在ru逐渐增加并接近1.0的过程中，u均呈现振动式发展，试样发生了循环液化。（3）数值试样与试验试样均呈现深部土层u消散快，浅部土层u消散慢的规律，这在物理机理上与底部边界不排水而顶部自由面排水有关，鉴于这一特点，液化后土层中的u消散路径为自下而上，消散过程中深部土层的u必然会传递至浅层并造成浅层u的累积，因此试样会呈现出深部土层u消散快，浅部土层u消散慢的规律。（4）与离心试样相比，数值试样u消散过程更为明显，这是由于离心机模型试验输入的德阳波在加速度峰值时刻后仍有较高水平的加速度次峰值输入，导致u维持在较高的水平。
进一步对比分析数值模拟与离心机模型试验结果，将数值试样和试验试样在液化前（即u增长阶段）不同时刻u随深度的变化曲线分别绘于图5a）和图5b），在u振荡上升段，取一个振动周期内u的平均值。对比分析图5a）与图5b）可知，虽然数值试样与试验试样发生液化的具体时间不同，但数值试样表现出与试验试样一致的u变化规律：（1）在液化前，各深度u均随时间的增加而增大，直至达到初始竖向有效应力水平（即发生液化）；（2）浅部土层u增长较快，随着深度的增加u增长逐渐变慢，说明浅部土层先发生液化，深部土层后发生液化。
3.2归一化拖曳力
激振引起的试样宏观u在细观机理上对应于流体对颗粒的拖曳力大小。定义归一化拖曳力fs为颗粒受到的流体竖向拖曳力与该颗粒有效自重的比
值，当fs≥1，在细观上即表现为试样发生液化。图6给出了数值试样4个深度土层fs随时间的变化规律，fs<0表示颗粒受到的合力方向发生了改变，颗粒开始产生向下的加速度。
分析图6，并结合3.1节u的规律分析可知：（1）fs的变化规律对应于u的变化规律，fs出现峰值的时刻与试样的u出现峰值的时刻基本一致，不同深度土层的fs峰值均接近或略大于1，表明试样整体已发生液化。（2）在整个激振过程中，伴随u的增长和消散过程，fs经历了先增大、后减小，并逐渐趋于稳定的过程。（3）由第3.1节分析可知，u的消散是从深部土层向浅部土层发展的，因此试样液化后，深部土层fs越早进入残余稳定状态，其数值上越接近于0，而浅部土层由于受到深部土层向上排水产生的残余孔隙水压力累积效应的影响，在u消散过程中还经历了fs缓慢增大的过程，激振结束时（50 s）仍有一定数值的残余fs存在。
4其他监测结果分析
4.1有效应力路径
利用布置在试样内部的量测圈可以测得试样内部6个有效应力分量
σ′xx，σ′yy，σ′zz，σ′xy，σ′yz和σ′zx，进一步可计算得到平均有效主应力p′和广义剪应力q。
p′=13（σ′xx+σ′yy+σ′zz） （7）
q=12
（（σ′xx-σ′yy）2+（σ′yy-σ′zz）2+
（σ′zz-σ′xx）2+6（σ′xy+σ′yz+σ′zx））12（8）
图7绘出了数值试样不同深度土层在p′～q应力空间内的有效应力路径。分析图7可知：从激振开始到数值试样液化前，第1～4层土体的有效应力路径均从初始应力状态向p′～q应力空间的原点移动，应力路径自上而下行进，p′和q值均逐渐减小，并趋近于零；当试样达到液化状态后，随着u的循环振荡，试样发生循环液化，应力路径在原点附近来回折返，这一现象在浅部土层中表现得尤为明显。值得注意的是，无论应力路径如何折返，4个深度土层的应力路径均没有逾越试样的应力路径包线（即试样强度包线）。汪闻韶[18]研究了砂土液化与极限平衡和破坏之间的关联，指出砂土液化过程会出现循环活动性，应力路径会发生来回折返现象，但应力路径始终不会逾越砂样的极限平衡强度包线。因此，本文数值试样呈现的有效应力路径变化规律符合实际砂样的液化特性。
4.2颗粒位移与表面沉降
为分析液化过程中内部颗粒的运动情况，在距试样顶面1.1 m深度处选取标志颗粒进行追踪。该颗粒初始坐标为X=2.02 m，Y=0.84 m，Z=6.40 m，位于第1层土体内。图8a）绘出了整个激振过程中该标志颗粒的三维空间运动轨迹，图中实心圆点为颗粒运动的起始位置，空心圆点为颗粒运动的最终位置。分析图8a）可知：在激振过程中该标志颗粒的X方向坐标持续增加，说明试样发生了沿激振方向的累积水平位移；
在激振过程中该标志颗粒的
Y方向坐标的变化随机性较强，没有表现出一定的规律；
在激振过程中该标志颗粒的
Z方向坐标值经历了先增大后减小的过程，反映了液化过程中颗粒上浮、液化后随着u消散颗粒逐渐沉积的运动过程。
图8b）为整个激振过程中试样平均高度随时间的变化规律。由图8b）可知：在試样发生液化的时刻（10 s左右），试样平均高度达到最大，由于流体对颗粒拖曳力的作用，试样内的颗粒发生了较为明显的上浮，与初始高度相比整个试样高度增加了约0.17 m；之后，随着u的消散，试样的平均高度迅速减小，即试样发生了表面沉降，与初始时刻相比，最终的表面沉降量约为0.33 m，这一规律反映了饱和砂土自由场地的液化后震陷特征。
a）颗粒运动轨迹
b）试样平均高度变化
4.3孔隙比
孔隙比变化规律能够反映试样液化过程中的体积变化特征，利用设置在试样内部的量测圈（见图1a））可以测得4个深度土层孔隙比随时间的变化规律，见图9。分析图9可知：（1）各层土体的初始孔隙比略有差别，由于制样时的重力沉积效应，深部土层的初始孔隙比略小于浅部土层的。（2）在液化发生前，各层土体的孔隙比均有不同程度的增加，其中第3层和第4层土体孔隙比的增加较为明显，其原因与流体压力从基底输入导致深部颗粒的上浮较为明显有关。（3）当试样处于液化后u消散阶段时，各层土体的孔隙比均逐渐减小，与初始孔隙比相比，各层土体在液化后的最终孔隙比均小于初始孔隙比，表明试样在振动液化后变得更为密实。
4.4平均配位数
平均配位数的变化能够反映试样在液化过程中内部结构稳定性的演化，图10给出了4个深度土层平均配位数随时间的变化。
由图10可知：（1）试样发生液化时，各层土体的平均配位数均从初始时刻的4左右迅速降低至1左右，在这一过程中，浅部第1层土体平均配位数的降低最快，深部土层的平均配位数降低稍慢，其原因与浅部土层先液化、深部土层后液化的规律有关。（2）在u消散阶段，各层土体的平均配位数均得到了一定程度的恢复，深部第3层和第4层土体的平均配位数在激振开始30 s后趋于稳定，而浅部第1层和第2层土体的平均配位数在整个u消散过程中都在持续增加，因为深部土层u消散快、浅部土层u消散慢，所以此处反映的不同深度平均配位数的变化规律与宏观u的消散规律是对应的。（3）液化后，各深度土层的最终配位数均略大于初始配位数，这一规律也反映了试样振动液化后变密的特征。
5结论
基于PFC3D颗粒离散元与计算流体动力学耦合模拟方法，研究了真实地震荷载激励下饱和砂土自由场地液化过程的宏细观力学响应，并与已有的离心机模型试验结果进行了对比，得到的主要结论有：
（1）整个试样发生液化的时刻对应于基底地震加速度峰值时刻，因为深部土层的初始竖向有效应力高于浅部土层的，所以浅部土层先于深部土层发生液化，浅部土层的峰值超静孔隙水压力比ru也更大；液化后超静孔隙水压力u的消散先从深部土层开始，逐渐向浅部土层发展，浅部土层的残余ru大于深部土层的；数值试样呈现出的上述u变化规律均与福建中细砂自由场地地震液化的离心机模型试验结果一致。
（2）在试样液化过程中，各深度土层的有效应力路径均从初始应力状态向p′～q（p′为平均有效主应力，q为广义剪应力）应力空间的原点移动，p′和q值均逐渐减小并趋于零；当试样达到液化状态后，随着u的循环振荡，应力路径在原点附近来回折返，但应力路径不会逾越试样的极限平衡强度包线，这一现象符合实际砂样振动液化过程中的有效应力路径特征。
（3）对应于u的消散是从深部土层向浅部土层发展的宏观液化规律，在试样液化后，深部土层的残余归一化拖曳力更接近于0，且平均配位数也更早进入残余稳定状态，而浅部土层中仍存在一定数值的残余拖曳力作用，且平均配位数随着u的消散逐渐增加；对应于自由场地的宏观液化震陷特征，激振结束后各深度土层的最终孔隙比均小于初始孔隙比，最终配位数均略大于初始配位数，且试样整体发生了约0.33 m的表面沉降。
参考文献：
[1]SEED H B， LEE K L. Liquefaction of saturated sands during cyclic loading[J]. Journal of the Soil Mechanics and Foundations Division， 1966， 92（6）： 105134.
[2]莊海洋， 胡中华， 王瑞， 等. 饱和南京细砂初始液化后特大流动变形特性试验研究[J]. 岩土工程学报， 2016， 38（12）： 21642174.
[3]LU Xilin， HUANG Maosong. Static liquefaction of sands under isotropically and K0consolidated undrained triaxial conditions[J]. Journal of Geotechnical and Geoenvironmental Engineering， 2015， 141（1）： 04014087.
[4]SITHARAM T G. Discrete element modelling of cyclic behavior of granular materials[J]. Geotechnical and Geological Engineering， 2003， 21： 297329.
[5]GONG Guobin， THORNTON C， CHAN A H C. DEM simulations of undrained triaxial behavior of granular material[J]. Journal of Engineering Mechanics， 2012， 138（6）： 560566.
[6]史旦达， 周健， 刘文白， 等. 初始组构影响砂土液化势的细观数值模拟研究[J]. 水利学报， 2011， 42（7）： 766774.
[7]SUZUKI K， KUHN M R. Discrete element simulations of cyclic biaxial shear of a granular material with oval particles[J]. International Journal of Geomechanics， 2014， 14（3）： 06014005.
[8]HAKUNO M， TARUMI Y. A granular assembly simulation for the seismic liquefaction of sand[J]. Proceedings of Japan Society of Civil Engineers， 1988， 5（2）： 129138.
[9]NAKASA H， TAKEDA T， ODA M. A simulation study on liquefaction using DEM[C]//Earthquake Geotechnical Engineering. Rotterdam： Balkema， 1999， 637642.
[10]ZEGHAL M， El SHAMY U. A continuumdiscrete hydromechanical analysis of granular deposit liquefaction[J]. International Journal of Numerical and Analytical Methods in Geomechanics， 2004， 28： 13611383.
[11]刘洋， 周健， 付建新. 饱和砂土流固耦合细观数值模型及其在液化分析中的应用[J]. 水利学报， 2009， 40（2）： 250256.
[12]JACKSOM R. Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[J]. Chemical Engineering Science， 1997， 52（15）： 24572469.
[13]ERGUN S. Fluid flow through packed columns[J]. Chemical Engineering Progress， 1952， 43（2）： 8994.
[14]杜延龄， 韩连兵. 土工离心模型试验技术[M]. 北京： 中国水利水电出版社， 2010.
[15]苏栋， 李相崧. 砂土自由场地地震响应的离心机试验研究[J]. 地震工程与工程振动， 2006， 26（2）： 166170.
[16]陈国兴， 金丹丹， 常向东， 等. 最近20年地震中场地液化现象的回顾与土体液化可能性的评价准则[J]. 岩土力学， 2013， 34（10）： 27372755.
[17]梁孟根， 梁甜， 陈云敏. 自由场地液化响应特性的离心机振动台试验[J]. 浙江大学学报（工学版）， 2013， 47（10）： 18051814.
[18]汪闻韶. 土体液化与极限平衡和破坏的区别和联系[J]. 岩土工程学报， 2005， 27（1）： 110.
（编辑赵勉）
摘要：为研究饱和砂土自由场地的地震液化过程，基于颗粒离散元与计算流体动力学耦合方法，在PFC3D程序中通过二次开发实现对自由场地基底不规则地震波的输入，研究在实测Kobe地震波激励下该液化过程的宏细观力学响应，并与已有的福建中细砂离心机模型试验结果进行对比。数值模拟结果表明：试样发生液化的时刻对应于基底地震加速度峰值时刻，浅部土层先于深部土层发生液化，超静孔隙水压力的消散是从深部土层逐渐向浅部土层发展的；伴随液化过程，各深度土层的有效应力路径均从初始应力状态向应力空间的原点移动；数值试验呈现出的上述规律均与实际砂土离心机模型试验呈现出的规律保持一致。数值试样宏观的液化特征在微观机理上可由流体对颗粒的归一化拖曳力、局部孔隙比、平均配位数、颗粒位移的变化规律来相应表征。
关键词：
饱和砂土； 自由场地； 地震荷载； 液化； 离散元模拟
中图分类号： TU441.5
文献标志码： A
Discrete element simulation on seismic liquefaction
of saturated freefield sand deposits
SHI Danda1， WANG Fei2
（1. College of Ocean Science and Engineering， Shanghai Maritime University， Shanghai 201306， China；
2. East China Architectural Design & Research Institute， Shanghai 200070， China）
Abstract：
In order to research the seismic liquefaction of saturated freefield sand deposits， using the coupled method of the particle discrete element method and computational fluid dynamics， both macro and micro mechanical responses of the liquefaction are explored subjected to the real Kobe earthquake wave excitation， where the input of irregular seismic waves on freefield base is realized by selfdeveloped procedures in PFC3D. The numerical results are compared with the centrifugal model test results on Fujian sand， which is a uniformly graded medium to fine sand. It is found that the exact instant for sample liquefaction corresponds to the peak time of the base seismic acceleration. Liquefaction happens in shallow sand layers prior to deep sand layers， whilst the excess pore water pressure dissipates earlier in deep layers than in shallow layers. During liquefaction， the effective stress path of each sand layer moves from the initial stress condition to the coordinate origin in the stress space. The above laws obtained by numerical tests are consistent with the laws obtained by the real sand centrifugal model test. The macroscopic liquefaction behaviors are correspondingly characterized by the normalized fluid drag force on particles， the porosity， the average coordination number and the particle displacement in the micro mechanism.
Key words：
saturated sand； free field； earthquake loading； liquefaction； discrete element simulation
0引言
在地震災害中，饱和砂土地基的液化是造成大量建筑物震害的重要原因。关于砂土液化宏观机理的研究，自SEED等[1]于1966年提出“初始液化”概念以来已近半个世纪，积累的成果较多[23]。然而，关于砂土液化细观机理的研究，由于试验条件和研究手段等方面的制约，目前成果相对较少。离散元方法（Discrete Element Method， DEM）由于制样便捷、细观组构信息丰富且容易获取等优点，在粒状土静、动力特性及其机理研究方面发挥了重要作用。
目前采用DEM模拟砂土液化通常有两种方法。一种称为常体积法，SITHARAM[4]、GONG等[5]、史旦达等[6]、SUZUKI等[7]均采用该方法模拟了砂土的振动液化特性，该方法的优点是计算效率高，缺陷是试样本身未设置流体单元，不能考虑固液二相细观相互作用。另一种为流固耦合方法，即颗粒离散元与计算流体动力学（Computational Fluid Dynamics， CFD）耦合，HAKUNO等[8]最早提出可用于饱和砂土液化分析的流固耦合方法，但其提出的在孔隙尺度上的流体效应模拟计算效率很低。为提高计算效率，NAKASA等[9]对HAKUNO等提出的方法进行了优化，流体采用比孔隙尺度稍大的流体网格来模拟，每个流体网格内能够容纳若干个颗粒，流体流速与压力梯度之间满足Darcy定律。ZEGHAL等[10]在前人研究基础上提出了一种更为实用的耦合方法，流体运动满足均一化NavierStokes方程，流体对颗粒的作用力采用已有的半经验公式表示，目前应用较广的离散元程序PFC3D提供的流固耦合分析模块就采用了这种方法。刘洋等[11]基于PFC3D流固耦合模块，研究了饱和砂土的液化响应，但局限于规则的简谐荷载输入，并未考虑实际地震荷载，数值模拟结果也未与物理试验结果展开比较。
本文基于PFC3D流固耦合模块，通过二次开发实现基底不規则地震波的输入，研究实测Kobe地震波激励下饱和砂土自由场地的液化宏细观力学响应，并与典型的离心机试验结果进行对比。除分析超静孔隙水压力、有效应力路径等宏观响应外，重点探讨流体对颗粒的拖曳力、孔隙比、平均配位数和颗粒位移等细观力学响应，用来揭示液化过程中宏观液化规律与细观液化机理之间的宏细观关联。
1流固耦合原理简介
1.1液相流体
假定液相流体不可压缩且密度不变，其运动满足均一化NavierStokes连续方程和动量方程[12]，
nt+·（nvf）=0 （1）
ρf（nvf）t+·（nvfvf）=·（nτf）-
fi+nρfg（2）
式中：n为孔隙率；ρf为流体质量密度；vf和τf分别为流体平均速度矢量和平均应力张量；fi为单位体积内平均流
体颗粒拖曳力矢量；g为重力加速度；
为哈密顿算子。
1.2流体颗粒拖曳力
流体颗粒相互作用拖曳力可采用Ergun半经验公式[13]计算，
fi=150μf（1-n）2n·dp2（vf-vp）+
1.75（1-n）ρfvf-vpdp（vf-vp） （3）
式中：μf为流体黏滞系数；dp为流体单元内的颗粒平均粒径；vp为颗粒平均速度矢量。
1.3固相颗粒
固相颗粒采用离散元模拟，其运动满足牛顿第二定律，且考虑孔隙流体的作用，单个颗粒的平动和转动方程分别为
mpv·p=mp·g+cfc+fd
（4）
Ipω·p=crc×fc （5）
式中：mp和Ip分别为颗粒的质量和惯性矩；vp和ωp分别为颗粒的平动和转动速度矢量；fc为接触点处粒间作用力矢量；rc为颗粒中心至接触点的方向矢量；fd为流体施加于颗粒的作用力，包括浮力项和流体颗粒拖曳力项。
fd=-pf+fi1-nVp（6）
式中：Vp为颗粒体积；pf为孔隙流体平均压力梯度。
2试样制备及地震荷载输入
2.1数值试样制备
在PFC3D中可以通过对数值试样施加一定倍数的重力加速度来还原现场原型的场地条件，数值模型尺寸为140 mm（长）×60 mm（宽）×250 mm（高），对数值试样施加30倍重力加速度，还原至现场原型的场地尺寸为4.2 m（长）×1.8 m（宽）×7.5 m（高）。本文涉及的土层深度和分析得到的物理力学量值均为换算至场地原型的数值，换算方法满足离心试验相似性原理[14]。试样制备采用重力沉积法，颗粒为三维圆球颗粒，颗粒粒径范围5～10 mm，颗粒平均粒径为7.65 mm，初始孔隙比为0.727，参考福建标准砂的最大孔隙比emax=0.848和最小孔隙比emin=0.519，计算得到数值试样的相对密度约为37%。生成的颗粒数量为5 198个，数值试样如图1a）所示，试样顶部为自由面，底部为刚性墙，侧面采用周期边界，地震荷载通过底部（基底）墙体输入。
试样内布置固定位置的流体网格，如图1b）所示。流体网格在X，Y，Z三个方向上的尺寸分别为28 mm，20 mm和21 mm，每个流体网格内约有颗粒32个。试样内部设置4个量测圈，用于监测孔隙比、平均配位数等细观组构量变化。量测圈的设置见图1a），量测圈布置在试样竖向中心轴位置，第1，2，3，4号量测圈的中心位置距试样顶面的高度分别为0.13h，0.29h，0.54h和0.79h，h为试样高度。第1，2，3，4号量测圈分别对应于第1，2，3，4层土层，且分别对应于图1b）中Z方向第11，9，6，3号流体网格。
数值试样中颗粒和流体的细观参数设置见表1。颗粒之间的接触采用线性接触模型，流体黏滞系数为5.02 Pa·s。DEM采用中心差分格式，计算得到的平均时步为2.5×10-6 s。为保证计算效率，每执行100步离散元计算后执行1步CFD计算，即CFD计算时步为DEM计算时步的100倍，为2.5×10-4 s。
2.2地震荷载输入
地震荷载通过基底墙体输入，采用实测Kobe地震波（见图2a）），加速度峰值为0.32g，持续时间为50 s。数值模拟中在基底墙体不能直接输入加速度，需要通过对加速度时程曲线积分获得速度时程曲线（见图2b）），再通过给基底墙体赋予速度来实现。对模型施加的加速度为30g，根据相似性原理，模型基底激振持续时间的相似比为1/30；原型激振
持续时间为50 s，换算至模型激振持续时间即为1.67 s。在给基底墙体输入速度的过程中，在速度时程曲线输入文件中每隔0.000 4 s读取一个数据点（选择较小的时间间隔能确保输入速度曲线的精度），通过PFC3D内嵌的Fish语言编程将速度值赋给基底墙体，完成基底地震荷载输入。
3超静孔隙水压力及归一化拖曳力分析
3.1超静孔隙水压力
超静孔隙水压力比ru，即为激振引起的试样超静孔隙水压力u与初始竖向有效应力之比，其数值大小及随时间变化的规律已被广泛用于试样液化分析中。本文数值模拟中，激振引发的u可以通过监测流体网格中心节点处的流体压力得到。4个深度土层的监测结果见表2，4个深度土层的ru随时间t的变化曲线见图3a）。
a）数值模拟结果
b）离心机模型试验结果[17]
图3数值模拟与离心机模型试验超静孔压比曲线对比
分析图3a）并结合表2可知：（1）整个试样发生液化的时刻约为10 s时，液化发生的时刻与输入的地震加速度峰值时刻（见图2a））保持一致；浅部土层ru峰值最大，随着深度的增加ru值略有减少，但数值上仍超过0.9，可以认为整个试样均发生了液化；从液化发生的时刻看，浅部土层先发生液化，深部土层后发生液化；已有的室内试验和液化现场实例调查均表明[1516]，对于饱和砂层的液化分布，随着埋深的增加，土层所受的初始竖向有效应力逐渐增加，因此深部土层较难达到初始液化。（2）第10～18 s，各深度土层的ru均出现迅速下降，u消散较快。（3）大概25 s后，试样的u消散过程趋于稳定，残余ru随土层深度的增加而减小，第1层土层的残余ru约为0.48，到第4层土层时，其残余ru减小为0.16左右。
梁孟根等[17]进行了饱和砂土自由场地液化响应的离心机振动台试验，试验砂料为福建中细砂，模型箱内的试样尺寸为730 mm（长）×330 mm（宽）×400 mm（高），加速度为50g，基底输入德阳波250 s，加速度峰值为0.392g，加速度峰值出现的时刻为
20 s。试样内沿高度方向自上而下共布置5个孔压计（P1～P5），为了与数值试样4个土层深度位置尽量保持对应，选取P1，P2，P3和P5的量测结果进行对比分析（P1，P2，P3和P5孔壓计距试样顶面的深度依次为0.16h，0.33h，0.50h和0.85h，h为试样高度）。该离心机模型试验输入的为德阳波，本文数值模拟采用Kobe地震波，在加速度峰值等特性上存在一定差异，因此仅从定性规律角度对比分析离心机模型试验和数值模拟的u变化规律。
图4为福建中细砂与本文数值试样的级配对比曲线。分析图4可得：福建中细砂平均粒径d50=0.16 mm，不均匀因数Cu=1.6，曲率因数Cc=0.96；数值试样d50=8.1 mm，Cu=1.422和Cc=1.051。与福建中细砂相比，数值试样除
适当放大颗粒平均粒径外（减少一定的颗粒数量，保证数值模拟可以实施），Cu和Cc值均
与福建中细砂的相近，可以认为数值试样的级配分布基本与福建中细砂的一致。
图3b）给出了文献[17]离心机模型试验实测的ru随时间变化的曲线。对比分析图3a）与图3b）可得：（1）从各深度ru曲线的整体发展形态看，数值试样和离心试样均经历了u增长、试样发生液化和液化后u消散的过程。（2）数值试样和离心试样发生液化的时刻均对应于输入加速度峰值的时刻（离心机模型试验加速度峰值时刻为20 s），在ru逐渐增加并接近1.0的过程中，u均呈现振动式发展，试样发生了循环液化。（3）数值试样与试验试样均呈现深部土层u消散快，浅部土层u消散慢的规律，这在物理机理上与底部边界不排水而顶部自由面排水有关，鉴于这一特点，液化后土层中的u消散路径为自下而上，消散过程中深部土层的u必然会传递至浅层并造成浅层u的累积，因此试样会呈现出深部土层u消散快，浅部土层u消散慢的规律。（4）与离心试样相比，数值试样u消散过程更为明显，这是由于离心机模型试验输入的德阳波在加速度峰值时刻后仍有较高水平的加速度次峰值输入，导致u维持在较高的水平。
进一步对比分析数值模拟与离心机模型试验结果，将数值试样和试验试样在液化前（即u增长阶段）不同时刻u随深度的变化曲线分别绘于图5a）和图5b），在u振荡上升段，取一个振动周期内u的平均值。对比分析图5a）与图5b）可知，虽然数值试样与试验试样发生液化的具体时间不同，但数值试样表现出与试验试样一致的u变化规律：（1）在液化前，各深度u均随时间的增加而增大，直至达到初始竖向有效应力水平（即发生液化）；（2）浅部土层u增长较快，随着深度的增加u增长逐渐变慢，说明浅部土层先发生液化，深部土层后发生液化。
3.2归一化拖曳力
激振引起的试样宏观u在细观机理上对应于流体对颗粒的拖曳力大小。定义归一化拖曳力fs为颗粒受到的流体竖向拖曳力与该颗粒有效自重的比
值，当fs≥1，在细观上即表现为试样发生液化。图6给出了数值试样4个深度土层fs随时间的变化规律，fs<0表示颗粒受到的合力方向发生了改变，颗粒开始产生向下的加速度。
分析图6，并结合3.1节u的规律分析可知：（1）fs的变化规律对应于u的变化规律，fs出现峰值的时刻与试样的u出现峰值的时刻基本一致，不同深度土层的fs峰值均接近或略大于1，表明试样整体已发生液化。（2）在整个激振过程中，伴随u的增长和消散过程，fs经历了先增大、后减小，并逐渐趋于稳定的过程。（3）由第3.1节分析可知，u的消散是从深部土层向浅部土层发展的，因此试样液化后，深部土层fs越早进入残余稳定状态，其数值上越接近于0，而浅部土层由于受到深部土层向上排水产生的残余孔隙水压力累积效应的影响，在u消散过程中还经历了fs缓慢增大的过程，激振结束时（50 s）仍有一定数值的残余fs存在。
4其他监测结果分析
4.1有效应力路径
利用布置在试样内部的量测圈可以测得试样内部6个有效应力分量
σ′xx，σ′yy，σ′zz，σ′xy，σ′yz和σ′zx，进一步可计算得到平均有效主应力p′和广义剪应力q。
p′=13（σ′xx+σ′yy+σ′zz） （7）
q=12
（（σ′xx-σ′yy）2+（σ′yy-σ′zz）2+
（σ′zz-σ′xx）2+6（σ′xy+σ′yz+σ′zx））12（8）
图7绘出了数值试样不同深度土层在p′～q应力空间内的有效应力路径。分析图7可知：从激振开始到数值试样液化前，第1～4层土体的有效应力路径均从初始应力状态向p′～q应力空间的原点移动，应力路径自上而下行进，p′和q值均逐渐减小，并趋近于零；当试样达到液化状态后，随着u的循环振荡，试样发生循环液化，应力路径在原点附近来回折返，这一现象在浅部土层中表现得尤为明显。值得注意的是，无论应力路径如何折返，4个深度土层的应力路径均没有逾越试样的应力路径包线（即试样强度包线）。汪闻韶[18]研究了砂土液化与极限平衡和破坏之间的关联，指出砂土液化过程会出现循环活动性，应力路径会发生来回折返现象，但应力路径始终不会逾越砂样的极限平衡强度包线。因此，本文数值试样呈现的有效应力路径变化规律符合实际砂样的液化特性。
4.2颗粒位移与表面沉降
为分析液化过程中内部颗粒的运动情况，在距试样顶面1.1 m深度处选取标志颗粒进行追踪。该颗粒初始坐标为X=2.02 m，Y=0.84 m，Z=6.40 m，位于第1层土体内。图8a）绘出了整个激振过程中该标志颗粒的三维空间运动轨迹，图中实心圆点为颗粒运动的起始位置，空心圆点为颗粒运动的最终位置。分析图8a）可知：在激振过程中该标志颗粒的X方向坐标持续增加，说明试样发生了沿激振方向的累积水平位移；
在激振过程中该标志颗粒的
Y方向坐标的变化随机性较强，没有表现出一定的规律；
在激振过程中该标志颗粒的
Z方向坐标值经历了先增大后减小的过程，反映了液化过程中颗粒上浮、液化后随着u消散颗粒逐渐沉积的运动过程。
图8b）为整个激振过程中试样平均高度随时间的变化规律。由图8b）可知：在試样发生液化的时刻（10 s左右），试样平均高度达到最大，由于流体对颗粒拖曳力的作用，试样内的颗粒发生了较为明显的上浮，与初始高度相比整个试样高度增加了约0.17 m；之后，随着u的消散，试样的平均高度迅速减小，即试样发生了表面沉降，与初始时刻相比，最终的表面沉降量约为0.33 m，这一规律反映了饱和砂土自由场地的液化后震陷特征。
a）颗粒运动轨迹
b）试样平均高度变化
4.3孔隙比
孔隙比变化规律能够反映试样液化过程中的体积变化特征，利用设置在试样内部的量测圈（见图1a））可以测得4个深度土层孔隙比随时间的变化规律，见图9。分析图9可知：（1）各层土体的初始孔隙比略有差别，由于制样时的重力沉积效应，深部土层的初始孔隙比略小于浅部土层的。（2）在液化发生前，各层土体的孔隙比均有不同程度的增加，其中第3层和第4层土体孔隙比的增加较为明显，其原因与流体压力从基底输入导致深部颗粒的上浮较为明显有关。（3）当试样处于液化后u消散阶段时，各层土体的孔隙比均逐渐减小，与初始孔隙比相比，各层土体在液化后的最终孔隙比均小于初始孔隙比，表明试样在振动液化后变得更为密实。
4.4平均配位数
平均配位数的变化能够反映试样在液化过程中内部结构稳定性的演化，图10给出了4个深度土层平均配位数随时间的变化。
由图10可知：（1）试样发生液化时，各层土体的平均配位数均从初始时刻的4左右迅速降低至1左右，在这一过程中，浅部第1层土体平均配位数的降低最快，深部土层的平均配位数降低稍慢，其原因与浅部土层先液化、深部土层后液化的规律有关。（2）在u消散阶段，各层土体的平均配位数均得到了一定程度的恢复，深部第3层和第4层土体的平均配位数在激振开始30 s后趋于稳定，而浅部第1层和第2层土体的平均配位数在整个u消散过程中都在持续增加，因为深部土层u消散快、浅部土层u消散慢，所以此处反映的不同深度平均配位数的变化规律与宏观u的消散规律是对应的。（3）液化后，各深度土层的最终配位数均略大于初始配位数，这一规律也反映了试样振动液化后变密的特征。
5结论
基于PFC3D颗粒离散元与计算流体动力学耦合模拟方法，研究了真实地震荷载激励下饱和砂土自由场地液化过程的宏细观力学响应，并与已有的离心机模型试验结果进行了对比，得到的主要结论有：
（1）整个试样发生液化的时刻对应于基底地震加速度峰值时刻，因为深部土层的初始竖向有效应力高于浅部土层的，所以浅部土层先于深部土层发生液化，浅部土层的峰值超静孔隙水压力比ru也更大；液化后超静孔隙水压力u的消散先从深部土层开始，逐渐向浅部土层发展，浅部土层的残余ru大于深部土层的；数值试样呈现出的上述u变化规律均与福建中细砂自由场地地震液化的离心机模型试验结果一致。
（2）在试样液化过程中，各深度土层的有效应力路径均从初始应力状态向p′～q（p′为平均有效主应力，q为广义剪应力）应力空间的原点移动，p′和q值均逐渐减小并趋于零；当试样达到液化状态后，随着u的循环振荡，应力路径在原点附近来回折返，但应力路径不会逾越试样的极限平衡强度包线，这一现象符合实际砂样振动液化过程中的有效应力路径特征。
（3）对应于u的消散是从深部土层向浅部土层发展的宏观液化规律，在试样液化后，深部土层的残余归一化拖曳力更接近于0，且平均配位数也更早进入残余稳定状态，而浅部土层中仍存在一定数值的残余拖曳力作用，且平均配位数随着u的消散逐渐增加；对应于自由场地的宏观液化震陷特征，激振结束后各深度土层的最终孔隙比均小于初始孔隙比，最终配位数均略大于初始配位数，且试样整体发生了约0.33 m的表面沉降。
参考文献：
[1]SEED H B， LEE K L. Liquefaction of saturated sands during cyclic loading[J]. Journal of the Soil Mechanics and Foundations Division， 1966， 92（6）： 105134.
[2]莊海洋， 胡中华， 王瑞， 等. 饱和南京细砂初始液化后特大流动变形特性试验研究[J]. 岩土工程学报， 2016， 38（12）： 21642174.
[3]LU Xilin， HUANG Maosong. Static liquefaction of sands under isotropically and K0consolidated undrained triaxial conditions[J]. Journal of Geotechnical and Geoenvironmental Engineering， 2015， 141（1）： 04014087.
[4]SITHARAM T G. Discrete element modelling of cyclic behavior of granular materials[J]. Geotechnical and Geological Engineering， 2003， 21： 297329.
[5]GONG Guobin， THORNTON C， CHAN A H C. DEM simulations of undrained triaxial behavior of granular material[J]. Journal of Engineering Mechanics， 2012， 138（6）： 560566.
[6]史旦达， 周健， 刘文白， 等. 初始组构影响砂土液化势的细观数值模拟研究[J]. 水利学报， 2011， 42（7）： 766774.
[7]SUZUKI K， KUHN M R. Discrete element simulations of cyclic biaxial shear of a granular material with oval particles[J]. International Journal of Geomechanics， 2014， 14（3）： 06014005.
[8]HAKUNO M， TARUMI Y. A granular assembly simulation for the seismic liquefaction of sand[J]. Proceedings of Japan Society of Civil Engineers， 1988， 5（2）： 129138.
[9]NAKASA H， TAKEDA T， ODA M. A simulation study on liquefaction using DEM[C]//Earthquake Geotechnical Engineering. Rotterdam： Balkema， 1999， 637642.
[10]ZEGHAL M， El SHAMY U. A continuumdiscrete hydromechanical analysis of granular deposit liquefaction[J]. International Journal of Numerical and Analytical Methods in Geomechanics， 2004， 28： 13611383.
[11]刘洋， 周健， 付建新. 饱和砂土流固耦合细观数值模型及其在液化分析中的应用[J]. 水利学报， 2009， 40（2）： 250256.
[12]JACKSOM R. Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[J]. Chemical Engineering Science， 1997， 52（15）： 24572469.
[13]ERGUN S. Fluid flow through packed columns[J]. Chemical Engineering Progress， 1952， 43（2）： 8994.
[14]杜延龄， 韩连兵. 土工离心模型试验技术[M]. 北京： 中国水利水电出版社， 2010.
[15]苏栋， 李相崧. 砂土自由场地地震响应的离心机试验研究[J]. 地震工程与工程振动， 2006， 26（2）： 166170.
[16]陈国兴， 金丹丹， 常向东， 等. 最近20年地震中场地液化现象的回顾与土体液化可能性的评价准则[J]. 岩土力学， 2013， 34（10）： 27372755.
[17]梁孟根， 梁甜， 陈云敏. 自由场地液化响应特性的离心机振动台试验[J]. 浙江大学学报（工学版）， 2013， 47（10）： 18051814.
[18]汪闻韶. 土体液化与极限平衡和破坏的区别和联系[J]. 岩土工程学报， 2005， 27（1）： 110.
（编辑赵勉）
