考虑动态流体网格的颗粒-流体耦合算法

doi:10.16183/j.cnki.jsjtu.2020.317

考虑动态流体网格的颗粒-流体耦合算法

何金辉, 李明广^,, 陈锦剑, 夏小和

上海交通大学土木工程系, 上海 200240

Particle-Fluid Coupling Algorithm Considering Dynamic Fluid Mesh

HE Jinhui, LI Mingguang^,, CHEN Jinjian, XIA Xiaohe

Department of Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: 李明广,男,副研究员,电话(Tel.):18801900699;E-mail:lmg20066028@sjtu.edu.cn.

责任编辑: 陈晓燕

收稿日期: 2020-09-30

基金资助:

国家自然科学基金面上项目(41977216)
国家自然科学基金面上项目(51978399)
上海市科委项目(BI0100054)

Received: 2020-09-30

作者简介 About authors

何金辉(1995-),男,河南省洛阳市人,硕士生,从事岩土流固耦合方面研究

摘要

传统颗粒-流体耦合算法难以考虑流体动态边界问题,在模拟大变形案例时因流固边界不匹配易造成计算误差,影响结果的准确性.针对该问题,引入流体网格动态更新方法,推导动网格下的达西渗流方程和颗粒-流体相互作用方程,在离散元商业软件PFC^2D基础上,开发了考虑流体动网格的颗粒-流体耦合算法.将该算法用于模拟饱和土不排水剪切双轴试验,通过与常体积法的结果对比,验证了该方法的有效性.最后,采用该算法模拟了不同围压下的不排水双轴试验,计算结果的规律与室内试验具有较好的一致性.该算法考虑了流体动网格的问题,在模拟大变形下的三轴压缩试验、一维固结试验等案例时可获得相适应的流固边界,有助于提高模拟精度,因此可为类似研究提供参考.

关键词： 流固耦合; 动网格; 常体积法; 双轴试验

Abstract

It is generally difficult to consider the fluid dynamic boundary problem in the traditional particle-fluid coupling algorithm, causing calculation errors owing to the mismatching of the fluid-solid boundary and affecting the accuracy of the results in the modeling of large deformation issues. In view of this problem, the dynamic updating method of fluid mesh is introduced and Darcy’s seepage equation and the particle-fluid interaction equation in dynamic mesh are derived. Based on the discrete element commercial software PFC^2D, the particle-fluid coupling algorithm considering dynamic fluid mesh is developed. The proposed algorithm is applied to simulate the undrained shear biaxial test of saturated soil. The comparison result with the constant volume method verifies the effectiveness of the developed algorithm. Finally, the algorithm is used to simulate the undrained biaxial tests at different confining pressures. The law of the calculated results agrees well with that of the laboratory tests. By considering the problem of fluid dynamic boundary, the developed algorithm can obtain the fluid-solid boundary matching in the simulation of triaxial compression test and one-dimensional consolidation test or other cases in large deformations, which can help to improve the simulation accuracy and offer a theoretical reference for similar studies.

Keywords： fluid-solid coupling; dynamic mesh; constant volume method; biaxial test

PDF (1573KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

何金辉, 李明广, 陈锦剑, 夏小和. 考虑动态流体网格的颗粒-流体耦合算法[J]. 上海交通大学学报（自然版）, 2021, 55(6): 645-651 doi:10.16183/j.cnki.jsjtu.2020.317

HE Jinhui, LI Mingguang, CHEN Jinjian, XIA Xiaohe. Particle-Fluid Coupling Algorithm Considering Dynamic Fluid Mesh[J]. Journal of shanghai Jiaotong University, 2021, 55(6): 645-651 doi:10.16183/j.cnki.jsjtu.2020.317

由于天然土具有散粒性及含水性,采用微观流固耦合方法进行岩土材料不排水特性分析已成为越来越多学者的选择.目前,微观耦合法主要有基于离散-连续耦合理论的管域法和粗网格法^[1].管域法按照初始颗粒的分布划分流体域,不同流体域通过管道进行流体交换^[2].该方法简单直观,易于实现,近年来得到了广泛应用.Zhang等^[3]采用管域法模拟了黄土中裂隙注浆及其发展过程,并同模型试验结果进行了对比分析,发现最终都会形成“Y”形注浆裂缝.Zeng等^[4]通过管域法研究了横向各向同性岩体中水压致裂的机理,并对注水速率、地应力比和倾角等参数进行了敏感性分析.Wang等^[5]在管域法基础上,根据流动守恒原则重新建立了流体压力更新准则,并通过稳态饱和介质渗流验证了方法的可靠性.尽管管域法可简化耦合模拟过程并获得较为准确的结果,但流域结构的分布只与颗粒初始位置有关,而不会根据颗粒位移进行更新,因此仅适用于小变形的情形.

粗网格法对颗粒群划分初始流体网格,通过求解N-S方程来更新流体信息并计算流体对颗粒的作用力^[6].因计算方便,许多学者采用该方法研究宏观现象机理,推动了相关理论的发展.如Zeghal等^[7]采用粗网格法研究了饱和无黏性沉积土在动态激励作用下的液化机理,揭示了若干沉积质液化的微观力学机理和响应模式.Tao等^[8]通过粗网格法对管涌的微观力学机制进行了研究,弥补了相关理论的空白.王学红^[6] 通过粗网格法研究了振动桩下沉的微观力学机理及环境响应规律,对完善相关理论、指导现场施工具有重要意义.粗网格法中的流体网格始终是固定的,不会随颗粒群边界的变化而移动,因此在大变形案例模拟中具有一定局限性.

随着微观耦合方法应用场景的扩大化、复杂化,考虑大变形下流固边界的匹配性具有重要的研究价值和现实意义.目前,已有学者基于传统耦合法进行了改进,考虑了流固边界相适应的因素,并进行了验证.如Zhang等^[9]通过把动网格节点的速度引入N-S方程,从而考虑动边界的影响;Zhang等^[10]在研究潜蚀问题时,根据大颗粒初始分布建立四面体流体网格,后续模拟过程中结合颗粒位置进行网格重新划分,从而实现流体网格的实时更新.当前动网格方法虽解决了部分应用场景下流固边界不匹配的问题,但在推广过程中仍存在一定局限性.比如,N-S方程迭代运算效率低,方程因简化而未考虑部分项(比如流体-颗粒相互作用力),四面体流体网格算法空间搜索耗时长,应用场景受限制(需大小颗粒共存).因此,亟需一种原理简单、易于实现、适用范围广的考虑动边界的微观双向耦合算法.

本文在离散元商业软件PFC^2D的基础上,引入流体网格动态更新方法,推导动网格下的达西渗流方程和颗粒-流体相互作用方程,并通过Python语言和C语言混编的方法将算法嵌入软件,开发了可考虑流体动边界的微观耦合模块.将该算法用于模拟饱和土不排水剪切双轴试验,并通过和常体积法的结果对比验证了算法的有效性.该算法扩大了微观耦合法的适用范围,在大变形算例的模拟中有助于提高仿真精度,同时通过技术优化,在一定程度上解决了运算效率低的问题,因此对微观尺度上的大变形流固耦合模拟具有一定的参考价值.

1 颗粒-流体耦合算法原理及程序实现

1.1 颗粒场分析

在不考虑流场作用下,颗粒一般受到体力、阻尼力及粒间相互作用力的影响.本文考虑了流体相和颗粒相之间的双向耦合,因此,颗粒还受到流体对其施加的作用力.在流场中,颗粒一般受到流体静态作用力和流体动态作用力的影响^[11].本文采用该耦合算法模拟双轴试验,因加载过程中流体和颗粒速度均较慢,因而可以忽略流体动态力的影响.根据牛顿第二定律,颗粒受力情况可以表示为

(1)

m_{p} a_{p} = \sum_{c} f_{c} + f_{static} + f_{d}

(2)

I_{p} α_{p} = \sum_{c} r_{c} \times f_{c} + M_{d}

式中:m_p,I_p 分别为颗粒p的质量和转动惯量;a_p,α_p 分别为颗粒p的加速度和角加速度; $\sum_{c} f_{c}$ 为颗粒p所受接触力的合力;f_static 为颗粒受到的流体静态作用力;r_c为由颗粒p中心指向其各个接触点的向量;f_d 和M_d为由阻尼作用引起的力和力矩.

1.2 流体场分析

从流体角度进行分析时,主要考虑体积应变引起孔压增量和不同孔隙之间发生渗流作用两个因素.每经历一个耦合间隔时,对流体网格位置及流体信息进行更新,然后返回对颗粒的作用力矩阵.本文参考文献[9]中的动网格实现方法,对网格的速度和位置进行更新,如下式:

(3)

u_{i}^{m} = u_{i}^{1} + \frac{m - 1}{a - 1} (u_{i}^{a} - u_{i}^{1})

(4)

x_{j}^{n} = x_{j}^{1} + \frac{n - 1}{b - 1} (x_{j}^{b} - x_{j}^{1})

式中: i,j=1,2,分别代表x,y方向;u_i^m 为第i个方向上第m个节点沿i方向的速度分量;u_i¹,u_i^a 分别为第i个方向起始边界节点和终止边界节点沿i方向的速度分量;a为第i个方向上节点总数目;x_jⁿ 为第j个方向上第n个节点单位耦合间隔位移增量沿j方向的分量;x_j¹,x_j^b分别为第j个方向起始边界节点和终止边界节点单位耦合间隔位移增量沿j方向的分量;b为第j个方向上节点总数目.

对于更新之后的网格及此时的颗粒分布状态,根据文献[12]中孔隙结构变化引起孔压累积的原理方法进行孔压增量计算.首先对各个网格有:

(5)

\bar{v} = \frac{1}{N} \overset{N}{\sum_{k = 1}} v^{k}

(6)

\bar{X} = \frac{1}{N} \overset{N}{\sum_{k = 1}} X^{k}

式中: $\bar{v}$ , $\bar{X}$ 分别为某流体网格内所有颗粒的平均速度及平均位置;N为网格内颗粒数目;v^k为第k个颗粒的速度向量;X^k为第k个颗粒的位置向量.对此N个颗粒,可求颗粒k的相对速度 $v_{rel}^{k}$ 和相对位置 $X_{rel}^{k}$ :

(7)

v_{rel}^{k} = v^{k} - \bar{v}

(8)

X_{rel}^{k} = X^{k} - \bar{X}

根据弹性力学理论,预估相对速度 $v_{pre}^{k}$ 可以由速度梯度张量G和相对位置 $X_{rel}^{k}$ 求得:

(9)

\begin{matrix} v_{pre}^{k} = G X_{rel}^{k} \end{matrix}

根据式(7)和(9),预测相对速度v_pre和实际相对速度v_rel间的误差δ可以表示为

(10)

\begin{matrix} δ = \overset{N}{\sum_{k = 1}} (v_{pre}^{k} - v_{rel}^{k})^{2} \end{matrix}

为使误差 $δ 最小, δ 应满足 \frac{\partial δ}{\partial G} = 0,$ 由此可得:

(11)

\begin{matrix} \overset{N}{\sum_{p = 1}} G_{rs} x_{rel, s}^{p} x_{rel, t}^{p} = \overset{N}{\sum_{p = 1}} v_{rel, r}^{p} x_{rel, t}^{p} \end{matrix}

式中: G_rs 为速度梯度张量矩阵G的第r行第s列元素;x_rel,s^p,x_rel,t^p 分别为颗粒p相对位置的第s个和第t个分量;v_rel,s^p 为颗粒p相对速度的第s个分量.根据式(11)可求得G,因此耦合间隔时间Δt内的体积应变ε_v计算式为

(12)

ε_{v} = {\overset{\cdot}{ε}}_{v} Δt = \frac{1}{2} (G_{ww} + G'_{ww}) Δ t = G_{ww} Δt

式中:G_ww、G'_ww为矩阵G、G'主对角元素的和,G'=G^T.因此,孔隙水压力增量Δp可以表示为

(13)

\begin{matrix} Δ p = B_{f} ε_{v} \end{matrix}

式中:B_f为流体体积模量.在空间流体网格之间,由于流体压力梯度的存在,会发生流体的渗流作用.为了避免繁重的计算任务,本文不考虑采用N-S方程进行模拟,而是采用达西渗流定律进行计算.如图1所示,每个流体单元在压差作用下,都会与周围其他单元发生流体交换(I,J表示流体网格在二维空间内的位置坐标,取1,2,…).任意流体单元因压力梯度而与相邻单元发生流体交换(流入/流出)的速度可以表示为^[13]

(14)

\begin{matrix} v_{gh} = \frac{2 K K_{h}}{K + K_{h}} \frac{p - p_{h}}{(U_{g} - U_{gh}) γ_{w}} \end{matrix}

式中: g=1,2表示二维空间某方向;h=1,2表示g方向坐标轴的负向或正向;v_gh 为某网格g方向上h侧的边界渗流速度;K和K_h 分别为该网格和g方向h侧相邻网格的渗透系数;p和p_h 分别为上述相邻两网格的孔隙水压力;U_g 和U_gh 分别为相邻两网格中心位置g方向上的坐标分量;γ_w为流体的容重.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 某网格与相邻网格间渗流作用示意图

Fig.1 Schematic diagram of seepage action between a certain mesh and adjacent meshes

因为网格是时刻变化的,所以需要根据网格的速度对流体交换速度进行修正,获得真实的渗流流速.结合式(3),(14),可得边界修正渗流速度为

(15)

\begin{matrix} v'_{gh} = v_{gh} + u_{g}^{p_{0}} \end{matrix}

式中: $u_{g}^{p_{0}}$ 为g方向上该边界所对应的第p₀ 个节点的速度.根据修正后的渗流速度,可计算出单位耦合间隔内的总的渗流量Q (以流入为正)为

(16)

\begin{matrix} Q = \overset{2}{\sum_{g = 1}} \overset{2}{\sum_{h = 1}} v'_{gh} S_{gh} Δ t \end{matrix}

式中: S_gh 为某网格g方向上h侧边界的面积.渗流引起的孔压改变量p'_darcy 可以用渗流量Q计算,如下式:

(17)

\begin{matrix} p'_{darcy} = B_{f} \frac{Q}{V_{cell}} \end{matrix}

式中:V_cell为网格容积.由上可得新的孔隙水压力p_new为

(18)

\begin{matrix} p_{new} = p_{old} + Δp + p'_{darcy} \end{matrix}

式中:p_old为该流体网格上一时刻的总孔压.

1.3 两相耦合场分析

如1.1节所述,本文只考虑流体对颗粒的静态作用力.如图2所示,当颗粒完全浸入网格时,颗粒处于平衡状态;而当颗粒处于边界位置时,颗粒处于不平衡状态,将受到流体对其的作用力.由此可知,流体网格内某不平衡状态颗粒所受的静态作用力为^[14]

(19)

\begin{matrix} f_{static} = p_{new} A_{cr} n \end{matrix}

式中:A_cr为流体网格内流场对不平衡状态颗粒的等效作用面积;n为颗粒所受流体作用力方向上的单位向量.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 流体网格中颗粒位置及受力情况示意图

Fig.2 Schematic diagram of positions of particles and their stress states in a fluid mesh

1.4 颗粒-流体耦合算法程序实现

该耦合算法的实现思路为:① 进行力学计算,在未达到程序终止条件的前提下,当该力学进程时间Δt_cur等于设定的耦合间隔时间Δt时,根据颗粒边界的速度和位移,对网格位置及节点速度进行更新.② 在更新后的流体网格下,求解因颗粒运动引起的体积应变,并更新孔隙水压力.③ 根据动网格下的修正达西渗流方程对流场进行更新.④ 判断流场内的颗粒浸润状态,计算流体对颗粒的作用力,并以外力形式施加于各个颗粒,然后进入下一个循环.具体流程如图3所示.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 耦合流程示意图

Fig.3 Schematic diagram of coupling process

2 颗粒-流体耦合算法验证

2.1 试样制备

首先在宽度为10 cm、高度为20 cm的容器内,按照初始孔隙率0.16生成试样,生成的颗粒数目为 1004.试样的颗粒级配曲线示意图如图4所示.图中:d_p为颗粒直径,w为试样中直径小于某数值的颗粒所占的质量分数.采用线性接触本构模拟密实粒状土不排水剪切特性,本构微观参数取值如表1所示.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 颗粒级配曲线图

Fig.4 Diagram of grain size distribution

表1 颗粒参数

Tab.1 Particle parameters

颗粒参数	取值	颗粒参数	取值
颗粒密度/(kg·m^-3)	2500	切向刚度/(MN·m^-1)	100
最小颗粒直径/mm	3.6	摩擦因数	0.5
最大颗粒直径/mm	5.58	渗透性系数/(μm·s^-1)	10
法向刚度/(MN·m^-1)	100	试样初始孔隙率	0.16

新窗口打开| 下载CSV

试样生成后,对材料施加100 kPa的有效应力,从而实现试样的等向固结.之后,保持100 kPa的围压,在耦合条件下进行不排水双轴试验,其中流体模量为2.0 GPa,密度为 1000 kg/m³.双轴试验模型示意图如图5所示,加载过程中对上下刚性墙体施加5 mm/s的轴向加载速度,同时通过伺服机制保持侧向围压恒定,记录加载过程中偏应力、超孔隙水压力及轴向应变等用于后续的对比和分析.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 双轴试验模型示意图

Fig.5 Schematic diagram of biaxial test model

2.2 常体积法原理

虽然常体积法不能考虑土水之间相互作用机制,但因其计算效率高、模拟结果基本能反映不排水剪切规律,常被用于耦合算法的验证.

当采用常体积法时,通过控制试样体积不变来模拟材料的不排水行为^[15].因此,当轴向加载速度为v_a时,计算得到的径向速度v_r应保证试样在任意时间间隔Δτ内体积保持不变,即体积应变ε_v为0,具体为

(20)

\begin{matrix} ε_{v} = \frac{2 v_{a} Δτ W_{x} - 2 v_{r} Δτ W_{y}}{W_{x} W_{y}} = 0 \end{matrix}

式中: $W_{x}, W_{y} 分别为容器 x, y$ 方向的长度.由式(20)可以得出:

(21)

\begin{matrix} v_{r} = \frac{W_{x}}{W_{y}} v_{a} \end{matrix}

不排水剪切试验中,侧向总应力(围压) σ通过伺服保持为定值,即σ=u+σ'=C(u和σ'分别为孔隙水压力和有效应力,C表示某常数),因此u可以表示为u=σ-σ'.常体积法模拟不排水剪切试验时,初始时刻u=0,此时σ=σ'₀ (σ'₀ 为初始固结围压),加载阶段可获得侧墙的有效应力σ'₂,因此可直接求得每一时刻的u=σ'₀-σ'₂.

2.3 结果对比验证

图6所示为常体积法和耦合法的超孔隙水压力对比图,图中:ε_a为轴向应变,p_exc为侧墙附近产生的超孔隙水压力.可以看出,两种方法下的孔压规律较为一致,数值比较接近.初始时刻均生成正孔压,然后试样发生体积剪胀,孔压逐渐消散,变为负孔压.该趋势与较密实砂土的不排水剪切行为一致.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 耦合法和常体积法下的超孔隙水压力对比

Fig.6 Comparisons of excess pore water pressure in coupling method and constant volume method

图7所示为两种方法下的偏应力曲线,图中q为偏应力.由图可得,轴向应变为4%以内时,两种方法下的偏应力数值基本没有差别.继续加载,轴向应变超过4%时,二者数据虽有差异但偏差不大,均表现为随着加载过程而偏应力不断增大的趋势,同时轴向应变在10%附近,两者趋于相近.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 耦合法和常体积法下的偏应力结果对比

Fig.7 Comparisons of deviatoric stress results in coupling method and constant volume method

事实上,常体积法理论中流体是不可压缩的,也未考虑土水的相互作用,通过间接方式考虑流体信息的计算和更新,这种简化使得其模拟结果与实验数据不能高度吻合,但其结果显示的规律性同实验较为一致.从图6、7可以看出,本文提出的微观耦合法同常体积法数值比较接近,模拟规律一致,都体现出密实砂土在不排水剪切时的“硬化”趋势.同时该耦合法结合了流固耦合及动网格技术,模拟过程更为贴合实际,因此可为类似研究提供参考.

3 不同围压影响下的双轴试验结果分析

为了进一步验证方法的有效性,采用耦合法模拟了围压分别为100,200,300 kPa时的饱和土不排水剪切双轴试验,并对偏应力、孔隙水压力及有效应力比的变化情况进行了对比分析.

图8所示为不同围压下的偏应力变化情况.从图8可以看出,当围压增大时,相同应变时偏应力更大.原因为大的围压对土体的侧向膨胀限制作用更加明显,颗粒间的摩擦力更大,表现为土体模量较大,因此剪切强度也就越大.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 不同围压下偏应力结果对比

Fig.8 Comparison of deviatoric stress results at different confining pressures

图9所示为不同围压下的超孔压发展情况.可以看出,不同围压下试样的超孔压变化趋势均为先升高后下降,表明3组试样在应变较小时孔隙受到压缩,产生正的孔隙水压力,随着加载的进行,试样发生了剪胀,孔压逐渐降低,出现负值.对比3组试样的孔压曲线可知,当围压增大时,正孔压幅值增大,峰值点延后,负孔压值较小.因为随着围压的增大,侧墙对试样剪胀的抑制作用更明显,所以试样受压缩程度较高,因而正孔压的累积作用也就较强.

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 不同围压下的超孔隙水压力结果对比

Fig.9 Comparisons of excess pore water pressure results at different confining pressures

图10所示为不同围压下有效应力比η'的变化趋势图.三者有效应力比峰值均出现在轴向应变为1%左右,之后应力比逐渐降低.在1%附近时,100,200及300 kPa围压(后文默认为此顺序)对应的有效摩擦角分别为31.9°,31.3°及30.3°.3组试样在轴向应变9%附近的摩擦角差别相对较大,分别对应为18.0°,23.1°,22.3°.由图8可以看到,100 kPa围压的试样在轴向应变8%~9%时偏应力有一个突变,此时可能为该试样局部力链发生断裂,颗粒随机重新分布.而在加载终止时,三者的有效应力比趋于一致,其对应的摩擦角分别为20.3°,22.0° 及21.4°.综上,围压变化对摩擦角影响不大,这与文献[13]的研究结论一致.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 不同围压下有效应力比结果对比

Fig.10 Comparisons of effective stress ratio results at different confining pressures

4 结论

本文研究了一种考虑动态流体网格的颗粒-流体耦合算法,将算法嵌入离散元商业软件PFC^2D中,在大变形微观耦合算例模拟中可提高仿真精度.将该算法用于模拟饱和土不排水剪切双轴试验,并通过和常体积法的结果对比验证了算法的有效性.最后,采用该算法研究了不同围压影响下粒状土双轴剪切特性,模拟规律与室内试验结果较为吻合.主要结论如下:

(1) 在引入流体网格动态更新方法的前提下,改进并推导动网格下的达西渗流方程和颗粒-流体相互作用方程,在离散元商业软件PFC^2D基础上开发了考虑动态流体网格的颗粒-流体耦合模块.该算法原理简单,易于实现,弥补了传统微观耦合方法的不足,为大变形微观耦合模拟提供了参考,但在流场可视化、算法效率等方面仍需进一步的改进和研究;

(2) 将算法用于模拟饱和土不排水剪切双轴试验,并通过和常体积法的结果对比,验证了算法的有效性.结果显示,两种方法模拟得到的偏应力曲线和孔压曲线规律一致,数值较为接近.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

蒋明镜, 张望城.

一种考虑流体状态方程的土体CFD-DEM耦合数值方法

[J]. 岩土工程学报, 2014, 36(5):793-801.

[本文引用: 1]

JIANG

Mingjing

, ZHANG

Wangcheng

Coupled CFD-DEM method for soils incorporating equation of state for liquid

[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5):793-801.

[本文引用: 1]

[2]

刘广, 荣冠, 侯迪, 等.

颗粒-流体耦合算法与饱和土不排水剪切特性研究

[J]. 岩土力学, 2016, 37(1):210-218.

[本文引用: 1]

LIU

Guang

, RONG

Guan

, HOU

, et al.

Fluid-particle coupled model and a numerical investigation on undrained shear behavior of saturated soil

[J]. Rock and Soil Mechanics, 2016, 37(1):210-218.

[本文引用: 1]

[3]

ZHANG

Z L

, SHAO

Z S

, FANG

X B

, et al.

Research on the fracture grouting mechanism and PFC numerical simulation in loess

[J]. Advances in Materials Science and Engineering, 2018, 2018:1-7.

[本文引用: 1]

[4]

ZENG

Y W

, XIA

Numerical simulation of hydraulic fracturing in transversely isotropic rock masses based on PFC-2D

[J]. Journal of Vibroengineering, 2019, 21(4):833-847.

DOI:10.21595/jve.2018.19962 URL [本文引用: 1]

[5]

WANG

, DONG

, CHEN

Seepage simulation using pipe network flow model in a discrete element system

[J]. Computers and Geotechnics, 2017, 92:201-209.

DOI:10.1016/j.compgeo.2017.08.010 URL [本文引用: 1]

[6]

王学红

高频液压振动锤沉桩机理的颗粒离散元模拟

[D]. 福州: 福州大学, 2011.

[本文引用: 2]

WANG

Xuehong

Discrete element modelling of mechanisms of high frequency hydraulic vibratory pile driving

[D]. Fuzhou: Fuzhou University, 2011.

[本文引用: 2]

[7]

ZEGHAL

, EL

SHAMY U

A continuum-discrete hydromechanical analysis of granular deposit liquefaction

[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2004, 28(14):1361-1383.

DOI:10.1002/(ISSN)1096-9853 URL [本文引用: 1]

[8]

TAO

, TAO

Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method

[J]. Acta Geotechnica, 2017, 12(3):573-592.

DOI:10.1007/s11440-016-0516-y URL [本文引用: 1]

[9]

ZHANG

, JIANG

, THORNTON

A coupled CFD-DEM method with moving mesh for simulating undrained triaxial tests on granular soils

[J]. Granular Matter, 2020, 22(1):1-13.

DOI:10.1007/s10035-019-0969-4 URL [本文引用: 2]

[10]

ZHANG

, WANG

, LIU

, et al.

Modeling of fluid-particle interaction by coupling the discrete element method with a dynamic fluid mesh: Implications to suffusion in gap-graded soils

[J]. Computers and Geotechnics, 2020, 124:103617.

DOI:10.1016/j.compgeo.2020.103617 URL [本文引用: 1]

[11]

SHAFIPOUR

, SOROUSH

Fluid coupled-DEM modelling of undrained behavior of granular media

[J]. Computers and Geotechnics, 2008, 35(5):673-685.

DOI:10.1016/j.compgeo.2007.12.003 URL [本文引用: 1]

[12]

LIU

, RONG

, PENG

, et al.

Numerical simulation on undrained triaxial behavior of saturated soil by a fluid coupled-DEM model

[J]. Engineering Geology, 2015, 193:256-266.

DOI:10.1016/j.enggeo.2015.04.019 URL [本文引用: 1]

[13]

KHALILI

, MAHBOUBI

Discrete simulation and micromechanical analysis of two-dimensional sa-turated granular media

[J]. Particuology, 2014, 15:138-150.

DOI:10.1016/j.partic.2013.07.005 URL [本文引用: 2]

[14]

OKADA

, OCHIAI

Coupling pore-water pre-ssure with distinct element method and steady state strengths in numerical triaxial compression tests under undrained conditions

[J]. Landslides, 2007, 4(4):357-369.

DOI:10.1007/s10346-007-0092-1 URL [本文引用: 1]

[15]

KEISHING

, HANLEY

K J

Improving constant-volume simulations of undrained behaviour in DEM

[J]. Acta Geotechnica, 2020, 15(9):2545-2558.

DOI:10.1007/s11440-020-00949-1 URL [本文引用: 1]

一种考虑流体状态方程的土体CFD-DEM耦合数值方法

2014

... 由于天然土具有散粒性及含水性,采用微观流固耦合方法进行岩土材料不排水特性分析已成为越来越多学者的选择.目前,微观耦合法主要有基于离散-连续耦合理论的管域法和粗网格法^[1].管域法按照初始颗粒的分布划分流体域,不同流体域通过管道进行流体交换^[2].该方法简单直观,易于实现,近年来得到了广泛应用.Zhang等^[3]采用管域法模拟了黄土中裂隙注浆及其发展过程,并同模型试验结果进行了对比分析,发现最终都会形成“Y”形注浆裂缝.Zeng等^[4]通过管域法研究了横向各向同性岩体中水压致裂的机理,并对注水速率、地应力比和倾角等参数进行了敏感性分析.Wang等^[5]在管域法基础上,根据流动守恒原则重新建立了流体压力更新准则,并通过稳态饱和介质渗流验证了方法的可靠性.尽管管域法可简化耦合模拟过程并获得较为准确的结果,但流域结构的分布只与颗粒初始位置有关,而不会根据颗粒位移进行更新,因此仅适用于小变形的情形. ...

一种考虑流体状态方程的土体CFD-DEM耦合数值方法

2014

颗粒-流体耦合算法与饱和土不排水剪切特性研究

2016

颗粒-流体耦合算法与饱和土不排水剪切特性研究

2016

Research on the fracture grouting mechanism and PFC numerical simulation in loess

2018

Numerical simulation of hydraulic fracturing in transversely isotropic rock masses based on PFC-2D

2019

Seepage simulation using pipe network flow model in a discrete element system

2017

高频液压振动锤沉桩机理的颗粒离散元模拟

2011

... 粗网格法对颗粒群划分初始流体网格,通过求解N-S方程来更新流体信息并计算流体对颗粒的作用力^[6].因计算方便,许多学者采用该方法研究宏观现象机理,推动了相关理论的发展.如Zeghal等^[7]采用粗网格法研究了饱和无黏性沉积土在动态激励作用下的液化机理,揭示了若干沉积质液化的微观力学机理和响应模式.Tao等^[8]通过粗网格法对管涌的微观力学机制进行了研究,弥补了相关理论的空白.王学红^[6] 通过粗网格法研究了振动桩下沉的微观力学机理及环境响应规律,对完善相关理论、指导现场施工具有重要意义.粗网格法中的流体网格始终是固定的,不会随颗粒群边界的变化而移动,因此在大变形案例模拟中具有一定局限性. ...

... [6] 通过粗网格法研究了振动桩下沉的微观力学机理及环境响应规律,对完善相关理论、指导现场施工具有重要意义.粗网格法中的流体网格始终是固定的,不会随颗粒群边界的变化而移动,因此在大变形案例模拟中具有一定局限性. ...

高频液压振动锤沉桩机理的颗粒离散元模拟

2011

A continuum-discrete hydromechanical analysis of granular deposit liquefaction

2004

Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method

2017

A coupled CFD-DEM method with moving mesh for simulating undrained triaxial tests on granular soils

2020

... 随着微观耦合方法应用场景的扩大化、复杂化,考虑大变形下流固边界的匹配性具有重要的研究价值和现实意义.目前,已有学者基于传统耦合法进行了改进,考虑了流固边界相适应的因素,并进行了验证.如Zhang等^[9]通过把动网格节点的速度引入N-S方程,从而考虑动边界的影响;Zhang等^[10]在研究潜蚀问题时,根据大颗粒初始分布建立四面体流体网格,后续模拟过程中结合颗粒位置进行网格重新划分,从而实现流体网格的实时更新.当前动网格方法虽解决了部分应用场景下流固边界不匹配的问题,但在推广过程中仍存在一定局限性.比如,N-S方程迭代运算效率低,方程因简化而未考虑部分项(比如流体-颗粒相互作用力),四面体流体网格算法空间搜索耗时长,应用场景受限制(需大小颗粒共存).因此,亟需一种原理简单、易于实现、适用范围广的考虑动边界的微观双向耦合算法. ...

... 从流体角度进行分析时,主要考虑体积应变引起孔压增量和不同孔隙之间发生渗流作用两个因素.每经历一个耦合间隔时,对流体网格位置及流体信息进行更新,然后返回对颗粒的作用力矩阵.本文参考文献[9]中的动网格实现方法,对网格的速度和位置进行更新,如下式: ...

Modeling of fluid-particle interaction by coupling the discrete element method with a dynamic fluid mesh: Implications to suffusion in gap-graded soils

2020

Fluid coupled-DEM modelling of undrained behavior of granular media

2008

... 在不考虑流场作用下,颗粒一般受到体力、阻尼力及粒间相互作用力的影响.本文考虑了流体相和颗粒相之间的双向耦合,因此,颗粒还受到流体对其施加的作用力.在流场中,颗粒一般受到流体静态作用力和流体动态作用力的影响^[11].本文采用该耦合算法模拟双轴试验,因加载过程中流体和颗粒速度均较慢,因而可以忽略流体动态力的影响.根据牛顿第二定律,颗粒受力情况可以表示为 ...

Numerical simulation on undrained triaxial behavior of saturated soil by a fluid coupled-DEM model

2015

... 对于更新之后的网格及此时的颗粒分布状态,根据文献[12]中孔隙结构变化引起孔压累积的原理方法进行孔压增量计算.首先对各个网格有: ...

Discrete simulation and micromechanical analysis of two-dimensional sa-turated granular media

2014

... 式中:B_f为流体体积模量.在空间流体网格之间,由于流体压力梯度的存在,会发生流体的渗流作用.为了避免繁重的计算任务,本文不考虑采用N-S方程进行模拟,而是采用达西渗流定律进行计算.如图1所示,每个流体单元在压差作用下,都会与周围其他单元发生流体交换(I,J表示流体网格在二维空间内的位置坐标,取1,2,…).任意流体单元因压力梯度而与相邻单元发生流体交换(流入/流出)的速度可以表示为^[13] ...

... 图10所示为不同围压下有效应力比η'的变化趋势图.三者有效应力比峰值均出现在轴向应变为1%左右,之后应力比逐渐降低.在1%附近时,100,200及300 kPa围压(后文默认为此顺序)对应的有效摩擦角分别为31.9°,31.3°及30.3°.3组试样在轴向应变9%附近的摩擦角差别相对较大,分别对应为18.0°,23.1°,22.3°.由图8可以看到,100 kPa围压的试样在轴向应变8%~9%时偏应力有一个突变,此时可能为该试样局部力链发生断裂,颗粒随机重新分布.而在加载终止时,三者的有效应力比趋于一致,其对应的摩擦角分别为20.3°,22.0° 及21.4°.综上,围压变化对摩擦角影响不大,这与文献[13]的研究结论一致. ...

Coupling pore-water pre-ssure with distinct element method and steady state strengths in numerical triaxial compression tests under undrained conditions

2007

... 如1.1节所述,本文只考虑流体对颗粒的静态作用力.如图2所示,当颗粒完全浸入网格时,颗粒处于平衡状态;而当颗粒处于边界位置时,颗粒处于不平衡状态,将受到流体对其的作用力.由此可知,流体网格内某不平衡状态颗粒所受的静态作用力为^[14] ...

Improving constant-volume simulations of undrained behaviour in DEM

2020

... 当采用常体积法时,通过控制试样体积不变来模拟材料的不排水行为^[15].因此,当轴向加载速度为v_a时,计算得到的径向速度v_r应保证试样在任意时间间隔Δτ内体积保持不变,即体积应变ε_v为0,具体为 ...

〈

〉