任意解流流固耦合数值方法及在砂土渗流分析中应用

王胤; 陶奕辰; 程旷; 杨庆

doi:10.11779/CJGE202111015

任意解流流固耦合数值方法及在砂土渗流分析中应用 English Version

王胤^1,,
陶奕辰¹,
程旷¹,
杨庆¹

大连理工大学海岸和近海工程国家重点实验室,辽宁　大连　116024

基金项目:

国家自然科学基金项目 51879035

国家自然科学基金项目 51890912

详细信息

作者简介:
王胤(1982— ),男,教授,博士生导师,主要从事海洋土力学、海洋结构基础及相关岩土数值方法方面的教学和科研工作。E-mail：y.wang@dlut.edu.cn

中图分类号: TU43
计量
- 文章访问数: 359
- HTML全文浏览量: 39
- PDF下载量: 348
出版历程
- 收稿日期: 2020-12-24
- 网络出版日期: 2022-12-01
- 刊出日期: 2021-10-31

Arbitrary resolved-unresolved CFD-DEM coupling method and its application to seepage flow analysis in sandy soil

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China

摘要

摘要: 基于欧拉-拉格朗日连续体与非连续体耦合理论进行岩土工程流固耦合问题分析是一种较新颖和盛行的方法。针对该理论下的全解流（Fully-resolved）耦合与非解流（Un-resolved）耦合方法各自的缺陷,在已有的半解流（Semi-resolved）流固耦合数值方法（将全解流与非解流联合）基础上,通过引入修正的高斯权函数,建立了新的任意解流流固耦合方法（Arbitrary Resolved-Unresolved CFD-DEM coupling method）。该任意解流流固耦合方法能够较好地解决全解流方法中由于对粗颗粒周围流场精细化所带来的计算量过大问题；同时,成功地解决了流体网格内细颗粒较大时无法获得局部平均化变量问题；因此,该方法能够对具有一定粒径级配砂土土体的流固耦合问题开展模拟分析。通过室内砂土向上渗流试验,对所建立的任意解流流固耦合方法的准确性和有效性进行了验证；进一步地,采用该任意解流流固耦合模型,从细观层面上分析和研究了砂土渗流过程中水力梯度、土体变形随渗流速度变化规律。建立的任意解流流固耦合方法能够为岩土工程土体渗流问题的研究提供新的方法和手段。
- 流固耦合 /
- 砂土渗流 /
- 土颗粒 /
- 拖曳力 /
- 孔隙率
Abstract: The Euler-Lagrange coupling scheme based on the continuous and discrete theories has been becoming increasingly popular in numerical analysis of fluid-particle interaction. In this study, by introducing the modified Gaussian weighting function, a new arbitrary resolved-unresolved CFD-DEM coupling method (ARU CFD-DEM) is proposed based on the authors’ previous developed semi-resolved coupling approach by combing the fully-resolved and un-resolved coupling methods. This ARU CFD-DEM method is powerful to relieve the overload in computation due to refining the flow field around the coarse particles in the fully-resolving method. At the same time, it is also able to solve the difficulty in computing the local averaging variables when fine particles with large diameter exist in fluid grids. By doing so, the ARU CFD-DEM is able to simulate the fluid-particle interaction in sand mass which contains a wide range of particle diameters. By comparing with the results of upward seepage flow tests in sand, the accuracy and effectiveness of ARU CFD-DEM model is verified. Furthermore, the hydraulic gradient-flow velocity relationship and soil deformation-flow velocity relationship in the upward seepage flow are analyzed on the particle-scale by the ARU CFD-DEM. The proposed ARU CFD-DEM model can provide a new tool for investigating the fluid-particle interaction in the seepage flow in geotechnical engineering.
- CFD-DEM coupling /
- seepage flow /
- sand particles /
- drag force /
- void fraction

HTML全文

0. 引言

自然界土体由多相物质组成,对于饱和土而言,其由固体颗粒与粒间液体两相组成。在岩土工程问题分析中,不排水条件下土体的力学特性必然要考虑流-固耦合问题。考虑土颗粒与水相互作用的流固耦合分析方法很多,其中将固相视为离散颗粒,将液相考虑为连续体的欧拉-拉格朗日流固耦合方法在解决流体与固体间相互作用问题中比较适用,越来越多地被国内外学者应用到岩土力学问题分析中^[1-5]。

根据是否精确求解颗粒周围绕流流场,欧拉-拉格朗日耦合方法可分为全解流（Fully-resolved）和非解流（Un-resolved）两类。全解流方法精细求解每个固体颗粒周围绕流流场,并通过流体应力在颗粒表面积分来直接获得流体对于固体颗粒的作用力。根据流体建模方法不同,目前常用的全解流方法包括全解流CFD-DEM耦合方法和LBM-DEM耦合方法。由于全解流耦合方法对流体网格的分辨率要求较高,即为满足计算精度,要求流体网格尺寸至少小于固体颗粒直径的1/8～1/10^[6],因此导致该方法计算量较大。目前全解流耦合方法还只用于小尺度颗粒-流体耦合模拟,且主要用于固体颗粒数量较少的颗粒沉降、理想颗粒堆积体内渗流的岩土问题分析中^[7]。非解流方法无需精细求解固体颗粒周围的绕流流场,而是基于局部平均化理论,通过求解局部平均化的流体动量方程来获得流体局部区域内的平均化场量^[8]；并通过局部平均化场量和相间相互作用力模型,在局部平均化尺度上处理固体颗粒与流体间的双向耦合关系。目前在非解流耦合模拟中,常用的方法主要包括局部平均化的非解流CFD-DEM耦合方法和SPH-DEM非解流耦合方法。Tsuji等^[9]首次采用局部平均化非解流耦合方法模拟了向上渗流引起颗粒堆积体流化现象,随后该方法被广泛应用于颗粒-流体耦合数值分析中。与全解流耦合方法相比,非解流方法对流体网格分辨率要求较低；为达到平均化效应,要求每个流体网格内能够容纳3～4个固体颗粒^[3]。因此,非解流耦合方法在降低颗粒周围流场精细化程度同时,大大减少了计算量,显著提高计算效率。

针对以上两种耦合方法的缺陷,笔者先前已将全解流耦合和非解流耦合方法联合,提出了半解流（semi-resolved）CFD-DEM流固耦合方法^[10],该方法在模拟砂土渗流中能够同时兼顾流场求解精细化和耦合计算效率。半解流CFD-DEM耦合方法基本思想是,以一定尺寸划分流体网格,将土体颗粒分成粗颗粒（通常下粒径为流体网格尺寸8～10倍或以上）和细颗粒（粒径小于流体网格尺寸）两类。对于粗颗粒而言,流体网格满足所需的分辨率,采用全解流方法精细化求解周围流场,并直接积分求解流体与土体颗粒相互作用力；而对于细颗粒,流体网格达到平均化效应,采用非解流的局部平均化理论求解流场及其与细颗粒的相互作用。然而,半解流耦合方法适用于颗粒粒径差别较大的断级配砂土渗流模拟分析,粗颗粒与细颗粒粒径比α_s为8～10以上为佳。但对于存在连续级配的砂土,采用半解流方法求解比较困难。

鉴于如上问题,本文基于半解流CFD-DEM耦合方法,通过引入修正高斯权函数,进一步提出一种新的任意解流流固耦合方法（Arbitrary Resolved- Unresolved CFD-DEM coupling method, ARU CFD- DEM）,很好地解决了局部平均化非解流方法中流体网格内存在中间粒径颗粒（即粗颗粒与其粒径比α_s小于8～10,或细颗粒粒径大于流体网格尺寸）无法求解孔隙率的问题。本文首先介绍任意解流耦合方法的基本原理和控制方程；随后,阐述该方法的数值实现过程及程序框架；最后,通过砂土向上渗流问题模拟,验证所提出方法的有效性和准确性,并进一步从细观层面上分析砂土渗流过程中水力梯度、土体变形随渗流速度变化规律。

1. 任意解流流固耦合基本原理

1.1 流体控制方程

联合全解流的虚流体方法与非解流的局部平均化理论,在笔者已发展的半解流耦合算法^[10]基础上,可将通用的任意解流流体控制方程表示为

$\frac{\partial ε}{\partial t} + \nabla \cdot (ε u_{f}) = 0 (Ω (t))$ ,

(1)

$ρ_{f} [\frac{\partial ε u_{f}}{\partial t} + \nabla \cdot (ε u_{f} u_{f})] = - ε \nabla p + ε \nabla \cdot τ + F_{f p} + ε f_{F D} (i n Ω (t))$ ,

(2)

式中,ε为孔隙率,t为计算时间, $u_{f}$ 为流体速度,ρ_f为流体密度,p,τ分别为流体压力与黏性力,F_fp为流体与细颗粒间相互作用力,f_FD表征约束粗颗粒虚流体具有刚体运动而产生的附加相间作用力,可通过下式隐式求解获得

$u_{r} (x, t) = {U^{'}}_{c p} (y_{c p}, t) + {ω^{'}}_{c p} (y_{c p}, t) \times (x - y_{c p}) (\forall x \in Ω_{c p} (t)) ，$

(3)

式中,u_r为粗颗粒虚拟流体区域Ω_cp中任一点x的刚体运动速度, ${U^{'}}_{c p}$ , ${ω^{'}}_{c p}$ 分别为粗颗粒形心的平动和转动速度,y_cp表示粗颗粒中心点位置。式（2）中,当计算流体域中只存在细颗粒时,f_FD消失；而当流体域中只存在粗颗粒时,F_fp为零。

1.2 固体颗粒控制方程

固体颗粒运动的计算采用传统离散元方法即牛顿第二定律,除考虑颗粒间接触力和体积力外,还引入流体与颗粒间的相互作用力f_h：

$M_{i} \frac{d U_{i}}{d t} = \sum_{c} f^{c} + ρ_{i} V_{i} g + {f^{'}}_{h}$ ,

(4)

$J_{i} \cdot \frac{d ω_{i}}{d t} + ω_{i} \times (J_{i} \cdot ω_{i}) = \sum_{c} m^{c} + m_{h}$ 。

(5)

式中M_i,V_i, $ρ_{i}$ 分别为固体颗粒i的质量、体积和密度；J_i,U_i, $ω_{i}$ 分别为固体颗粒i的转动惯量、平动速度和转动速度；g为重力加速度； $f^{c}$ 为粒间接触力,包括法向力与切向力；m^c为固体颗粒所受到力矩,包括滚动力矩和由于接触力 $f^{c}$ 所引起的力矩。 $f^{c}$ 和 $m^{c}$ 通过Hertz-Mindline非线性粒间接触模型^[11]进行确定。m_h表示作用在粗颗粒上附件相间作用力所引起的转动力矩。 ${f^{'}}_{h}$ 表示流体作用在固体颗粒（包括细颗粒和粗颗粒）上的相间作用力。对于细颗粒, ${f^{'}}_{h}$ 包括液体的压差力即浮力 ${f^{'}}_{b}$ （=-ρr_fV_fpg）和拖曳力 ${f^{'}}_{d}$ ：

${f^{'}}_{d h} = - ρ_{f} V_{f p} g + (- \nabla p) V_{f p} + (\nabla \cdot τ) V_{f p} + {f^{'}}_{d}$ ,

(6)

式中,V_fp为细颗粒体积。对于固体土颗粒,近似考虑成具有形同体积的球形颗粒,等效半径为d_fp,其所受到的液体浮力可表示为

${f^{'}}_{b} = \frac{1}{6} π ρ_{f} d_{f p}^{3} g$ 。

(7)

拖曳力 ${f^{'}}_{d}$ 通过改进的Di Felice^[12]半经验公式确定：

${f^{'}}_{d} = {f^{'}}_{d}_{, 0} f (ε) = \frac{1}{8} C_{D} ρ_{f} π d_{f p}^{2} | u_{f} - {U^{'}}_{f p} | (u_{f} - {U^{'}}_{f p}) f (ε)$ ,

(8)

式中, ${U^{'}}_{f p}$ 为细颗粒运动速度, $f (ε)$ 为与孔隙率有关的影响参量,表达计算区域内周围固体颗粒对某一特定颗粒受流体拖曳力作用的影响效果,其可进一步表示为

$f (ε) = {\begin{cases} 1 (单颗粒体系) \\ ε^{- m} (多颗粒体系) \end{cases}$ ,

(9)

式中, $ε^{- m}$ 为颗粒间相互影响的孔隙率修正系数,一般为指数形式,Di Felice通过试验结果拟合确定指数m为

$m = 3.7 - 0.65 \exp [- \frac{{(1.5 - \lg_{10} R e_{p})}^{2}}{2}]$ ,

(10)

$R e_{p} = \frac{ε ρ_{f} d_{p} | U^{f} - U^{P} |}{μ_{f}}$ 。

(11)

本文对原始Di Felice^[12]模型进行了改进,通过引入三维颗粒形状系数Ψ,建立考虑砂土颗粒形状特性的单颗粒拖曳力系数C_D0,通过大量单颗粒在液体中沉降试验结果,拟合获得适用于任意形状颗粒的拖曳力系数为^[13-14]

$C_{D 0} = 0.945 \frac{C_{D, s p h}}{Ψ^{θ}} R e^{- 0.01}$ ,

(12)

式中,θ为雷诺数影响系数,可确定为0.641Re_p^0.153, $μ_{f}$ 为流体动力黏滞系数,单位为Pa·s。颗粒形状系数Ψ为

$Ψ = \frac{Φ}{X}$ ,

(13)

其中,Φ为颗粒比表面积,X为颗粒圆形度,其为颗粒最大投影面的周长与相等面积圆的周长之比。该形状系数考虑了颗粒三维形态,能够很好地描述不规则颗粒形状特征。Ψ值介于0～1,其值越大,颗粒形状越接近圆球体；反之,颗粒形状越不规则。颗粒形状系数的确定可进一步参考文献[13,14]。

对于粗颗粒, ${f^{'}}_{h}$ 和m_h可分别通过对作用在粗颗粒表面的流体应力进行体积积分获得

${f^{'}}_{h} = ∭_{CP} (- \nabla p_{d} + \nabla \cdot τ) d V - ρ_{f} g V_{c p}$ ,

(14)

$m_{h} = ∭_{CP} r \times (- \nabla p + \nabla τ) d V$ ,

(15)

式中,p_d为动水压力,V_cp为粗颗粒体积,r为相对于粗颗粒中心处的位置向量。

1.3 高斯分布函数及孔隙率求解

本文提出一种网格无关权函数—高斯权函数,与半解流流固耦合方法进行结合。传统高斯函数在统计学中用来描述正态分布。当作为局部平均化过程中的权函数时,高斯函数能够描述位于y点颗粒与x点处流体质点的相互作用随它们之间距离增加而逐渐减弱的效果。在半解流计算中需满足归一化条件和紧支性条件以保证局部平均化过程中的质量和动量守恒；同时,还要满足筛选条件和湮灭条件以避免细颗粒与虚流体之间的非物理相互作用。由于篇幅所限,以上4个条件可参考文献[10]。然而,传统高斯函数的支撑域为全空间,无法满足上述条件。因此,本文对高斯函数进行修正,在有限大小的球形空间域Ω_sd内对其进行截断以及满足以上条件处理（图1）,

$g_{G} (x - y) = {\begin{cases} δ (x - y) \begin{matrix} \end{matrix} (x \in Ω_{c p}) \\ \frac{H_{(Ω / Ω_{c p}) \cap Ω_{s d}} (y)}{N_{G} (x)} \exp (- \frac{| x - y |^{2}}{b_{w}^{2}}) (x \in Ω / Ω_{c p}) \end{cases} ，$

(16)

图 1 高斯权函数支撑域示意图

Figure 1. Diagram of supporting domain of Gaussian-based weighting function

下载: 全尺寸图片幻灯片

式中, $δ (x - y)$ 为狄拉克函数, $H_{(Ω / Ω_{c p}) \cap Ω_{s d}} (y)$ 为在 $(Ω / Ω_{c p}) \cap Ω_{s d} (x)$ （图1（b））上的单位阶跃函数,b_w为高斯函数的带宽,N_G(x)表达为

$N_{G} (x) = \underset{(Ω / Ω_{c p}) \cap Ω_{s d}}{∭} \exp (- \frac{| x - y |^{2}}{b_{w}^{2}}) d V_{y}$ 。

(17)

关于高斯函数的球形支撑域半径的选取,要保证截断后的高斯函数能够最大限度地保留高斯函数特征,高斯函数在支撑域内的累计权重足够大。这方面已作了相应的研究,选取支撑域半径为带宽的2倍以上即可满足高斯函数在支撑域内的累计权重接近1^[15]。

2. 耦合模型数值实现与验证

2.1 高斯分布函数离散与场变换

在任意解流流固耦合方法数值实现中,需要对流体网格一定体积V_m内欧拉场变量及颗粒相关的拉格朗日场变量进行转换。在细颗粒-流体耦合时,需基于修正的高斯权函数g_G(x_m-y)进行场变量转换。

首先,进行高斯权函数的空间离散化。通过将修正高斯权函数（式（16））的x直接由流体网格形心点x_m替换对高斯权函数进行空间离散。考虑网格形心相对于粗颗粒区域Ω_cp内、外部移动时,为避免g_G(x_m-y)发生突变,对δ(x_m-y)和 $\frac{H_{(Ω / Ω_{c p}) \cap Ω_{s d} (x_{m})} (y)}{N_{G} (x_{m})} \exp (- \frac{| x_{m} - y |^{2}}{b_{w}^{2}})$ 进行线性插值处理,得到网格V_m内修正高斯权函数离散化形式为

$\begin{matrix} g_{G} (x_{m} - y) = h (x_{m}, t) δ (x_{m} - y) + [1 - h (x_{m}, t)] \frac{H_{(Ω / Ω_{c p}) \cap Ω_{s d}} (y)}{N_{G} (x_{m})} \cdot \\ \exp (- \frac{| x_{m} - y |^{2}}{b_{w}^{2}}) = h (x_{m}, t) H_{m} (y) \frac{1}{V_{m}} + [1 - h (x_{m}, t)] \cdot \\ \frac{H_{(Ω / Ω_{c p}) \cap Ω_{s d}} (y)}{N_{G} (x_{m})} \exp (- \frac{| x_{m} - y |^{2}}{b_{w}^{2}}) \end{matrix}$ ,

(18)

式中,H_m(y)为流体网格内单位阶跃函数,h(x_m, t)为流体网格内粗颗粒体积分数,

$h (x_{m}, t) = ∭_{m} H_{Ω_{c p}} (y) d V_{y} = \frac{V_{c p o}}{V_{m}}$ ,

(19)

式中, $∭_{V_{m}} 为$ 在流体网格 $V_{m}$ 内进行体积积分, $V_{c p o}$ 为粗颗粒与流体网格重叠部分体积, $H_{Ω_{c p}}$ 为粗颗粒区域内的单位阶跃函数。

然后,基于高斯权函数进行拉格朗日-欧拉场变换。将与细颗粒相关的拉格朗日场变量,即细颗粒体积V_fp、拖曳力 ${f^{'}}_{d}$ (y_fp,t)和速度 ${\bar{U}}^{'}_{f p}$ (y_fp,t)变换为流体网格内的欧拉场变量,即孔隙率ε(x_m,t)、相间作用力F_fp(x_m,t)、和细颗粒局部平均速度 ${\bar{U}}^{'}_{f p} (x_{m}, t)$ ,

$\begin{matrix} ε (x_{m}, t) = 1 - [1 - h (x_{m}, t)] \cdot \\ \frac{\sum_{f p \in Ω \cap Ω_{s d}} [\exp (- \frac{| x_{m} - y_{f p} |^{2}}{b_{w}^{2}}) V_{f p}]}{\sum_{m^{'} \in Ω_{s d}} [h (x_{m}, t) \exp (- \frac{| x_{m} - y_{m^{'}} |^{2}}{b_{w}^{2}}) V_{m^{'}}]} \end{matrix}$ ,

(20)

$\begin{matrix} F_{f p} (x_{m}, t) = - [1 - h (x_{m}, t)] \cdot \\ \frac{\sum_{f p \in Ω \cap Ω_{s d}} [\exp (- \frac{| x_{m} - y_{f p} |^{2}}{b_{w}^{2}}) {f^{'}}_{d} (y_{f p}, t)]}{\sum_{m^{'} \in Ω_{s d}} [h (x_{m}, t) \exp (- \frac{| x_{m} - y_{m^{'}} |^{2}}{b_{w}^{2}}) V_{m^{'}}]} \end{matrix}$ ,

(21)

$\begin{matrix} {\bar{U}}^{'}_{f p} (x_{m}, t) = \frac{1 - h (x_{m}, t)}{1 - ε (x_{m}, t)} \cdot \\ \frac{\sum_{f p \in Ω \cap Ω_{s d}} [\exp (- \frac{| x_{m} - y_{f p} |^{2}}{b_{w}^{2}}) {U^{'}}_{f p} (y_{f p}, t) V_{f p}]}{\sum_{m^{'} \in Ω_{s d}} [h (x_{m}, t) \exp (- \frac{| x_{m} - y_{m^{'}} |^{2}}{b_{w}^{2}}) V_{m^{'}}]} \end{matrix}$ 。

(22)

2.2 流体运动方程数值求解

对任意解流耦合方法中的流体运动方程进行数值求解,将粗颗粒区域内虚流体流动约束为刚性运动,求解粗颗粒周围扰流流场。本研究中,为避免直接采用全解流算法而引起流体对静止粗颗粒发生穿透问题,在全解流算法基础上考虑细颗粒存在对求解粗颗粒周围绕流流场的影响,在笔者之前已发展的半解流算法中对该问题进行了处理^[10]。

在流体运动方程（式（1）,（2））的数值求解中,采用隐式欧拉格式进行时间离散,利用中心差分和迎风差分混合格式将对流项和扩散项进行空间离散^[16]。

$\frac{ε^{*} - ε^{n}}{Δ t} + \nabla \cdot (ε^{*} u_{f}^{*}) = 0$ ,

(23)

$\begin{array}{l} ρ_{f} [\frac{ε^{*} u_{f}^{*} - ε^{n} u_{f}^{n}}{Δ t} + \nabla \cdot (ε^{*} u_{f}^{*} u_{f}^{*})] \\ = - ε^{*} \nabla p^{*} + ε^{*} \nabla \cdot τ^{n} + F_{f p}^{*} + ε^{*} f_{F D}^{*} \end{array}$ ,

(24)

式中,上角标*为当前时间步有待确定的变量,上角标n为前一时间步已确定的变量。在进行流体运动方程求解时, $ε^{*}$ , $u_{f}^{*}$ 和 $F_{f p}^{*}$ 分别基于式（20）～（22）并根据当前时间步内更新的颗粒信息直接求得。 $f_{F D}^{*}$ 通过式（3）隐式计算获得。根据本文提出的半解流算法对离散后的流体运动方程（式（23）,（24））进行求解,概述以下4个主要步骤：

（1）预测阶段：忽略粗颗粒求解流体运动方程式,假定粗颗粒附加相间作用力 $f_{F D}^{*}$ 取为0,式（24）可简化为

$\begin{matrix} A_{P} {(u_{f})}_{P}^{*} + \sum_{N} A_{N} {(u_{f})}_{N}^{*} = - {(ε \nabla \frac{p_{f}}{ρ_{f}})}_{P}^{*} + \\ H [{(u_{f})}_{P}^{n}] + {(\frac{β_{f}}{ρ_{f}} u_{f p})}_{P}^{*} + {(\frac{ε u_{f}}{Δ t})}_{P}^{n} \end{matrix}$ ,

(25)

式中,

$H [{(u_{f})}_{P}^{n}] = ε_{P}^{*} \nabla \cdot {[\frac{μ_{f}}{ρ_{f}} (\nabla^{T} u_{f} - \frac{2}{3} t r (\nabla^{T} u_{f}) I)]}_{P}^{n}$ ,

(26)

其中,下标P和N分别表示该变量为当前流体网格和周围相邻网格中的变量,Σ_N表示遍历当前流体网格P的所有相邻网格N进行求和。A是由动量方程中的时间导数项、对流项、扩散项和相间作用力项离散而得的系数矩阵。为处理式（23）,（25）中的压力速度耦合,采用PISO（pressure implicit with splitting of operators）^[17]算法对其进行求解,得到临时速度场 $u_{f}^{*}$ 和压力场 $p_{}^{*}$ 。

（2）映射处理：为考虑粗颗粒对流体的作用,将粗颗粒区域内的虚流体流速显式地映射为粗颗粒刚体运动速度,

$u_{f}^{* *} = h u_{r} + (1 - h) u_{f}^{*}$ ,

(27)

式中, $u_{f}^{* *}$ 为映射后的流体流速,u_r为粗颗粒的刚体运动速度,通过式（3）计算获得,h为流体网格内粗颗粒体积分数（式（19））。

（3）连续性修正：经过映射后,需要对粗颗粒区域内、外部流体流速进行修正,使其均满足连续性方程（式（23））。当采用传统的全解流耦合算法修正时,由于粗颗粒周围细颗粒的存在会引起流速项粗颗粒内部发生不合理的绕流问题；特别地,当粗颗粒静止时,流体会对粗颗粒产生不合理的穿透现象。为解决该问题,本文提出半解流修正算法,将采用下式进行流体流速连续性修正：

$ε u_{f}^{n + 1} = ε u_{f}^{* *} - (1 - h) \nabla φ$ ,

(28)

式中,通过引入体积分数(1-h)来保证粗颗粒区域内的流速在修正前后始终维持为粗颗粒刚体运动速度,ϕ为标量场^[18]。

（4）压力重分布：对流体压力场进行重分布使其满足动量方程。由于在连续性修正中将 $u_{f}^{* *}$ 修正为 $u_{f}^{n + 1}$ ,打破了动量平衡,因而需要将压力场 $p_{}^{* *}$ 修正为 $p_{}^{n + 1}$ 以满足如下动量方程,

$\begin{array}{l} A_{P} {(u_{f})}_{P}^{n + 1} + \sum_{N} A_{N} {(u_{f})}_{N}^{n + 1} \\ = - {(ε \nabla \frac{p}{ρ_{f}})}_{P}^{n + 1} + H [{(u_{f})}_{P}^{n}] + {(\frac{β_{f}}{ρ_{f}} u_{f p})}_{P}^{*} + {(\frac{ε u_{f}}{Δ t})}_{P}^{n} + {(ε f_{C F})}_{P}^{n + 1} ， \end{array}$

(29)

式中,

$p_{}^{n + 1} = p_{}^{*} + Δ p$ ,

(30)

其中,Δp可通过其泊松方程求解获得

$\nabla \cdot (\frac{ε}{ρ_{f}} \nabla Δ p) = \nabla \cdot A_{P} {[(1 - h) \nabla φ]}_{P} + \nabla \cdot \sum_{N} A_{N} {[(1 - h) \nabla φ]}_{N} 。$

(31)

经过上述4个求解步骤,获得当前计算步最终的流体速度场 $u_{f}^{n + 1}$ 和压力场 $p_{}^{n + 1}$ ,数值求解过程见图2。

图 2 任意解流耦合数值求解流程图

Figure 2. Flow chart of numerical solution of ARU CFD-DEM coupling

下载: 全尺寸图片幻灯片

2.3 数值求解平台

基于开源离散元程序LIGGGHTS^[19]和开源CFD类库OpenFOAM^[20],本文自主开发了任意解流耦合求解器——ARUcfd-dem Solver。在该求解器中,基于离散元程序求解固体颗粒运动（式（4）,（5））；在OpenFOAM原有求解器icoFoam的基础上进行开发,编程实现任意解流耦合过程中离散元DEM与计算流体动力学CFD的信息交换,进而实现双向耦合计算。

2.4 基准渗流试验模拟

本文基于高斯函数进行场变量变换的算法,虽对于颗粒信息进行场变换及之后的修正步骤亦引入了额外的计算负担,但相比传统全解流方法,其计算效率仍得到了大幅提高。下面通过一个小规模基准性渗流试验模拟,验证本文所采用的任意解流方法相比较全解流方法的计算效率提升情况。基准渗流模拟所采用材料属性以及相关设置见表1,2。

表 1 基准渗流模拟所采用材料属性

Table 1. Material properties in the benchmark of seepage flow

颗粒属性						流体属性
颗粒密度ρ_s/(kg·m^-3)	滑动摩擦系数μ_r	杨氏模量E/Pa	回弹系数e	泊松比 $ν$	滚动摩擦系数μ_f	密度ρ_f/(kg·m^-3)	动力黏度μ_f/(kg·m^-1·s^-1)
2650	0.84	2×10⁷	0.9	0.2	0.26	1000	1×10^-3

下载: 导出CSV

| 显示表格

表 2 基准渗流模拟相关设置

Table 2. Computational settings in the benchmark of seepage flow

计算设置		任意解流设置		全解流设置
DEM时间步长Δt_DEM/s	CFD时间步长Δt_CFD/s	流体网格尺寸L_m/mm	高斯权函数带宽b_w/mm	流体网格尺寸L_m/mm
2×10^-5	2×10^-4	1.32, 0.8, 0.5	8.0	0.5

下载: 导出CSV

| 显示表格

建立向上渗流试验模拟模型,考虑多种不同粒径球形颗粒（细颗粒径为4～6 mm,粗颗粒经为13.2 mm）,最大粒径比为3.3,颗粒总数为782个。模拟区域为长方体,尺寸如图3所示。长方体底面为流体入口,顶面为流体出口,四周为墙体无滑动边界。采用任意解流方法时,对模拟区域流体采用立方体网格进行剖分,采用3种网格密度（即网格数量N_m为34200,156250,640000；粗颗粒直径与网格长度比分别为10,17,27）；采用全解流方法时,同样采用立方体网格对流体区域进行剖分,只采用1种网格密度（即N_m= 640000,粗颗粒直径与网格长度比为27）。

图 3 基准渗流模型及边界条件

Figure 3. CFD-DEM model for benchmark of seepage flow

下载: 全尺寸图片幻灯片

模拟分为两个阶段,t为0～0.5 s底部入口水流速度U_in从0 cm/s 增加至0.5 cm/s,之后保持底面入口水流速度不变,模拟至t=2 s,直到在土体内获得稳定的渗流场,产生水力梯度值为i=0.41。

为了对比模拟的计算效率,两种方法模拟所耗CPU时间对比如图4所示。可以看出全解流方法所用CPU时间为T_CPU=14.5 h；采用相同网格密度时（N_m= 640000）,任意解流方法用时为T_CPU=0.5 h,大幅降低计算时间。在任意解流方法中可进一步采用较低的网格密度（即N_m为34200,156250）进行模拟计算获得相近的模拟结果,而所用CPU时间更少,T_CPU分别为0.15,0.3 h。由此可见,当采用本文发展的任意解流方法进行渗流模拟时,在获得同样模拟结果时,其计算效率较全解流方法有较大提升。

图 4 全解流与任意解流用时对比

Figure 4. Comparison in CPU time between full resolved and ARU method

下载: 全尺寸图片幻灯片

3. 砂土渗流流固耦合分析

本节将采用新发展的任意解流流固耦合数值方法对砂土（或树脂颗粒）中向上渗流现象进行模拟,通过室内试验结果对流固数值模型进行验证,进一步使用流固耦合模型对渗流特性进行细观分析。在流固耦合模拟计算中,将固体土颗粒分为粗颗粒和细颗粒两类：粗颗粒粒径为流体网格尺寸10倍或以上；而对于细颗粒粒径,通过采用新耦合模型的处理方法,其粒径区间更加广泛,可大于流体网格尺寸。

3.1 向上渗流试验

向上渗流试验装置布置如图5所示。试验段中布置两层土颗粒,底部土基由具有一定级配的石英砂（粒径区间d_bs=1～1.5 mm）堆积而成,填筑高度为40 mm,控制孔隙率为0.4；试验中将石英砂染成红色以为便于观察细颗粒迁移现象。土基上部覆盖为粗颗粒滤层,分别采用两种不同单一粒径（D_f为10,5 mm）球形玻璃颗粒,即粗颗粒与最大细颗粒粒径之比α_s分别为6.7,3.3,滤层堆积高度均为40 mm。在试验段的顶部和底部分别安装测压管,用于测定试样两端的水压力值。使用高清摄像机记录试样段中土颗粒迁移和土体变形情况,并采用PIV技术获取颗粒运移轨迹和速度。试验过程中,水流从装置底部均匀流入,通过调节变频自吸式水泵功率,逐渐增加水流流速,直至试验段中上部粗颗粒滤层发生较大沉降或隆起破坏。在施加每一级水流时,待水流发展稳定后,通过测压管记录整个土体试验段（粗粒滤层40 mm+细粒土基40 mm）的水压力差值ΔP_s,并根据单位时间内流过试样单位横截面积的出口流量计算试样中表观水流速度U_sup。

图 5 向上渗流试验布置图

Figure 5. Set-up of upward seepage tests

下载: 全尺寸图片幻灯片

3.2 流固耦合数值模型

采用任意解流CFD-DEM流固耦合方法,完全按照上述试验条件建立三维数值模型,模拟向上渗流过程,如图6所示。

图 6 向上渗流数值模型及边界条件

Figure 6. CFD-DEM model of seepage tests

下载: 全尺寸图片幻灯片

在数值模型中,将滤层颗粒建模为粗颗粒,土基颗粒为细颗粒,且滤层和土基颗粒的粒径与试验保持一致。流体计算区域采用立方体网格进行离散,网格尺寸S_m选取为滤层粗颗粒粒径的1/10,对于两种不同粒径滤层颗粒,相应的网格尺寸S_m分别为1,0.5 mm。模型底部设置流速边界条件,水流由试样底部均匀流入,顶部出口设置为零压力边界条件,侧壁设置为无滑动边界条件。在任意解流流固耦合方法中,采用修正高斯权函数计算流体网格孔隙率等耦合所需的欧拉场变量,将高斯权函数带宽b_w设置为细颗粒最大粒径的两倍,即b_w=3.0 mm。为保证足够的计算精度和效率,确定DEM计算步长Δt_DEM=2×10^-5 s；CFD计算步长为Δt_CFD=2×10^-4 s。数值模型材料及计算参数设置见表3～5。

表 3 滤层-土基数值模型中的颗粒属性

Table 3. Material properties in numerical simulations of seepage test

参数	数值	参数	数值
砂颗粒密度ρ_s/(kg·m^-3)	2650	玻璃-玻璃滑动摩擦μ_{s_gg}	0.1545
玻璃珠密度ρ_g/(kg·m^-3)	2450	砂-玻璃滑动摩擦μ_{s_sg}	0.3
砂颗粒杨氏模量E_s/Pa	2×10¹⁰	砂-砂滚动摩擦μ_{r_ss}	0.26
玻璃珠杨氏模量E_g/Pa	5×10¹⁰	玻璃-玻璃滚动摩擦μ_{r_gg}	0.045
砂颗粒泊松比 $ν_{s}$	0.2	砂-玻璃滚动摩擦μ_{r_sg}	0.1
玻璃珠泊松比 $ν_{g}$	0.25	砂颗粒回弹系数e_s	0.9
砂-砂滑动摩擦μ_{s_ss}	0.84	玻璃珠回弹系数e_g	0.9

下载: 导出CSV

| 显示表格

表 4 滤层-土基数值模型中的流体属性

Table 4. Computational settings in numerical simulations of seepage test

密度ρ_f/(kg·m^-3)	动力黏度μ_f/(kg·m^-1·s^-1)	DEM时间步长Δt_DEM/s
1000	1×10^-3	2×10^-5

下载: 导出CSV

| 显示表格

表 5 滤层-土基数值模型中的计算设置

Table 5. Computational settings in numerical simulations of seepage test

CFD时间步长Δt_CFD/s	流体网格尺寸L_m/mm	高斯权函数带宽b_w/mm
2×10^-4	1, 0.5	3.0

下载: 导出CSV

| 显示表格

3.3 模拟结果与分析

图7展示了流固耦合模型与渗流试验获得的水力梯度随表观渗流速度变化曲线。在表观渗流速度U_sup增加的初始阶段（U_sup<0.95 cm/s）,模型预测水力梯度随表观渗流速度的增加而增大,呈线性关系,经线性拟合,表明数值预测的斜率与试验结果斜率接近；当渗流速度增加到一定值时（U_sup=0.95 cm/s,即过渡渗流速度）,无论是数值模型预测还是试验结果都发现,水力梯度先呈现缓慢增长的非线性变化规律,表明土体渗透性提高；紧随其后便迅速降低,我们称其为临界水力梯度,此时土体发生显著变形或破坏。对于粗颗粒与最大细颗粒粒径比α_s=6.7土样,耦合数值模型预测达到临界水力梯度时的临界渗流速度为U_sup=1.32 cm/s,其与渗流试验获得的值U_sup=1.44 cm/s接近,从而初步验证了本文所发展的耦合模型的准确性和有效性。对于粗颗粒与最大细颗粒粒径比α_s=3.3土样,水力梯度随渗流速度变化具有相似的性质,但其达到临界水力梯度时渗流速度U_sup=1.65 cm/s（数值预测）和U_sup=1.50 cm/s（试验结果）,较粒径比α_s=6.7的土样结果均有所增加。这与粗粒滤层颗粒粒径减小进而增加反滤效应有关。

图 7 水力梯度、滤层变形随表观渗流速度变化关系

Figure 7. Relationship between hydraulic gradient and coarse-particle layer deformation with superficial flow velocity

下载: 全尺寸图片幻灯片

进一步考察随着表观渗流速度的增加,粗粒滤层和细粒土基的变形情况。如图8（a）所示,对于粗细粒径比α_s=6.7土样,当表观渗流速度初始增加阶段（未达到临界水流速度）,可发现少量的土基细颗粒逐渐侵入滤层中；当渗流速度达到临界值（U_sup=1.32 cm/s）时,即达到临界水力梯度,大量细颗粒侵入甚至贯穿粗颗粒滤层；此时,粗粒滤层也发生了较大的下沉,如图8中颗粒位置变化所示,取滤层右下方具有较大下沉量颗粒作为参考,其最大下沉量可达12 mm（图7（a））。数值模拟与渗流试验获得的现象十分接近,再次验证了本文数值模型的准确性。对于粗颗粒与最大细颗粒粒径比α_s=3.3土样,不同地是,当渗流速度达到临界值之后,粗颗粒滤层迅速隆起；下部细粒土基出现较大膨胀变形,隆起量见图7（b）,可达约23 mm并逐渐呈现流化状态,如图8（b）所示。

图 8 不同阶段粗颗粒滤层变形随表观渗流速度变化情况

Figure 8. Evolution of coarse-particle layer deformation with superficial flow velocity

下载: 全尺寸图片幻灯片

进一步分析数值模拟结果发现,对于粗细粒径比α_s=6.7土样,在渗流达到过渡水流速度后,进入和穿越滤层的土基细颗粒堆积体孔隙率相对较大,堆积状态较为松散,如图9（b）所示,该现象很好地揭示了渗流达到过渡速度时（U_sup=0.95 cm/s）,土体渗透性提高的细观机理；当渗流速度进一步增加,超越临界渗流速度（U_sup=1.32 cm/s）时,细颗粒完全侵入和贯穿粗颗粒滤层,引起滤层更大的沉降发生,土体发生变形破坏,并呈现出更大的孔隙率,如图9（c）所示,从而引起水力梯度陡降。

图 9 土层孔隙率分布随表观渗流速度变化情况（粗细颗粒粒径比α_s=6.7）

Figure 9. Evolution of porosity profile with superficial flow velocity (diameter ratio, α_s=6.7)

下载: 全尺寸图片幻灯片

粗颗粒滤层发生的沉降过程,可进一步通过颗粒接触力链分布的演化情况从细观层面上进行解释。图10展示了土基细颗粒与滤层粗颗粒间、土基细颗粒间的接触力链分布演化结果（粗细粒径比α_s=6.7土样）。力链的粗细表示粒间接触力的大小；力链越粗,表明粒间接触力越大。当渗流试验开始前土样处于静水压力状态下（U_sup=0）,滤层与土基交界面处存在稳定、均匀的接触力链（图10（a））,表明此时土基能够为滤层提供足够的支撑。随着渗流流速增加直至临界渗流速度（U_sup=1.32 cm/s）时,如图10（b）所示,滤层与土基交界面附近力链明显减弱,这是由于在向上渗流作用下,土基细颗粒与滤层粗颗粒粒径差异较大,土基颗粒无法在滤层孔隙中形成稳定的土拱支撑滤层粗颗粒。进一步地,土基颗粒侵入上部滤层,土基无法为滤层提供足够支撑,如图10（c）所示,此时交界面附近力链基本消失,引起滤层粗颗粒下沉失稳。

图 10 土基颗粒力链分布随表观渗流速度变化情况（粗细颗粒_.粒径比α_s=6.7）

Figure 10. Evolution of force chain in fine-particle layer with superficial flow velocity (diameter ratio, α_s=6.7)

下载: 全尺寸图片幻灯片

4. 结论

本文为解决传统流固耦合方法对流体网格细化要求高和对粗细颗粒粒径比存在限制的问题,基于已有的半解流流固耦合方法,提出并发展新的任意解流流固耦合方法及数值模型,并通过模拟与分析砂土中向上渗流问题对模型进行了验证,检验了模型的准确性,也展示了其在砂土渗流细观尺度分析上的优势,获得3点结论。

（1）基于所提出的任意权函数框架,对高斯权函数进行修正使其满足流固耦合计算的筛选条件和湮灭条件；将修正高斯权函数与半解流流固耦合方法相结合,建立了任意解流流固耦合方法并开发了相应的模型程序。

（2）采用任意解流流固耦合方法模拟了滤层粗颗粒与土基细颗粒试样向上渗流过程,通过与室内向上渗流试验结果进行对比,验证了流固耦合方法及模型的有效性和准确性,表明任意解流流固耦合方法能够合理反映不同粗细颗粒粒径比下水-土耦合相互作用特性。本文所探讨的粗细颗粒粒径比（6.7,3.3）能够解除之前半解流方法对粒径比的限制（8～10以上）,虽尚未采用连续级配粒径进行模拟,但可以说明对于连续级配粒径模拟是有效的,由于篇幅有限,将在以后研究中进行详细说明。

（3）任意解流耦合数值模型再现了室内渗流试验中水力梯度随渗流速度变化规律。同时,基于数值模拟能够从细观层面上对滤层变形和细颗粒迁移进行深入分析；模拟获得的孔隙率与粒间接触力链分布变化较好地描述了渗流作用对滤层变形与细颗粒迁移的影响规律。本文所建立的任意解流流固耦合数值方法为深入研究砂土渗流问题提供了有效工具和手段。

图 1 高斯权函数支撑域示意图

Figure 1. Diagram of supporting domain of Gaussian-based weighting function

下载: 全尺寸图片幻灯片

图 2 任意解流耦合数值求解流程图

Figure 2. Flow chart of numerical solution of ARU CFD-DEM coupling

下载: 全尺寸图片幻灯片

图 3 基准渗流模型及边界条件

Figure 3. CFD-DEM model for benchmark of seepage flow

下载: 全尺寸图片幻灯片

图 4 全解流与任意解流用时对比

Figure 4. Comparison in CPU time between full resolved and ARU method

下载: 全尺寸图片幻灯片

图 5 向上渗流试验布置图

Figure 5. Set-up of upward seepage tests

下载: 全尺寸图片幻灯片

图 6 向上渗流数值模型及边界条件

Figure 6. CFD-DEM model of seepage tests

下载: 全尺寸图片幻灯片

图 7 水力梯度、滤层变形随表观渗流速度变化关系

Figure 7. Relationship between hydraulic gradient and coarse-particle layer deformation with superficial flow velocity

下载: 全尺寸图片幻灯片

图 8 不同阶段粗颗粒滤层变形随表观渗流速度变化情况

Figure 8. Evolution of coarse-particle layer deformation with superficial flow velocity

下载: 全尺寸图片幻灯片

图 9 土层孔隙率分布随表观渗流速度变化情况（粗细颗粒粒径比α_s=6.7）

Figure 9. Evolution of porosity profile with superficial flow velocity (diameter ratio, α_s=6.7)

下载: 全尺寸图片幻灯片

图 10 土基颗粒力链分布随表观渗流速度变化情况（粗细颗粒_.粒径比α_s=6.7）

Figure 10. Evolution of force chain in fine-particle layer with superficial flow velocity (diameter ratio, α_s=6.7)

下载: 全尺寸图片幻灯片

表 1 基准渗流模拟所采用材料属性

Table 1 Material properties in the benchmark of seepage flow

颗粒属性						流体属性
颗粒密度ρ_s/(kg·m^-3)	滑动摩擦系数μ_r	杨氏模量E/Pa	回弹系数e	泊松比 $ν$	滚动摩擦系数μ_f	密度ρ_f/(kg·m^-3)	动力黏度μ_f/(kg·m^-1·s^-1)
2650	0.84	2×10⁷	0.9	0.2	0.26	1000	1×10^-3

下载: 导出CSV

表 2 基准渗流模拟相关设置

Table 2 Computational settings in the benchmark of seepage flow

计算设置		任意解流设置		全解流设置
DEM时间步长Δt_DEM/s	CFD时间步长Δt_CFD/s	流体网格尺寸L_m/mm	高斯权函数带宽b_w/mm	流体网格尺寸L_m/mm
2×10^-5	2×10^-4	1.32, 0.8, 0.5	8.0	0.5

下载: 导出CSV

表 3 滤层-土基数值模型中的颗粒属性

Table 3 Material properties in numerical simulations of seepage test

参数	数值	参数	数值
砂颗粒密度ρ_s/(kg·m^-3)	2650	玻璃-玻璃滑动摩擦μ_{s_gg}	0.1545
玻璃珠密度ρ_g/(kg·m^-3)	2450	砂-玻璃滑动摩擦μ_{s_sg}	0.3
砂颗粒杨氏模量E_s/Pa	2×10¹⁰	砂-砂滚动摩擦μ_{r_ss}	0.26
玻璃珠杨氏模量E_g/Pa	5×10¹⁰	玻璃-玻璃滚动摩擦μ_{r_gg}	0.045
砂颗粒泊松比 $ν_{s}$	0.2	砂-玻璃滚动摩擦μ_{r_sg}	0.1
玻璃珠泊松比 $ν_{g}$	0.25	砂颗粒回弹系数e_s	0.9
砂-砂滑动摩擦μ_{s_ss}	0.84	玻璃珠回弹系数e_g	0.9

下载: 导出CSV

表 4 滤层-土基数值模型中的流体属性

Table 4 Computational settings in numerical simulations of seepage test

密度ρ_f/(kg·m^-3)	动力黏度μ_f/(kg·m^-1·s^-1)	DEM时间步长Δt_DEM/s
1000	1×10^-3	2×10^-5

下载: 导出CSV

表 5 滤层-土基数值模型中的计算设置

Table 5 Computational settings in numerical simulations of seepage test

CFD时间步长Δt_CFD/s	流体网格尺寸L_m/mm	高斯权函数带宽b_w/mm
2×10^-4	1, 0.5	3.0

下载: 导出CSV

参考文献(20)

[1]	周健, 姚志雄, 张刚. 砂土渗流过程的细观数值模拟[J]. 岩土工程学报, 2007, 29(7): 977-981. doi: 10.3321/j.issn:1000-4548.2007.07.004 ZHOU Jian, YAO Zhi-xiong, ZHANG Gang. Mesomechanical simulation of seepage flow in sandy soil[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(7): 977-981. (in Chinese) doi: 10.3321/j.issn:1000-4548.2007.07.004
[2]	EL SHAMY U, ZEGHAL M, DOBRY R, et al. Micromechanical aspects of liquefaction-induced lateral spreading[J]. International Journal of Geomechanics, 2010, 10(5): 190-201. doi: 10.1061/(ASCE)GM.1943-5622.0000056
[3]	ZHAO J D, SHAN T. Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[J]. Powder Technology, 2013, 239: 248-258. doi: 10.1016/j.powtec.2013.02.003
[4]	蒋明镜, 张望城. 一种考虑流体状态方程的土体CFD-DEM耦合数值方法[J]. 岩土工程学报, 2014, 36(5): 793-801. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm JIANG Ming-jing, ZHANG Wang-cheng. Coupled CFD-DEM method for soils incorporating equation of state for liquid[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5): 793-801. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm
[5]	王胤, 艾军, 杨庆. 考虑粒间滚动阻力的CFD-DEM流-固耦合数值模拟方法[J]. 岩土力学, 2017, 38(6): 1771-1780. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201706027.htm WANG Yin, AI Jun, YANG Qing. A CFD-DEM coupled method incorporating soil inter-particle rolling resistance[J]. Rock and Soil Mechanics, 2017, 38(6): 1771-1780. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201706027.htm
[6]	HAGER A, KLOSS C, PIRKER S, et al. Parallel resolved open source CFD-DEM: method, validation and application[J]. The Journal of Computational Multiphase Flows, 2014, 6(1): 13-27. doi: 10.1260/1757-482X.6.1.13
[7]	WANG Z L, FAN J R, LUO K. Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles[J]. International Journal of Multiphase Flow, 2008, 34(3): 283-302. doi: 10.1016/j.ijmultiphaseflow.2007.10.004
[8]	ANDERSON T B, JACKSON R. Fluid mechanical description of fluidized beds. equations of motion[J]. Industrial & Engineering Chemistry Fundamentals, 1967, 6(4): 527-539.
[9]	TSUJI Y, KAWAGUCHI T, TANAKA T. Discrete particle simulation of two-dimensional fluidized bed[J]. Powder Technology, 1993, 77(1): 79-87. doi: 10.1016/0032-5910(93)85010-7
[10]	CHENG K, WANG Y, YANG Q. A semi-resolved CFD-DEM model for seepage-induced fine particle migration in gap-graded soils[J]. Computers and Geotechnics, 2018, 100: 30-51. doi: 10.1016/j.compgeo.2018.04.004
[11]	DI RENZO A, DI MAIO F P. Comparison of contact-force models for the simulation of collisions in DEM-based granular flow codes[J]. Chemical Engineering Science, 2004, 59(3): 525-541. doi: 10.1016/j.ces.2003.09.037
[12]	DI FELICE R. The voidage function for fluid-particle interaction systems[J]. International Journal of Multiphase Flow, 1994, 20(1): 153-159. doi: 10.1016/0301-9322(94)90011-6
[13]	WANG Y, ZHOU L X, WU Y, et al. New simple correlation formula for the drag coefficient of calcareous sand particles of highly irregular shape[J]. Powder Technology, 2018, 326: 379-392. doi: 10.1016/j.powtec.2017.12.004
[14]	王胤, 周令新, 杨庆. 基于不规则钙质砂颗粒发展的拖曳力系数模型及其在细观流固耦合数值模拟中应用[J]. 岩土力学, 2019, 40(5): 2009-2015. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201905046.htm WANG Yin, ZHOU Ling-xin, YANG Qing. New drag coefficient model for irregular calcareous sand particles and its application into fluid-particle coupling simulation[J]. Rock and Soil Mechanics, 2019, 40(5): 2009-2015. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201905046.htm
[15]	YANG Q, CHENG K, WANG Y, et al. Improvement of semi-resolved CFD-DEM model for seepage-induced fine-particle migration: Eliminate limitation on mesh refinement[J]. Computers and Geotechnics, 2019, 110: 1-18. doi: 10.1016/j.compgeo.2019.02.002
[16]	JASAK H. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows[D]. London: Imperial College, University of London, 1996.
[17]	ISSA R I. Solution of the implicitly discretised fluid flow equations by operator-splitting[J]. Journal of Computational Physics, 1986, 62(1): 40-65. doi: 10.1016/0021-9991(86)90099-9
[18]	SHIRGAONKAR A A, MACIVER M A, PATANKAR N A. A new mathematical formulation and fast algorithm for fully resolved simulation of self-propulsion[J]. Journal of Computational Physics, 2009, 228(7): 2366-2390.
[19]	KLOSS C, GONIVA C. LIGGGHTS - open source discrete element simulations of granular materials based on lammps[M]//Supplemental Proceedings. Hoboken: John Wiley & Sons, Inc., 2011: 781-788.
[20]	WELLER H. OpenFOAM: The open source CFD toolbox, Programmer’s Guide[EB/OL]. UK: CFD Direct Ltd, 2009, Http://www.openfoam.com.

施引文献(10)

期刊类型引用(5)

1.	肖劲卿，温松诚，郭源. 基于单孔冲刷测试的黏性土抗侵蚀性及各向异性试验研究. 岩土力学. 2025(01): 187-198 . 百度学术
2.	蔡国庆，刁显锋，杨芮，王北辰，高帅，刘韬. 基于CFD-DEM的流-固耦合数值建模方法研究进展. 哈尔滨工业大学学报. 2024(01): 17-32 . 百度学术
3.	毛佳，余健坤，邵琳玉，赵兰浩. 三维可变形圆化多面体离散单元法. 岩土力学. 2024(03): 908-916 . 百度学术
4.	曹洋，刘杨，张超宇，杨俊杰，李国政. 基于离散元法的盾尾同步注浆扩散及参数优化研究. 岩土工程学报. 2024(10): 2119-2128 . 本站查看
5.	王安庆. 基于流固耦合的海上风机桩基础受力性能分析. 电工技术. 2022(23): 19-23 . 百度学术