含微生物气泡高饱和砂土循环三轴剪切CFD-DEM模拟

2021-05-13蒋明镜沈振义谭亚飞鸥

地震工程学报 2021年2期

蒋明镜,沈振义,刘俊,杨涛,谭亚飞鸥

(1. 天津大学建筑工程学院土木系, 天津 300072; 2. 天津大学水利工程仿真与安全国家重点实验室, 天津 30007;3. 同济大学地下建筑与工程系, 上海 2000921; 4. 同济大学土木工程防灾重点实验室, 上海 200092;5. 上海理工大学环境与建筑学院, 上海 200093)

0 引言

砂土液化是指土体受到振动作用,孔隙水在压力作用下力图向上排出,土颗粒受重力影响向下沉落,使土颗粒处于悬浮状态,导致土体抗剪强度完全消失的现象。砂土液化现象是地震灾害的主要形式之一,通常会引起建筑物基础的不均匀沉降及结构的破坏,甚至造成严重灾害和人员伤亡。如1975年的海城地震、1976年的唐山地震和2008年的汶川地震都引起了大区域砂土的液化,许多建筑物由于基础的不均匀沉降遭到不同程度的破坏。目前,减少地基液化主要有换填、强夯地基土的物理处理方法和使用人工化学材料的化学处理方法,但物理方法能量效率低下、化学材料污染环境。因此,探寻新的提高砂土地基抗液化能力的方法具有重要意义。微生物灌浆技术作为一种防止砂土液化的方法,为地基处理提供了新的选择。

国内外学者对砂土液化的研究,大多集中在室内试验和数值模拟。在室内试验方面,Xia和Hu[1]对四组不同饱和度和五组施加不同反压的试样进行了循环三轴试验。结果表明:饱和度对抗液化能力有明显的影响,反压也对抗液化能力有影响。但此试验是通过反压技术来提高试样的饱和度,未能准确探究饱和度对液化的影响。Okamura和Soga[2]通过控制饱和度、初始围压和初始孔隙压力压三个参数,进行了一系列三轴试验。结果表明:饱和度对抗液化强度有显著影响。Tsukamoto[3]对粉砂进行了不排水抗剪强度的研究。结果表明:部分饱和和非饱和粉砂的抗液化能力较饱和粉砂明显提高。He等[4]通过反硝化反应生成的氮气制取含微生物气泡的高饱和砂土试样,进行了振动模拟台试验。结果表明:微生物产生的气泡可以提高土体的抗液化程度。然而,对于微生物气泡改善土体工程力学性质的微观机理认识尚浅,有必要对其微观加固机理进行详细研究。

针对室内试验在试样制取耗时长、成本高和重复性差上的不足,许多学者开始采用数值模拟方法对砂土液化微观机理进行研究。离散单元法最早由Cundall和Strack[5]提出,其用微观信息解释了离散颗粒材料的宏观力学特性。随后发展出考虑颗粒粗糙接触的二维抗转动模型[6]。目前,颗粒模型进入三维时代,考虑抗转动和扭转的三维模型[7]也被提出。计算流体动力学是流体力学的一个分支,运用特定的数值方法和算法求解与流体流动相关的问题[8]。Tsuji等[9]1993年首次将CFD嵌入DEM源程序,模拟固体颗粒中气体流动行为。结果表明:气流速度在气泡出现和压力的波动频率改变时与理论解较为吻合。Hoomans等[10]用该方法模拟气泡的形成,并采用硬球模型代替软球模型,对气体和颗粒之间的相互作用采用N-S方程描述,模拟结果与方程求得的结果能很好吻合。Wei等[11]采用离散单元法对不排水循环荷载作用下的砂土液化过程进行了数值模拟,通过对配位数等的变化的分析,探究了砂土液化机理。蒋明镜和张望城[12]在传统的CFD流体运动控制方程和大量室内试验结果的基础上,引入了考虑流体弱可压缩性的流体状态方程,建立了流体密度与压力时间的关系,利用FISH语言将CFD与PFC2D离散元计算耦合起来,用以对饱和砂土的固结不排水剪切试验进行模拟。对于微生物岩土技术抗液化能力的试研究,许多学者的研究重点依旧集中于生物胶结,对微生物气泡的重视程度不足,鲜有微生物气泡数值模拟的成果。因此,对含有生物气泡非饱和砂土开展CFD与DEM耦合的液化分析很有必要。

本文基于CFD-DEM耦合程序,运用PFC3D5.0商业软件,模拟了不同生物气泡含量的高饱和砂土的固结不排水循环剪切。在模拟过程中,通过跟踪试样的力链分布、液化振次、孔压比、轴向应变和配位数在加载过程中的变化情况,从宏微观角度分析了微生物处理砂土固结不排水循环剪切作用下的试样力学特性。将模拟结果与室内试验相对比,探究微生物气泡对砂土抗液化性能的影响及其作用机理。

1 CFD-DEM流-固耦合计算基本原理

目前,流固耦合的方法主要有3种。笔者采用的CFD-DEM耦合方法考虑了不排水试验中的边界移动、含气泡流体的可压缩性、以及土体试样在循环荷载作用下加卸载刚度的不同。此方法适用于微生物气泡引起流体压缩系数变化的不排水循环三轴模拟。

1.1 颗粒运动方程

颗粒运动方程与经典DEM方法一致,均遵循牛顿第二定律:

(1)

(2)

1.2 流体运动控制方程

(1) 连续性方程

根据质量守恒定律推导得出的流动控制方程称为连续方程。单位时间内,通过体积单元的质量净通量等于体积元内质量的变化量,即

(3)

式中:εf为一个流体计算单元中的孔隙率;ρf为流体密度,▽为拉普拉斯算子;uf和ug分别为流体运动速度和颗粒运动速度。

(2) 动量方程

动量方程遵循牛顿第二定律,其代表的物理意义是作用于控制体上力的总和等于控制体的质量乘以控制体运动时的加速度。其可由公式(4)表达:

-εf▽P+▽·(εfμf▽uf)

(4)

式中:P为流体压力;μf为流体动力粘滞系数。

(3) 能量方程

将热力学第一原理应用于流体模型的控制方程称为能量方程。其代表的物理意义是控制体内的能量变化率等于流入控制体内的净热流量加上体积力和表面力对控制体做功的功率。可由公式(5)表达:

▽·(Puf)+εfρfFuf

(5)

(4) 液体状态方程

在不考虑温度和压力耦合的情况下,能量方程是满足的。而只有动量方程和连续性方程是无法闭合求解的。为了形成一个封闭的方程组,本文引入了流体的状态方程。一般情况下,水的可压缩性很低,可以忽略不计。而对于含气泡的非饱和液体,其可压缩性就会随着气泡含量的增加而提高。根据液体增压前后质量保持不变的性质,得到如下液体压缩系数的表达式:

(6)

本文通过液体的压缩系数建立起CFD的状态方程,如下:

ρf-ρf0=ρfC(p-p0)

(7)

式中:ρf0为流体初始密度;p0为流体初始压力。

微视频也是随着网络发展而产生的。我们都知道现在社会的生活节奏是相当的快，人们一天几乎没有很好的时间去享受自己的空闲时间，也没有时间去利用自己的娱乐时间。长的视频需要花大量的时间去看，而且看过之后心里也总是在想着剧情，容易扰乱自己的思想，不能专心工作。也正是因为这样的情况，我们的段视频应运而生，它不需要花费很多的时间，简单地几分钟就可以浏览完[4]。那我们的商务英语翻译课程就可以制作成简短的视频，如果内容长，我们就可以做成相关系列的东西，便于人们利用闲碎的时间去学习，只需简短的几分钟，便可以轻松搞定。再有闲暇时间还可以拿出来重新，温习一下，让知识得到巩固的效果。这就是利用做简单的微视频的好处。

1.3 CFD网格移动与网格内孔隙比计算

在模拟过程中,试样颗粒被六道刚性墙包围,上下墙的速度保持不变,侧墙的速度根据伺服系统得到。因此,CFD边界条件采用速度边界,即靠近边界的液体速度与墙的速度一致。加载过程中,由于试样形态发生变化,CFD的边界网格需要与墙同时移动,内部网格则采用简单的均匀变形方式移动,以x方向为例,内部网格的运动控制方程为:

(8)

1.4 耦合软件及流程

本文DEM计算通过PFC5.0进行,安装在Windows系统下,采用FISH语言和Python语言混合编程。而CFD计算通过OpenFOAM进行,安装在Linux系统下,采用Python语言编程。两者通过基于TCP的sockets编程实现CFD-DEM的数据交换。耦合计算流程示意图如图1所示。

图1 CFD-DEM耦合流程示意图Fig.1 Process schematic of CFD-DEM coupling

2 微生物处理过程的时效性

笔者采用文献[16]中的气泡生成方法,制备得到微生物处理2天、3天和5天后的含不同气泡含量的高饱和砂土试样。在将试样中的气泡体积与饱和度转换为数值模拟中的流体压缩系数时,不考虑水表面张力的影响,采用刘侃等[17]提出的含气孔隙流体压缩系数的表达式,如下所示:

(9)

式中:下标0代表初始状态;Va为气泡体积;Vw为水的体积;uw为孔隙水压;Pa为标准大气压;ua为气泡压力;Sr为饱和度;C0为饱和情况下的压缩系数。最终的参数计算结果如表1所示。

3 循环三轴CFD-DEM耦合模拟

为研究砂土抗液化能力,进行砂土固结不排水循环三轴室内试验时,制取试样后,施加围压使试样处于等应力状态,通过激振系统在竖向施加循环荷载,直至试样达到破坏状态。但室内试验无法获取试样的力链分布、配位数等微观数据,数值模拟可以直观试样破坏的整个过程,便于探究砂土液化的宏微观机理。

表1 数值模拟试样微生物处理过程时效性的参数

3.1 模拟步骤

含微生物气泡的高饱和砂土循环三轴模拟分成样、预压、固结、循环剪切四个过程。前三个过程采用DEM完成,最后的循环剪切过程采用CFD-DEM耦合。

(1) 成样。为与室内试验[18]对比,本文离散元模拟对象为渥太华50-70砂,级配采用Ottawa 50-70砂级配,初始孔隙比为0.73。模型微观参数如表2所列,离散元试样的制备采用分层欠压法[19],分10层制备,每层2 000颗粒。

表2 模型微观参数

(2) 预压。制得均匀试样后,固定四道侧墙,通过伺服系统对上下墙施加大小为12.5 kPa的应力进行预压。

(3) 固结。预压稳定后,通过伺服系统对试样进行等向固结。固结试样的平均有效主应力p=100 kPa,偏应力q=0。

(4) 循环剪切。此阶段中加入了CFD模块进行耦合模拟,CFD参数借鉴文献[20]。剪切时DEM时步被固定为1×10-6s/cycle,DEM与CFD的交换频率为25(即DEM计算25次,CFD计算1次),CFD时步为2.5×10-5s/cycle。液体的初始密度为1 000 kg/m3,压缩系数为4.5×10-7。网格数目为6×6×12 (长×宽×高),网格边界与DEM中六道墙体的坐标一致。循环加载前施加100 kPa反压。通过上下墙加载进行剪切,当偏应力达到应力幅值qcyc时,上下墙反向继续运动,当偏应力在相反方向达到应力幅值qcyc则再次反向。侧墙仍然采用伺服系统控制,伺服应力为总围压减去孔压后的有效应力。

3.2 模拟方案

为研究不同气泡含量的高饱和砂土的抗液化强度,选取了微生物处理2天、3天及5天的高饱和试样进行对比,详细工况见表3。其中,循环应力比为循环荷载作用下的最小主应力与最大主应力的比值,其表征平均应力与偏应力的联合作用。为方便描述,对每个试样进行了编号,例如,US2代表微生物处理2天的砂土试样,US2-1/2/3/4代表对微生物处理2天砂土试样的不同循环应力比循环三轴试验。

表3 微生物处理高饱和砂土循环三轴数值模拟工况列表

4 DEM模拟结果

4.1 力链分布

图2为US2-3试样在循环剪切过程中特征时刻的力链演化。图中力链越粗代表颗粒接触力越大,力链颜色同时也代表接触力的大小,红色代表的接触力最大,接下来依次是黄色、绿色,蓝色代表接触力最小。如图所示,随着加载的进行,试样颗粒接触力越来越小。试样经过等向固结后,内部力链分布较为均匀。试样在受压时,竖向受力为主导,内部力链主要沿竖向分布。试样在受拉状态时,侧向受力为主导,内部力链主要沿水平向分布。模拟结果表明:在加载过程中,孔隙水压力累积,承受了很大一部分力的作用,导致土骨架的有效传力减小。当试样中的有效应力为零时,便达到了初始液化状态。试样内部力链的分布方向随拉压状态的转换而转换,沿主要受力方向分布。

图2 试样US2-3加载过程中特征时刻力链分布演化Fig.2 The force chains distribution of sample US2-3 during loading

4.2 抗液化强度

微生物处理不同天数的高饱和砂土的抗液化强度曲线如图3(a)所示。图 3(b)(c)为室内试验[1,21]中高饱和砂土循环三轴试验抗液化强度曲线。图2(a)中,循环应力比CSR=0.25时,微生物处理2天的试样在2次振动后即达到液化状态,处理3天的试样达到液化状态需要振动38次,而处理5天的试样需要87次振动达到液化状态。模拟结果表明:含微生物气泡的高饱和砂土的抗液化强度比处理前的饱和净砂有所提升。且处理天数越多,气泡含量越多,试样的抗液化强度越高。模拟结果与室内试验相吻合。

图3 微生物处理高饱和砂土试样循环应力比-液化振次曲线Fig.3 Cyclic stress ratio-liquefaction vibration times curves of highly-saturated sand treated by microorganisms

4.3 孔压比

孔压比是指在循环应力作用下试样的孔隙水压力增量与侧向有效固结应力的比值。图4为循环应力比CSR=0.25时不同气泡含量的试样孔压比发展曲线,同时给出饱和试样在CSR=0.15下的孔压比。如图所示,微生处理后的试样的CSR比饱和试样大,含气泡的试样孔压比曲线发展更为缓慢。饱和试样和微生物处理2天后的试样的孔压比在加载后都迅速上升至1,达到初始液化状态。试样经过微生物处理后的第3天和第5天都没有达到初始液化状态,其最终的孔压比分别在0.5和0.4附近波动,且达到累积应变破坏所需的加载时间大幅度增加。模拟结果表明:随着微生物处理的天数增加,试样的饱和度不断减小,孔压比的累积随之变弱。

图4 不同气泡含量试样孔压比发展曲线Fig.4 Pore pressure ratio development curves of samples containing different bubbles

4.4 轴向应变

图5所示为循环应力比为CSR=0.25时的不同微生物气泡含量的试样轴向应变发展曲线,图中还给出了饱和试样在CSR=0.15下的轴向应变曲线。从图中可以看出,微生物处理2天后的试样的轴向应变发展曲线与饱和试样较类似,试样达到初始液化状态时的最终应变都在2%以内。微生物处理3天和5天后的试样轴向应变不断累积增长,直至达到5%时发生累积应变破坏。模拟结果表明:随着微生物处理时间的增加,试样的饱和度逐渐降低,应变累积随之变慢。

图5 不同气泡含量试样轴向应变发展曲线Fig.5 Axial strain development curves of samples with different bubble contents

4.5 力学配位数

力学配位数是指单个颗粒所含有接触力的数目,通常要求每个接触的法向接触力大于零。其曲线可以表征加载过程中试样的微观变化。图6为循环应力比CSR=0.25时不同微生物气泡含量的试样力学配位数发展曲线,同时给出了饱和试样在CSR=0.15下的力学配位数曲线。从图中可以看出,加载前各试样的力学配位数都是4。微生物处理2天后的试样与饱和试样的力学配位数发展曲线较为相似,达到初始液化状态时力学配位数都下降至3.5附近。而微生物处理3天和5天后的试样不发生初始液化破坏,它们的力学配位数发展的整体规律为先上升后下降。结果表明:经微生物处理后,试样的力学配位数发展曲线发生了变化,且当气泡含量达到一定值后,试样发生累积应变破坏。