饱和砂土流固耦合细观数值模型及其在液化分析中的应用

　2009 年 2 月 水 　　利 　　学 　　报 SHUILI 　　XUEBAO 第 40 卷 　第 2 期 收稿日期 :2008201225 基金项目 :国家自然科学基金资助项目 (50808016) 作者简介 :刘洋 (1979 - ) ,男 ,江苏徐州人 ,博士 ,副教授 ,主要从事饱和砂土液化方面的研究工作。E2mail :ly2ocean @sohu. com 文章编号 :055929350 (2009) 0220250207 饱和砂土流固耦合细观数值模型及其在液化分析中的应用 刘 　洋1 ,周 　健2 ,付建新1 (11 北京科技大学 土木与环境工程学院 ,北京　100083 ; 21 同济大学 地下建筑与工程系 ,上海　200092) 摘要 :在饱和砂土连续力学模型的基础上 ,建立了一个饱和砂土固相颗粒和液相流体耦合的细观力学模型 ,固相颗 粒采用离散元的颗粒流理论模拟 ,液相流体通过求解平均 Navier2stokes 方程的计算流体动力学技术计算。该细观 模型用二维渗流问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 进行了验证 ,并将其应用于饱和砂土的液化分析。数值模拟结果显示 ,采用的模型可以描述 从低雷诺流到高雷诺流范围很大的流体运动。液化分析 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 明 :液化首先在试样的表层产生 ,接着向深部发展 ;而超 孔隙水压力则从底部开始消散 ,然后逐渐向上发展 ,直至完全消散。砂土液化过程中土体孔隙率不断变化 ,渗透系 数也随之改变。 关键词 :饱和砂土 ; 流固耦合 ; 离散单元法 ; 液化 中图分类号 :TU441 + 14 文献标识码 :A 饱和砂土是由固相颗粒集合体以及充填其中的流体组成的混合物 ,基于颗粒之间以及颗粒 - 流体 之间的相互作用 ,饱和砂土呈现出高度复杂的非线性力学特性。传统的饱和砂土的连续耦合模型通常 基于混合物理论、Biot 理论或微观力学运动方程的均一化来建立[1 - 2 ] 。这些不同形式的方程都需要某种 形式的本构关系来描述固相的应力应变关系。这些本构关系中 ,以帽子模型、多重屈服面塑性模型、边 界面塑性模型最为常用。上述用于模拟饱和砂土的连续力学模型都是基于固液两相的唯象描述。在场 方程和本构关系中有很多假设 ,如初始局部孔隙率通常认为与总体平均孔隙率相等 ,相之间的动量交换 遵循 Darcy 定律 ,而渗透系数是独立于应变的等等。实际上 ,Darcy 定律是低阶的 Navier2Stokes 方程 ,仅 当孔隙流是薄层流、惯性力可以忽略的情况下才有效[3 ] 。在高的水力梯度下将产生非薄层流 ,随着砂土 骨架大的体积变形产生的孔隙率变化也可能导致偏离 Darcy 定律。为了得到对饱和砂土液化现象比较 真实的模拟结果 ,Manzari 等[4 ]建议使用渗透函数 ,即土的渗透性随时间而增加。这种假定和修正提高 了这些模型的预测精度 ,但仍然存有很多问题。基于此 ,本文采用一个比较实用的饱和砂土的细观力学 模型 ,通过多尺度、多相耦合来反映离散骨架的变形、孔隙流体的流动以及它们之间的相互作用。 1 　饱和砂土的微观模拟 饱和砂土的固相颗粒 ,可以采用离散颗粒的集合来有效地模拟 ,在岩土力学中最常用的方法就是离 散单元法 (DEM) [5 ] 。对于液相流体的模拟 ,早期的研究主要是解慢速流绕球或圆柱体的 stokes 方程 ,得 到流体作用于单个球体或圆柱体的拽曳力。Meegoda 等[6 ] 通过有限元求解 Navier2Stokes 方程 ,建议了一 个各向异性介质的渗透性公式。在他的研究中 ,土被假定为一个由给定方式周期排列的空间周期介质。 光滑粒子流体动力学方法 (SPH) [7 ] 、格栅波尔兹曼方法 (LB) [8 - 9 ] 也用来研究多孔介质中的孔隙流问题。 Cundall[10 ]曾提出了一个简化的流体 - 颗粒耦合方法 ,在这个方法中 ,颗粒之间的孔隙被认为是通过颗 —052— 粒接触之间的管子连接的流域 ,通过建立网络来表征流域的体积和管子的长度和直径。周健等[11 ] 采用 了这种方法对渗流问题进行了模拟。这种方法对于少颗粒系统比较简便 ,但三维情况下孔隙体积的计 算将变得很繁琐。 从上述不同的方法来看 ,完全微观的方法的应用范围有限 ,其用处多限于研究孔隙中流体流动的局 部力学机制或用来验证已有的流体颗粒相互作用的经验模型。 2 　饱和砂土流固耦合的细观数值模型 211 　数值建模思路　基于上述分析 ,考虑在比孔隙尺度稍微粗糙一点的尺度上对流体进行模拟 ,这种 模拟采用空间域的平均流体速度来代替每一点的流体速度 ,这个域包含一定数目的颗粒但相对于整个 混合物空间来说又是很小的。流体施加于单个颗粒的力与这个局部的平均流体速度以及孔隙率相关 联。Tsuji [12 ]首先将这种离散元模拟和均一化流体公式相结合的方法用于分析流化床中泡沫的形成。 从那以后 ,这种数值模型被广泛地用于一些化学过程的模拟。本文即采用这个思路来建立饱和砂土固 液耦合的细观力学模型。 具体建模思路如下 :采用离散元法的颗粒流理论模拟饱和砂土中砂颗粒的力学行为 ,采用均一化流 体 CFD 计算技术模拟孔隙流体的运动。对于固相颗粒 ,通过求解运动方程模拟颗粒运动 ;对于液相水 , 通过求解平均的 Navier2Stokes 方程模拟孔隙流体的运动。这两组方程通过施加于颗粒表面非滑动的边 界条件相耦合。采用这种方法可以分析流体速度、流线以及压力分布 ,同时可以得到颗粒的运动和应力 分布 ,而这些细观参量在实验室很难或几乎不可能观察到的。 212 　液相水数值模型　平均的 N2S 方程提供了一种可计算的细观途径 ,能够定量地描述孔隙流体特性 以及流体2颗粒的相互作用 ,其形式如下[13 ] :9 n9 t + Δ·( n€vf ) = 0 (1) ρf 9 ( n€vf )9 t + Δ·( n€vf€vf ) = Δ( n€σf ) - ∑p ∫Spσf ·ng (| x - y | ) d sy + nρf f g (2) 式中 : n = n ( x) = ∫ RVf g ( r) d vy ,为 x 附近的局域孔隙率 ; g ( r) = g (| x - y | ) 为一个单调增加的权函数 , r 为局域半径 ;€vf 和€σf 为平均流体速度和平均应力张量 ;ρf 为流体质量密度 ; f g 为重力加速度矢量 ; n 为 颗粒表面的单位外矢量 ;下标 p 为平均体积中的颗粒 ; S p 为颗粒的表面积。 Jackson 认为式 (2)可以写成如下形式 : ρf 9 ( n€vf )9 t + Δ( n€vf€vf ) = n Δ€τf - €f i + nρf f g (3) 式中 :€f i 为平均单位体积流体 - 颗粒相互作用力 ; €τf 为流体平均应力张量 ,这个张量除了包括平均常规 应力 €σf 外还包括了类似于湍流中的雷诺应力项。 213 　流体 - 颗粒相互作用　流体 - 颗粒相互作用拽曳力可以用一些半经验的公式来计算[14 - 15 ] 。 当 n ≤018 时 ,采用 Ergun[14 ]的计算公式 :€f i = (1 - n) 150μf (1 - n) n…d2p + 1175ρf | €vf - €vp |…dp (€vf - €vp ) (4) 　　当 n ≥018 时 ,应用 Wen 和 Yu 的公式[15 ] :€f i = (1 - n) 34 Cd Re p n - 2165μfρf | €vf - €vp |…dp (€vf - €vp ) (5) 式中 :μf 为流体黏性系数 ; Rep = nρf df | €vf - €vp | Πμf ,为具有平均速度为 €vp 和平均半径为…dp 的散粒体的雷 诺系数 ; Cd 是拽曳力系数。 —152— 与连续介质模型中通常使用的两相流公式相比 ,式 (4) 、式 (5)中的拽曳力项代表了土的渗透性是骨 架孔隙率以及固相和孔隙流体性质的函数。即与目前的连续介质模型不同 ,在模拟过程中土渗透性的 时空改变是随着土和孔隙流体的性质而改变的。 214 　考虑流体作用的砂土颗粒数值模拟　饱和砂土颗粒仍采用离散元模拟 ,当有流体作用时 ,要考虑 流体颗粒相互作用。单个颗粒运动仍遵循牛顿第二定律 ,颗粒的平移速度 vp 和转动速度 …ωp 通过下式 计算 : mpÛvp = mpvg + ∑ c f c +∫ Sp σf ·nd Sp ; Ipω · - p = ∑ c rc ×f c +∫ Sp r ×(σf ·n) d Sp (6) 式中 : mp 、Ip 、Sp 分别为质量、惯性力矩以及颗粒表面积 ; f c 为接触处 c 的粒间相互作用力 ; rc 为连接颗 粒中心和接触力位置的矢量 ; r 为半径矢量。式 (6)的积分项包括了流体作用于颗粒 p 上的合力和合力 矩。按照牛顿第三定律 ,流体 - 颗粒相互作用力应该相应于式 (4) 、式 (5) 中的流体2颗粒相互作用项€f i 。 采用平均流体应力 €τf 和流体2颗粒相互作用力€f i ,对于给定体积 Vp 内的颗粒 p ,方程 (6)的积分项可以写 成 : ∫ Sp σf ·ndSp = - Δpf + €f i 1 - n Vp ;∫ Sp r ×(σf ·n) dSp = rg × Δ·€τf Vp (7) 式中 : ΔPfVp 为浮力 ; rg 为连接颗粒形心和质心的矢量。 215 　算法结构　流固耦合方程 (1) 、(3) 、(6)的解采用显式的结构。在求解耦合方程过程中的第 i 时步 ( i = 1 ,2 ,3 ⋯) ,在接触力、外部荷载以及流体拽曳力作用下颗粒的位移采用离散元技术计算。显式的中 心差分来积分固相的运动方程 (6) ,计算稳定所需的时步由粒间接触的刚度和阻尼系数决定。孔隙流体 的压力和速度由平均的N2S 方程求解 ,求解采用有限体积积分格式。基于 SIMPLE结构和向上差分格式 的 CFD 计算技术用来求解离散的平均 N2S 方程。采用空间交错有限体积结构以保证计算的稳定性 ,即 流体的平均孔隙水压力和固相的平均指标的值在差分中心位置 (体积或面积的中心位置) 计算 ,流体的 速度在体积 (或面积)边的中心位置交错计算。 图 1 　PFC2CFD 二维 数值模拟试样 3 　二维验证 311 　试样生成 　为了验证上述模型进行了数值模拟[10 ] ,模型大小 240mm ×840mm , 分成细砂和粗砂两种试样 ,平均颗粒半径分别为 015mm 和 3mm。每种类型的试样 由 3 种不同的孔隙率组成 ,试验颗粒粒径从最小粒径到最大粒径均匀分布 ,每个试 样的粒径范围不同 ,但平均粒径相同 ,模型参数如表 1 所示。 颗粒试样首先在 PFC2D随机生成器中生成 ,然后在重力作用下平衡 ,调整初始 孔隙率反复试验直到得到 3 种孔隙率分别为 0115、0120、0125 ,共生成 6 个标准试 样。图 1 是其中的一个代表性的试样结构 ,试样两边是刚性的墙体 ,下边界采用点 墙来模拟试验室的透水石材料。图中网格是 CFD 流体计算网格 ,网格的选取大小 适当 ,既包含一定数量的颗粒又不太大 ,流体计算域大于颗粒范围。 表 1 　PFC2CFD 数值试验基本输入参数 试样类型 试样尺寸Π mm2 平均粒径Π mm 摩擦系数 接触刚度 (法向Π切向)ΠMPa 流体密度Π(kgΠm3) 粘滞系数Π(Pa·S) 细砂 240 ×840 015 016 500Π500 1000 10 - 3 粗砂 240 ×840 3 016 500Π500 1000 10 - 3 312 　数值模拟过程 　对砂样分别在顶部和底部施加流体压力边界条件 ,使之产生压力差 ,在不同的水 力梯度条件下进行 PFC2CFD 耦合运算 ,使试样最终达到稳定渗流状态 ,平衡稳定的标准是计算所得的 平均流速稳定并趋于一定值。 —252— 313 　模拟结果分析 　模拟结果如图 2 和图 3 所示。从图中可见 ,在低雷诺系数 ( Rep < 1) 条件下 ,流速 较小 ,流体流速和水力梯度呈明显的线性关系 ,服从 Darcy 定律。当雷诺系数较高时 ,流体流动已经偏 离了层流条件。随着雷诺系数的增加 ,流体流动明显偏离了 Darcy 定律 ,流速和水力梯度成明显的非线 性关系。数值计算结果表明 ,在高雷诺系数非层流条件下 ,采用的 PFC2CFD 耦合模型能够模拟多孔介 质中的非 Darcy 流。从二维的验证来看 ,采用的饱和砂土流固耦合细观模型可以描述从低雷诺流到高 雷诺流范围很大的流体运动。 图 2 　低雷诺系数二维数值模拟结果 　 图 3 　高雷诺系数二维数值模拟结果 　 4 　饱和砂土液化分析 在地震荷载作用下 ,饱和砂土颗粒重组 ,土层结构发生变化会产生孔隙率的不均匀变化和孔隙水压 力的升高 ,当流体作用于颗粒的力达到有效应力时将产生液化或失稳。此时土体变成流态 ,孔隙水流动 也变得异常复杂 ,不再是简单的 Darcy 流。下面采用上述细观力学模型对饱和砂土液化过程中的颗粒 流体响应进行模拟分析。 图 4 　数值模拟结果数值模拟输入地震波时程曲线 411 　数值模型 　颗粒试样由 PFC随机生成器生成[10 ] ,然 后在给定重力场作用下平衡 ,模型试样尺寸 140mm × 60mm×250mm ,其中固相高 H = 200mm ,模型两侧采用周 期边界条件。颗粒平均直径 6mm ,采用 Heitz 接触模 型[10 ] ,剪切模量 290MPa ,基底用不透水的墙来模拟 ,表面 自由排水 ,基底输入正弦地震波如图 4 ,震动持时 15s。 412 　渗透系数变化　定义渗透系数 k = c…d2p n3 / (1 - n) 2 , 其中 c = 01003～010055 ,为与土性有关的常数 , dp为有限体积范围内颗粒的平均直径。砂土液化过程中 土体孔隙率不断变化 ,渗透系数也随之改变 ,如前所述 ,采用的细观模型能够反映这一变化。图 5 是不 同深度孔隙率和渗透系数的时程曲线。从图中可以看出 ,液化后浅部的土层显著地变密了。随着孔隙 率的变化 ,土体的局部渗透系数也随之变化。 413 　超孔隙水压力　连续介质模型一般将动孔压的增长消散变化与土体的体积变形相关联 ,提出某种 形式的孔压增长模式。饱和砂土液化过程是土体有效应力衰减的过程 ,振动过程中超孔隙水压力的反 应是土体有效应力变化最直接的因素 ,下面就超孔压变化进行分析。 图 6 是震动过程试样不同深度超孔隙水压力比 (超孔隙水压力与竖向有效应力之比) 的增长曲线。 从图中可见 ,随着震动的进行 ,超孔隙水压力比逐渐增加达到 110 左右 ,超孔隙水压力比从试样表面 (0115 H)到底部 (019 H) 依次趋近于 110 ,试样表面在震动 310s 后到达 110 ,底部约 710s 到达 110。超孔 隙水压力的演化表明 ,试样在模拟地震荷载作用下产生了液化 ,土体内部液化的发生不是同时的 ,模型 液化首先从浅部开始 ,逐渐向深部发展。 图 6 中可注意到超孔隙水压力在试样底部增加相对缓慢 ,这主要是因为在试样底部的砂层密实度 —352— 图 5 　不同深度孔隙率及渗透系数的时程曲线 　 较大所致。同时也可以看出在试样表面 (0115 H) ,除短暂时刻外 ,超孔压比均未达到 1 ,值大约在 012～ 018 左右。原因是数值试样表面采用的是自由排水边界条件 ,限制了表层超孔隙水压力的进一步增加。 图 7 是试样不同深度超孔隙水压力的消散曲线 ,从图中可以看出 ,与超孔压增长模式不同的是 ,超 孔压不是从表面而是从底部开始消散 ,然后逐渐向上发展 ,最后完全消散 ,整个试样超孔压的消散大约 在振动结束后的 45～50s 内完成。 图 6 　振动过程不同深度超孔隙水压力时程曲线 　 414 　平均配位数的演化　粒间接触的变化是颗粒受力变化的间接反映 ,液化时颗粒处于似悬浮状态 , 粒间接触平均配位数的变化也反映着液化过程的发展。不同深度位置的平均配位数采用相应的测量球 来测量 ,图 8 是在振动过程中深度分别为 0115 H、014 H、017 H、019 H 处平均配位数的演化曲线 ,从图中 可以看出 ,随着地震荷载的输入 ,各位置平均配位数都逐渐降低 ,0115 H 处平均配位数首先减小至 4 以 下 ,接着随着时间的增加 ,下面各层土体平均配位数依次下降至小于 4 的值 ,而 4 或者大于 4 是三维球 体颗粒稳定所必需的配位数。这说明自表层开始 ,试样产生了液化 ,各层土体液化发生的程度也有不 同 ,0115 H 处 ,配位数减小至最低 310 左右 ,维持 5s 左右又重新恢复值 410 以上。014 H 处土体平均配位 数减小地最为明显 ,最低至大约在 016 左右 ,在振动阶段基本维持在此水平 ,这说明此处土体基本处于 悬浮状态 ,颗粒之间基本脱离接触。017 H 处平均配位数降低至 112 左右。基底附近 019 H 处土层平均 配位数减小到 118 左右 ,并且很快就恢复到了大于 410 的值。试样不同深度细观配位数的演化也说明 试样液化首先在表层产生 ,然后向深部发展。 —452— 图 7 　不同深度超孔隙水压力消散曲线 　 图 8 　不同深度平局配位数的演化曲线 　 5 　结论 本文在分析饱和砂土连续力学模型的基础上 ,建立了一个描述饱和砂土固相颗粒和液相流体耦合 的细观力学模型 ,固相颗粒采用离散元法的颗粒流理论模拟 ,液相流体通过求解平均 N2S 方程的计算流 体动力学技术计算。利用建议的模型对二维渗流问题进行了验证 ,并将其应用于饱和砂土的液化分析 , 初步的研究得出以下几点结论 : (1)与传统的连续介质力学模型相比 ,采用的饱和砂土流固耦合的细观 模型可以模拟在流体作用下固相的细观力学响应。在数值模拟过程中介质渗透系数不是时空固定的 , 而是随着土和孔隙流体性质的改变而改变 ; (2)数值模拟结果表明 ,采用的饱和砂土流固耦合细观模型 可以描述从低雷诺流到高雷诺流范围很大的流体运动 ,可以用于包括渗流、饱和砂土液化等岩土工程问 题分析 ; (3)细观数值模拟结果显示 ,液化首先在试样的表层产生 ,接着向深层发展 ;而超孔隙水压力则 从底部开始消散 ,然后逐渐向上发展 ,直至完全消散。 参 　考 　文 　献 : [ 1 ] 　Lewis R W , Schrefler B A. The finite element method in the deformation and consolidation of porous media[M]. Wiley , —552— New York , USA ,1987. [ 2 ] 　Zienkiewicz O C , Chan A H C , Pastor M , Schrefler B A , Shiomi T. Computational Geomechanics with Special Reference to Earthquake Engineering[M]. Wiley : England , 1998. [ 3 ] 　Trussell R R , Chang M. Review of flow through porous media as applied to head loss in water filters [J ] . Journal of Environmental Engineering 1999 , 125 (11) :998 - 1006. [ 4 ] 　Mazari M , Arulanandan K. Numerical prediction for model no 1[ C]ΠΠArulananadan K, Scott R , editors ,Verification of numerical procedures foe the analysis of soil liquefaction problems. Proceedings of the International Conference on the Verification of numerical procedures for the analysis of soil liquefaction problems. Davis , California , 1993 : 179 - 785. [ 5 ] 　Cundall P A. Strack O D L. A discrete numerical model for granular assemblies[J ] . Geotechnique ,1979 , 29 (1) :47 - 65. [ 6 ] 　Meegoda M , King I , Arulanandan K. An expression for the permeability of anisotropic granular media [J ] . International Journal for Numerical and Analytical Methods in Geomechanics , 1989 , 13 : 575 - 598. [ 7 ] 　Potapov A , Hunt M , Campbell C. Liquid2solid flows using smoothed particle hydrodynamics and the discrete element method[J ] . Powder Technology. 2001 , 116 :204 - 213. [ 8 ] 　Ladd A. Numerical simulation of particulate suspensions via a discretized boltzmann2equation : 1. theoretical result [J ] . Journal of Fluid Mechanics , 1994 ,271 : 285 - 291. [ 9 ] 　Ladd A. Numerical simulation of particulate suspensions via a discretized boltzmann2equation : 2numerical result [J ] . Journal of Fluid Mechanics , 1994 ,271 : 311 - 317. [10 ] 　Particle Flow Code in 2 Dimensions , PFC2D Version 3. 1 , Itasca Consulting Group Inc , Minneapolis , 2004. [11 ] 　周健 ,张刚 ,孔戈. 渗流的颗粒流细观模拟[J ] . 水利学报 ,2006 , 37 (1) :28 - 32. [12 ] 　Tsuji Y, Kawaguchi T , Tanaka T. Discrete particle simulation of two2dimensional fluidized bed[J ] . Powder Technology , 1993 , 77 :79 - 87. [13 ] 　Jackson R. Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[J ] . Chemical Engineering Science , 1997 , 52 (15) :2457 - 2469. [14 ] 　Ergun S. Fluid flow through packed columns[J ] . Chemical Engineering Progress , 1952 ,43 (2) :89 - 94. [15 ] 　Wen C Y, Yu Y H. Mechanics of fluidization[J ] . Chemical Engineering Progress Symposium Series , 1966 ,62 (2) :100 - 113. Fluid2particle coupled model for saturated sand and its application to liquefaction analysis LIU Yang1 , ZHOU Jian2 , FU Jian2xin1 (1. Beijing University of Science and Technology , Beijing 　100083 , China ; 2. Tongji University , Shanghai 　200092 , China) Abstract :Based on the continuum model for analysis of saturated porous media , a fluid2particle coupled model was established to simulate the meso - scale pore water flow and micro2scale solid deformation of saturated soils. The particles were simulated by discrete element method , while the fluid flow was simulated by solving the average Navier2Stokes equations. The validity of this model was verified by a two2dimension flow case study and the model was applied to analyze the liquefaction phenomena of saturated sand. The numerical results indicate that the adopted model can be used to simulate the flow in a wide range of Reynolds number. The simulation results of liquefaction analysis show that the liquefaction begins from the surface and propagates downward ,while the dissipation of pore water pressure begins from the bottom and propagates upward until vanished. The spatial and temporal changes of soil porosity and the variation of permeability coefficient in the process of liquefaction accord with this phenomenon. Key words :saturated sand ; fluid2solid coupled model ; discrete element method ; liquefaction (责任编辑 :李 　琳) —652—