D0I:10.13374/1.issnl00103.2008.05.022 第30卷第5期 北京科技大学学报 Vol.30 No.5 2008年5月 Journal of University of Science and Technology Beijing My2008 砂沸现象的细观数值模拟及影响因素分析 刘洋)吴顺川)周健) 1)北京科技大学土木与环境工程学院，北京1000832)同济大学地下建筑与工程系，上海200092 摘要采用饱和砂土流固耦合的细观力学模型，通过颗粒流数值模拟了在临界水力梯度条件下饱和砂土的细观响应·讨论 了在建议模型框架内临界水力梯度、颗粒刚度等因素对模拟结果的影响，并从细观机理上做了进一步分析·模拟结果表明，在 临界流或超临界流条件下，试样最终的稳态响应分成两部分：一部分颗粒停留在试样顶部，其颗粒数量正比于施加的水力梯 度值；另一部分颗粒则在试样底部形成某种形式的对流运动，细观数值模型分析所得的临界水力梯度比连续介质力学分析得 到的值要小，这个差别反映了模型的动力特性·在临界流或超临界流条件下，砂沸时水的流动除局部区域外基本上仍属于层 流范围，但已偏离了达西定律. 关键词砂沸：颗粒流：细观模型：砂土液化 分类号TU441.+4 Micro-numerical simulation of sand boil and analysis of its influencing factors LIU Yang),WU Shunchuan).ZHOU Jian2) 1)School of Civil and Environmental Engineering.University of Science and Technology Beijing Beijing 100083.China 2)Department of Geotechnical Engineering.Tongji University,Shanghai 200092.China ABSTRACI Particle flow code(PFC)simulation was performed to investigate meso-response of sand soil under critical flow condi- tion using a coupled liquid particle model.The effects of hydraulic gradient amplitude and particle stiffness were discussed and the oc- currence of sand boil was analyzed from the viewpoint of micro"mechanism.The numerical results indicated that under critical or over- critical gradients.a layer of clustered particles formed near the top wall.and the remainder of the particles formed a suspension and underwent convective-like motion.The thickness of the top layer was proportional to the gradient amplitude.The gradient amplitude based on microanalysis was a little lower than that based on the traditional theory,which reflected the dynamic behaviour of the adopted model.The pore flow was mainly laminar but deviated from Darcy 's law under critical or over-critical gradients. KEY WORDS sand boil:particle flow code (PFC):micro-model:sand liquefaction 饱和砂土的孔隙水压力变化是土体对外荷载作 分布与地面位移的关系；Holzer研究了非地震荷 用的宏观响应，当饱和土体中的超孔隙水压力增长 载引起的砂沸现象。·关于砂沸细观机理方面的研究 到一定程度时引起液化，土体由固态转变成悬浮态· 目前还比较少 砂沸的产生是砂土液化中一种特殊的现象，当超孔 本文采用饱和砂土固液耦合细观力学模型，通 隙水压力逐渐增大，达到土体本身的重力时，在土体 过颗粒流P℉C数值实验观察在临界水力梯度条件 某些薄弱的部位就产生了砂沸现象，这时的孔隙水 下饱和砂土的细观响应，讨论了在建议模型框架内 压或水力梯度称为临界孔隙水压力或临界水力梯 临界水力梯度、颗粒刚度等因素对模拟结果的影响， 度，汪闻韶把砂沸作为一种液化机理，讨论了砂沸 并从细观机理上对砂沸现象做了进一步的分析， 过程中应力状态的演化-2)；Bardet研究了San Francisco地区Loma Prieta地震时的砂沸现象及其 1饱和砂土固液耦合细观数值模型 收稿日期：2007-03-22修回日期.2007-04-26 1.1砂沸的连续介质理论 基金项目：教育部新世纪优秀人才支持计划资助项目(No·NCET一 根据传统的土力学理论，渗透力是造成砂沸现 06-0084) 作者简介：刘洋(1979一)男，讲师，博士， 象的原因，在渗流场的作用下，当渗透力大于上覆 E-mail:imaginationly@163.com 土柱重量时，土体将被抬起，粘性土地基中将产生流
砂沸现象的细观数值模拟及影响因素分析 刘 洋1） 吴顺川1） 周 健2） 1） 北京科技大学土木与环境工程学院北京100083 2） 同济大学地下建筑与工程系上海200092 摘 要 采用饱和砂土流固耦合的细观力学模型通过颗粒流数值模拟了在临界水力梯度条件下饱和砂土的细观响应．讨论 了在建议模型框架内临界水力梯度、颗粒刚度等因素对模拟结果的影响并从细观机理上做了进一步分析．模拟结果表明在 临界流或超临界流条件下试样最终的稳态响应分成两部分：一部分颗粒停留在试样顶部其颗粒数量正比于施加的水力梯 度值；另一部分颗粒则在试样底部形成某种形式的对流运动．细观数值模型分析所得的临界水力梯度比连续介质力学分析得 到的值要小这个差别反映了模型的动力特性．在临界流或超临界流条件下砂沸时水的流动除局部区域外基本上仍属于层 流范围但已偏离了达西定律． 关键词 砂沸；颗粒流；细观模型；砂土液化 分类号 TU441∙＋4 Micro-numerical simulation of sand boil and analysis of its influencing factors LIU Y ang 1）W U Shunchuan 1）ZHOU Jian 2） 1） School of Civil and Environmental EngineeringUniversity of Science and Technology BeijingBeijing100083China 2） Department of Geotechnical EngineeringTongji UniversityShanghai200092China ABSTRACT Particle flow code （PFC） simulation was performed to investigate meso-response of sand soil under critical flow condition using a coupled liquid-particle model．T he effects of hydraulic gradient amplitude and particle stiffness were discussed and the occurrence of sand boil was analyzed from the viewpoint of micro-mechanism．T he numerical results indicated that under critical or overcritical gradientsa layer of clustered particles formed near the top walland the remainder of the particles formed a suspension and underwent convective-like motion．T he thickness of the top layer was proportional to the gradient amplitude．T he gradient amplitude based on microanalysis was a little lower than that based on the traditional theorywhich reflected the dynamic behaviour of the adopted model．T he pore flow was mainly laminar but deviated from Darcy’s law under critical or over-critical gradients． KEY WORDS sand boil；particle flow code （PFC）；micro-model；sand liquefaction 收稿日期：2007-03-22 修回日期：2007-04-26 基金项目：教育部新世纪优秀人才支持计划资助项目（No．NCET－ 06－0084） 作者简介：刘 洋（1979－）男讲师博士 E-mail：imaginationly＠163．com 饱和砂土的孔隙水压力变化是土体对外荷载作 用的宏观响应．当饱和土体中的超孔隙水压力增长 到一定程度时引起液化土体由固态转变成悬浮态． 砂沸的产生是砂土液化中一种特殊的现象．当超孔 隙水压力逐渐增大达到土体本身的重力时在土体 某些薄弱的部位就产生了砂沸现象这时的孔隙水 压或水力梯度称为临界孔隙水压力或临界水力梯 度．汪闻韶把砂沸作为一种液化机理讨论了砂沸 过程中应力状态的演化［1－2］；Bardet ［3］ 研究了 San Francisco 地区 Loma Prieta 地震时的砂沸现象及其 分布与地面位移的关系；Holzer ［4］ 研究了非地震荷 载引起的砂沸现象．关于砂沸细观机理方面的研究 目前还比较少． 本文采用饱和砂土固液耦合细观力学模型通 过颗粒流 PFC 数值实验观察在临界水力梯度条件 下饱和砂土的细观响应．讨论了在建议模型框架内 临界水力梯度、颗粒刚度等因素对模拟结果的影响 并从细观机理上对砂沸现象做了进一步的分析． 1 饱和砂土固液耦合细观数值模型 1∙1 砂沸的连续介质理论 根据传统的土力学理论渗透力是造成砂沸现 象的原因．在渗流场的作用下当渗透力大于上覆 土柱重量时土体将被抬起粘性土地基中将产生流 第30卷 第5期 2008年 5月 北 京 科 技 大 学 学 报 Journal of University of Science and Technology Beijing Vol．30No．5 May2008 DOI:10．13374／j．issn1001－053x．2008．05．022
.474 北京科技大学学报 第30卷 土，砂土地基则产生砂沸.作用在土体上的渗透力 学(CFD)程序与离散元法的颗粒流(PFC)模拟环境 和重力相等时的水力梯度称为临界水力梯度，记为 相结合来求解平均的流体运动方程，孔隙流体的压 i,连续介质力学理论分析结果为： 力和速度由平均的Navier Stokes方程求解，求解采 i=(Y-Y)/Y=Y/Y (1) 用有限体积积分格式.基于SIMPLE结构和向上差 式(1)也可以改写为： 分格式的CFD计算技术用来求解离散的平均 i=(d.-1)/(1+e) (2) Navier Stokes方程，计算采用空间交错有限体积结 式中，Y、Y为土的重度和浮重度，Y为水的重度， 构以保证计算的稳定性，关于饱和砂土流固耦合细 d。为土粒重度，e为孔隙比 观力学模型更详细的计算步骤可参见有关文献[9] 砂沸的产生主要是渗透压产生的渗透不稳定性 采用这种方法可以分析流体速度、流线以及压 引起的，渗透压力与土中的孔隙水压力分布有关. 力分布，同时可以得到颗粒的运动和应力分布，而这 这种孔隙水压力的变化可以是动力的（如地震作 些细观参量在实验室是很难或几乎不可能观察到 用)，也可能为非动力的原因（如地下水位的上升）， 的 因此式(1)和式(2)并未涉及砂土的应变特性（剪胀 2数值实验程序 或剪缩) 1.2饱和砂土流固耦合的细观力学模型 2.1数值试样生成 与传统的土力学理论不同，本文采用饱和砂土 试样采用PFC随机生成器生成，模拟容器大小 流固耦合的细观力学模型来模拟砂沸现象，Tsu] 为15mm×27mm×90mm,然后在重力作用下沉积 首先将离散元(DEM)模拟和均一化流体公式相结 达到初始平衡，这个过程一般需要较多的时步，模 合用于分析流化床中泡沫的形成，此后，这种数值 拟容器侧壁采用不透水的刚性墙边界条件，上下为 模型被广泛地用于一些化学过程的模拟可.本文 由P℉C线墙模拟的透水边界条件，图1是生成的模 采用Tsui提出的方法建立饱和砂土流固耦合的细 型图，表1是数值实验采用的指标参数.本构模型 观力学模型 采用Herthz接触模型，为了加速运算，颗粒的剪切 模型采用离散单元法的颗粒流理论(P℉C)模拟 模量采用了比实际土颗粒（如石英）小的值 饱和砂土中砂颗粒的力学行为，采用均一化流体 0 CFD计算技术模拟孔隙流体的运动[8)].对于固相砂 9 土，通过求解运动和动量方程模拟颗粒运动；对于液 相水，通过求解平均的Navier-Stokes方程模拟孔隙 流体的运动，这两组方程通过施加于颗粒表面非滑 动的边界条件相耦合· 流固耦合方程的解采用显式结构，在计算的第 I时步，在接触力、外部荷载以及流体拽曳力作用下 颗粒的位移采用DEM技术计算，采用显式的中心 差分来积分固相的运动方程，计算稳定所需的时步 图1砂沸现象模拟P℉C数值模型试样 由粒间接触的刚度和阻尼系数决定，计算流体动力 Fig-I Sand samples by PFC for sand boil simulation 表1颗粒流模型细观参数表 Table 1 Parameters for numerical simulation 流体 计算时步/s 粒径/摩擦因数/孔隙 剪切模 有限体积 试样尺寸 密度/ 粘滞系数/ mm 泊松比 比 量/MPa 网格 PFC CFD (kg'm3) (Pas) 15mm×27mm×90mm0.3-0.70.5/0.30.42 2.9 1000 10-3 3mmX3mm×3mm2.0X10-62.0X10-3 生成的颗粒试样在重力和流体共同作用下砂土 个容器高度的67%左右 颗粒沉积平衡，循环计算消除不平衡力，并使此时试 2.2数值模拟过程 样的超孔压为零.初始模拟的沉积砂土颗粒约占整 首先在模拟容器中充满水，接着施加一个向上
土砂土地基则产生砂沸．作用在土体上的渗透力 和重力相等时的水力梯度称为临界水力梯度记为 icr连续介质力学理论分析结果为： icr＝（γ－γw）／γw＝γ′／γw （1） 式（1）也可以改写为： icr＝（ ds－1）／（1＋e） （2） 式中γ、γ′为土的重度和浮重度γw 为水的重度 ds 为土粒重度e 为孔隙比． 砂沸的产生主要是渗透压产生的渗透不稳定性 引起的渗透压力与土中的孔隙水压力分布有关． 这种孔隙水压力的变化可以是动力的（如地震作 用）也可能为非动力的原因（如地下水位的上升） 因此式（1）和式（2）并未涉及砂土的应变特性（剪胀 或剪缩）． 1∙2 饱和砂土流固耦合的细观力学模型 与传统的土力学理论不同本文采用饱和砂土 流固耦合的细观力学模型来模拟砂沸现象．Tsuji ［5］ 首先将离散元（DEM）模拟和均一化流体公式相结 合用于分析流化床中泡沫的形成．此后这种数值 模型被广泛地用于一些化学过程的模拟［6－7］．本文 采用 Tsuji 提出的方法建立饱和砂土流固耦合的细 观力学模型． 模型采用离散单元法的颗粒流理论（PFC）模拟 饱和砂土中砂颗粒的力学行为采用均一化流体 CFD 计算技术模拟孔隙流体的运动［8］．对于固相砂 土通过求解运动和动量方程模拟颗粒运动；对于液 相水通过求解平均的 Navier－Stokes 方程模拟孔隙 流体的运动．这两组方程通过施加于颗粒表面非滑 动的边界条件相耦合． 流固耦合方程的解采用显式结构．在计算的第 I 时步在接触力、外部荷载以及流体拽曳力作用下 颗粒的位移采用 DEM 技术计算．采用显式的中心 差分来积分固相的运动方程计算稳定所需的时步 由粒间接触的刚度和阻尼系数决定．计算流体动力 学（CFD）程序与离散元法的颗粒流（PFC）模拟环境 相结合来求解平均的流体运动方程．孔隙流体的压 力和速度由平均的 Navier－Stokes 方程求解求解采 用有限体积积分格式．基于 SIMPLE 结构和向上差 分格式的 CFD 计算技术用来求解离散的平均 Navier－Stokes 方程计算采用空间交错有限体积结 构以保证计算的稳定性．关于饱和砂土流固耦合细 观力学模型更详细的计算步骤可参见有关文献［9］． 采用这种方法可以分析流体速度、流线以及压 力分布同时可以得到颗粒的运动和应力分布而这 些细观参量在实验室是很难或几乎不可能观察到 的． 2 数值实验程序 2∙1 数值试样生成 试样采用 PFC 随机生成器生成模拟容器大小 为15mm×27mm×90mm然后在重力作用下沉积 达到初始平衡这个过程一般需要较多的时步．模 拟容器侧壁采用不透水的刚性墙边界条件上下为 由 PFC 线墙模拟的透水边界条件图1是生成的模 型图表1是数值实验采用的指标参数．本构模型 采用 Herthz 接触模型．为了加速运算颗粒的剪切 模量采用了比实际土颗粒（如石英）小的值． 图1 砂沸现象模拟 PFC 数值模型试样 Fig．1 Sand samples by PFC for sand boil simulation 表1 颗粒流模型细观参数表 Table1 Parameters for numerical simulation 试样尺寸 粒径／ mm 摩擦因数／ 泊松比 孔隙 比 剪切模 量／MPa 流体 密度／ （kg·m －3） 粘滞系数／ （Pa·s） 有限体积 网格 计算时步／s PFC CFD 15mm×27mm×90mm 0∙3～0∙7 0∙5／0∙3 0∙42 2∙9 1000 10－3 3mm×3mm×3mm 2∙0×10－6 2∙0×10－3 生成的颗粒试样在重力和流体共同作用下砂土 颗粒沉积平衡循环计算消除不平衡力并使此时试 样的超孔压为零．初始模拟的沉积砂土颗粒约占整 个容器高度的67％左右． 2∙2 数值模拟过程 首先在模拟容器中充满水接着施加一个向上 ·474· 北 京 科 技 大 学 学 报 第30卷
第5期 刘洋等：砂沸现象的细观数值模拟及影响因素分析 .475 压力流，采用的方法是：在容器的底部施加一个逐 看出：在接下来的模拟过程中，留在上部的一部分颗 渐增加的压力值，直至达到临界水力梯度或更大，在 粒将一直留在顶部，这部分颗粒速度逐渐趋向于零； 此过程中始终保持上边界的压力为零.在每一个水 其余的大部分颗粒则开始下沉，并形成某种形式的 力梯度增量下流体和固相颗粒均达到稳定状态后再 对流运动， 施加下一个增量，每一步采用小的增量逐渐达到临 0.12 界或超临界条件，固相颗粒的指标特性通过沿轴向 0.10 的一系列测量球体（如图1）的平均值来反映，流体 0.08 的特性指标通过与某一测量球相对应位置的有限体 积网格中相关指标的值平均求得. 3数值模拟结果分析 0.02 e6 根据经典理论分析，产生砂沸的水力条件是：水 力梯度大于临界水力梯度，根据生成的试样，临界水 0.020 4 力梯度是1.17，P℉C数值模拟结果表明，当水力梯 时间s 度小于1.10时，颗粒仅产生很小的位移，体积变化 图3不同深度平均竖向流体速度时程曲线 也很小，这可以从孔隙率时程曲线看出（图2），图中 Fig-3 Curves of average vertical fluid velocity at different depths S1、S3、S5,S10代表不同深度的测量球（如图1）.此 时的孔隙水运动以层流为主，但已经偏离了达西 定律 1.0 S10 0.9 0.8 0.7 S3 0.5 图4不同时刻颗粒的演化图(XZ面) 0.3 Fig4 Evolution of particles at different time (XZ plane) 6 10 时间s 图2不同深度孔隙率时程曲线 Fig.2 Curves of porosity at different depths 图3是不同深度平均竖向流速时程曲线。当试 样水力梯度由1.10增加到1.12时，由图中看出，颗 粒被流体一颗粒的相互作用力带起，产生很大速度 和位移.颗粒的大位移运动使试样在空间上的孔隙 率产生突然、不均匀地增加，图3中L1、L3、L5、L10 是与测量球编号1,3,5,10相对应土层的位置，速 度的增加在图4~图6中也有反映. 001020304050.60.70.8091.0 图4是XZ面不同时刻(t=1,2,3,4,4.5,8s) 图5不同时刻孔隙率的演化图(XZ面) 颗粒的演化图，图5是对应时刻试样的孔隙率分布， Fig.5 Evolution of porosity at different time (XZ plane) 从图中可以看出，随着流体垂向速度的继续增加，上 部的颗粒到达顶部墙体，接着到达顶部的颗粒一部 图6和图7是在水力梯度为1.12时，给定时刻 分在重力和拽曳力的作用下反向沉积，颗粒到达顶 (t=1,2,3,4,4.5,8s)颗粒和流体的速度矢量图. 部墙体后，流体的速度也迅速降低，从图4中可以 孔隙流体择取高孔隙率的路径流动，使其流动过程
压力流．采用的方法是：在容器的底部施加一个逐 渐增加的压力值直至达到临界水力梯度或更大在 此过程中始终保持上边界的压力为零．在每一个水 力梯度增量下流体和固相颗粒均达到稳定状态后再 施加下一个增量每一步采用小的增量逐渐达到临 界或超临界条件．固相颗粒的指标特性通过沿轴向 的一系列测量球体（如图1）的平均值来反映．流体 的特性指标通过与某一测量球相对应位置的有限体 积网格中相关指标的值平均求得． 3 数值模拟结果分析 根据经典理论分析产生砂沸的水力条件是：水 力梯度大于临界水力梯度根据生成的试样临界水 力梯度是1∙17PFC 数值模拟结果表明当水力梯 度小于1∙10时颗粒仅产生很小的位移体积变化 也很小这可以从孔隙率时程曲线看出（图2）图中 S1、S3、S5、S10代表不同深度的测量球（如图1）．此 时的孔隙水运动以层流为主但已经偏离了达西 定律． 图2 不同深度孔隙率时程曲线 Fig．2 Curves of porosity at different depths 图3是不同深度平均竖向流速时程曲线．当试 样水力梯度由1∙10增加到1∙12时由图中看出颗 粒被流体－颗粒的相互作用力带起产生很大速度 和位移．颗粒的大位移运动使试样在空间上的孔隙 率产生突然、不均匀地增加图3中 L1、L3、L5、L10 是与测量球编号13510相对应土层的位置．速 度的增加在图4～图6中也有反映． 图4是 XZ 面不同时刻（ t＝12344∙58s） 颗粒的演化图图5是对应时刻试样的孔隙率分布． 从图中可以看出随着流体垂向速度的继续增加上 部的颗粒到达顶部墙体接着到达顶部的颗粒一部 分在重力和拽曳力的作用下反向沉积．颗粒到达顶 部墙体后流体的速度也迅速降低．从图4中可以 看出：在接下来的模拟过程中留在上部的一部分颗 粒将一直留在顶部这部分颗粒速度逐渐趋向于零； 其余的大部分颗粒则开始下沉并形成某种形式的 对流运动． 图3 不同深度平均竖向流体速度时程曲线 Fig．3 Curves of average vertical fluid velocity at different depths 图4 不同时刻颗粒的演化图（ XZ 面） Fig．4 Evolution of particles at different time （ XZ plane） 图5 不同时刻孔隙率的演化图（ XZ 面） Fig．5 Evolution of porosity at different time （ XZ plane） 图6和图7是在水力梯度为1∙12时给定时刻 （ t＝12344∙58s）颗粒和流体的速度矢量图． 孔隙流体择取高孔隙率的路径流动使其流动过程 第5期 刘 洋等： 砂沸现象的细观数值模拟及影响因素分析 ·475·
476 北京科技大学学报 第30卷 用的结果，反映了粒间接触力的减小 图9是不同深度平均配位数变化的时程曲线 平均配位数从开始的4.5左右减少到2左右（此时 颗粒作对流运动)·对于三维球体颗粒试样而言，配 位数大于4是稳定的结构 图6不同时刻试样的颗粒速度图(XZ面) Fig.6 Evolution of particle velocity at different time (XZ plane) 2 中所受阻力最小，能量消耗最少，同时其流动轨迹也 受颗粒运动变化的影响而不断变化，大部分区域的 流体流动虽然偏离了达西定律，但仍在层流范围内 流动，而有一些区域的流动则明显是紊流状态 时间s 图9不同深度配位数的时程曲线图 Fig.Curves of coordinate number at different depths 4影响因素分析 4.1水力梯度的影响分析 为了研究水力梯度的影响，进行了一系列的数 值模拟实验，水力梯度的值从1.0到1.3.当水力梯 度大于1.10时，所有数值模拟的结果在形式上都很 图7不同时刻试样的流体速度图(XZ面) 致，即试样内的颗粒处于某种悬浮状态和作某种 Fig-7 Evolution of fluid velocity at different time(XZ plane) 形式的对流运动，在试样的顶部会有一部分颗粒稳 定在此，这部分颗粒数量的多少和施加的水力梯度 图8是给定时间平均竖向有效应力以及孔隙水 大小有关系，水力梯度越大，停留在顶部的颗粒数量 压力一深度图.从图8中可以看出，平均有效应力随 就越多，图10是水力梯度为1.3试样稳态响应时 着水力梯度的增加逐渐减小，最终在到达临界水力 颗粒分布、孔隙率分布、颗粒流速以及水流速的模拟 梯度时趋于零.当水力梯度到达临界水力梯度后， 下部的土体首先产生了液化，接着向上发展。有效 结果.与图4相比，试样的稳态响应基本相同· 应力的减小是流体颗粒相互作用产生的拽曳力作 +1=05 +10s 010.45 0104s +1=1.0s 0.07 +1=1.0s 81=2.0s 01■2.0s 91=4.0s 1=4.0s E0.05 0.03 0.01 600 400 -200 200 400 图10水力梯度=1.30时的模拟结果 600 有效应力Pa 孔单水压力Pa Fig-10 Results of numerical simulation when i=1.30 图8有效应力和孔隙水压力一深度变化图 这里注意到临界流条件，在水力梯度达到1.10 Fig.8 Curves of effective stress and pore water pressure with depth 时，颗粒运动就开始加速，砂沸现象发生，这个数值
图6 不同时刻试样的颗粒速度图（ XZ 面） Fig．6 Evolution of particle velocity at different time （ XZ plane） 中所受阻力最小能量消耗最少同时其流动轨迹也 受颗粒运动变化的影响而不断变化．大部分区域的 流体流动虽然偏离了达西定律但仍在层流范围内 流动而有一些区域的流动则明显是紊流状态． 图7 不同时刻试样的流体速度图（ XZ 面） Fig．7 Evolution of fluid velocity at different time （ XZ plane） 图8 有效应力和孔隙水压力－深度变化图 Fig．8 Curves of effective stress and pore water pressure with depth 图8是给定时间平均竖向有效应力以及孔隙水 压力－深度图．从图8中可以看出平均有效应力随 着水力梯度的增加逐渐减小最终在到达临界水力 梯度时趋于零．当水力梯度到达临界水力梯度后 下部的土体首先产生了液化接着向上发展．有效 应力的减小是流体－颗粒相互作用产生的拽曳力作 用的结果反映了粒间接触力的减小． 图9是不同深度平均配位数变化的时程曲线． 平均配位数从开始的4∙5左右减少到2左右（此时 颗粒作对流运动）．对于三维球体颗粒试样而言配 位数大于4是稳定的结构． 图9 不同深度配位数的时程曲线图 Fig．9 Curves of coordinate number at different depths 4 影响因素分析 4∙1 水力梯度的影响分析 为了研究水力梯度的影响进行了一系列的数 值模拟实验水力梯度的值从1∙0到1∙3．当水力梯 度大于1∙10时所有数值模拟的结果在形式上都很 一致即试样内的颗粒处于某种悬浮状态和作某种 形式的对流运动．在试样的顶部会有一部分颗粒稳 定在此这部分颗粒数量的多少和施加的水力梯度 大小有关系水力梯度越大停留在顶部的颗粒数量 就越多．图10是水力梯度为1∙3试样稳态响应时 颗粒分布、孔隙率分布、颗粒流速以及水流速的模拟 结果．与图4相比试样的稳态响应基本相同． 图10 水力梯度 i＝1∙30时的模拟结果 Fig．10 Results of numerical simulation when i＝1∙30 这里注意到临界流条件在水力梯度达到1∙10 时颗粒运动就开始加速砂沸现象发生这个数值 ·476· 北 京 科 技 大 学 学 报 第30卷
第5期 刘洋等：砂沸现象的细观数值模拟及影响因素分析 .477. 比连续介质力学分析得到的值1.17要小约6%.这 为3.0GPa,较“软”的颗粒试样刚度为3.0MPa 个差别实际是反映了模拟的动力特性，与采用细观 图11是在临界水力梯度下两种试样达到稳定状态 模型的计算模式有关 时的模拟结果，从宏观响应来看，除局部区域有细 4.2颗粒刚度的影响分析 小差别外，颗粒刚度对砂沸现象的模拟结果影响不 为了研究颗粒刚度对数值模拟结果的影响，采 大·颗粒刚度的减小会使计算效率有很大的提高， 用两个结构相同但剪切模量不同的试样进行了临界 因此采用较“软”颗粒进行数值模拟，可以加速计算 水力梯度下的数值模拟，刚度较大的试样剪切模量 (a 图11不同颗粒刚度的模拟结果.(a)剪切模量3.0MPa:(b)剪切模量3.0GPa Fig.11 Results of numerical simulation at different particle shear moduli:(a)3.0 MPa:(b)3.0GPa (汪闻韶.土的动力强度与液化特性.北京：中国电力出版社， 5结论 1997) [2]Wang W S.Distinction and interrelation between liquefaction. (I)采用P℉C数值方法模拟砂沸现象，能够从 state of limit equilibrium and failure of soil mass.Chin Geotech 细观的角度分析液化时砂沸产生的细观机理 Eng,2005,27(1):1 (2)数值模拟结果显示，在临界流或超临界流 (汪闻韶.土体液化与极限平衡和破坏的区别和关系·岩土工 条件下，试样最终的稳态响应分成两部分：一部分颗 程学报，2005,27(1)：1) 粒停留在试样的顶部，这部分颗粒数量正比于施加 [3]Bardet J P,Kapuskar M.Liquefaction sand hoils in San Francisco 的水力梯度的大小；另一部分颗粒则在试样底部形 during 1989 Loma Prieta Earthquake.JGeotech Eng.1993.119 (3):543 成某种形式的对流运动， [4]Hokzer T L.Clark MM.Sand boils without earthquakes.Geolo- (③)采用的细观数值模型分析所得的临界水力 9y,1993.21(10):873 梯度比连续介质力学分析得到的值要小约6%，这 [5]Tsuji Y,KawaguchiT,Tanaka T.Discrete particle simulation of 个差别反映了模型的动力特性·在临界流或超临界 two-dimensional fluidized bed.Powder Technol,1993.77:79 流条件下，砂沸时水的流动除局部区域外基本上仍 [6]Hooms BP H.Kuipers JA M.Briels W J.Discrete particle sim- ulation of bubble and slug formation in a two dimensional gas flu 属于层流范围，但已偏离了达西定律. idized bed:a hard-sphere approach.Chem Eng Sci.1996.51 (4)颗粒刚度对相同结构特性的试样最终稳态 (1):99 响应影响不大；但采用较“软”的颗粒进行模拟，可以 [7]Kafui K D.Thornton C.Discrete particle continuum fluid model- 加速计算，提高计算效率. ing of gad-solid fluidized beds.Chem Eng Sei.2002.57(13): (5)需要指出的是：本文在模拟砂沸现象时，采 2395 [8]Jackson R.Locally averaged equations of motion for a mixture of 用的方法是通过建立的砂沸数值模型，在给定超孔 identical spherical particles and a Newtonian fluid.Chem Eng 隙水压力增量的条件下，使之逐渐达到临界流状态 Sci,1997,52(15):2457 从而产生砂沸；进一步的工作研究在地震荷载作用 [9]Liu Y.Meso scale Mechanism of Sand Liquefaction and Numer- 下，砂沸现象随机产生和发展的全过程， ical Simulation Dissertation ]Shanghai:Tongii University. 2006 参考文献 (刘洋.砂土液化破坏的细观力学机制与数值模拟[学位论 [1]Wang WS.Dynamic Strength and Liquefaction Properties of 文]上海：同济大学，2006) Soil.Beijing:China Electric Power Press,1997
比连续介质力学分析得到的值1∙17要小约6％．这 个差别实际是反映了模拟的动力特性与采用细观 模型的计算模式有关． 4∙2 颗粒刚度的影响分析 为了研究颗粒刚度对数值模拟结果的影响采 用两个结构相同但剪切模量不同的试样进行了临界 水力梯度下的数值模拟．刚度较大的试样剪切模量 为3∙0GPa较“软” 的颗粒试样刚度为3∙0MPa． 图11是在临界水力梯度下两种试样达到稳定状态 时的模拟结果．从宏观响应来看除局部区域有细 小差别外颗粒刚度对砂沸现象的模拟结果影响不 大．颗粒刚度的减小会使计算效率有很大的提高 因此采用较“软”颗粒进行数值模拟可以加速计算． 图11 不同颗粒刚度的模拟结果．（a） 剪切模量3∙0MPa；（b） 剪切模量3∙0GPa Fig．11 Results of numerical simulation at different particle shear moduli：（a）3∙0MPa；（b）3∙0GPa 5 结论 （1） 采用 PFC 数值方法模拟砂沸现象能够从 细观的角度分析液化时砂沸产生的细观机理． （2） 数值模拟结果显示在临界流或超临界流 条件下试样最终的稳态响应分成两部分：一部分颗 粒停留在试样的顶部这部分颗粒数量正比于施加 的水力梯度的大小；另一部分颗粒则在试样底部形 成某种形式的对流运动． （3） 采用的细观数值模型分析所得的临界水力 梯度比连续介质力学分析得到的值要小约6％这 个差别反映了模型的动力特性．在临界流或超临界 流条件下砂沸时水的流动除局部区域外基本上仍 属于层流范围但已偏离了达西定律． （4） 颗粒刚度对相同结构特性的试样最终稳态 响应影响不大；但采用较“软”的颗粒进行模拟可以 加速计算提高计算效率． （5） 需要指出的是：本文在模拟砂沸现象时采 用的方法是通过建立的砂沸数值模型在给定超孔 隙水压力增量的条件下使之逐渐达到临界流状态 从而产生砂沸；进一步的工作研究在地震荷载作用 下砂沸现象随机产生和发展的全过程． 参 考 文 献 ［1］ Wang W S．Dynamic Strength and L iquef action Properties of Soil．Beijing：China Electric Power Press1997 （汪闻韶．土的动力强度与液化特性．北京：中国电力出版社 1997） ［2］ Wang W S．Distinction and interrelation between liquefaction state of limit equilibrium and failure of soil mass．Chin J Geotech Eng200527（1）：1 （汪闻韶．土体液化与极限平衡和破坏的区别和关系．岩土工 程学报200527（1）：1） ［3］ Bardet J PKapuskar M．Liquefaction sand boils in San Francisco during1989Loma Prieta Earthquake．J Geotech Eng1993119 （3）：543 ［4］ Holzer T LClark M M．Sand boils without earthquakes．Geology199321（10）：873 ［5］ Tsuji YKawaguchi TTanaka T．Discrete particle simulation of two-dimensional fluidized bed．Pow der Technol199377：79 ［6］ Hooms B P HKuipers J A MBriels W J．Discrete particle simulation of bubble and slug formation in a two-dimensional gas fluidized bed：a hard-sphere approach．Chem Eng Sci199651 （1）：99 ［7］ Kafui K DThornton C．Discrete particle-continuum fluid modeling of gad-solid fluidized beds．Chem Eng Sci200257（13）： 2395 ［8］ Jackson R．Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid． Chem Eng Sci199752（15）：2457 ［9］ Liu Y．Meso-scale Mechanism of Sand L iquef action and Numerical Simulation ［ Dissertation ］．Shanghai：Tongji University 2006 （刘洋．砂土液化破坏的细观力学机制与数值模拟 ［学位论 文］．上海：同济大学2006） 第5期 刘 洋等： 砂沸现象的细观数值模拟及影响因素分析 ·477·