基于浸入边界法的高解析度CFD-DEM流固耦合方法

doi:10.16183/j.cnki.jsjtu.2022.095

基于浸入边界法的高解析度CFD-DEM流固耦合方法

毛佳, 肖景文, 赵兰浩^,, 底瑛棠

河海大学水利水电学院,南京 210098

A Resolved CFD-DEM Approach Based on Immersed Boundary Method

MAO Jia, XIAO Jingwen, ZHAO Lanhao^,, DI Yingtang

College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China

通讯作者: 赵兰浩,教授,博士生导师;E-mail:zhaolanhao@hhu.edu.cn.

责任编辑: 王一凡

收稿日期: 2022-04-1 修回日期: 2022-07-21 接受日期: 2022-07-27

基金资助:

国家自然科学基金资助项目(52009034)
霍英东青年教师基金(151073)
中央高校基本科研业务费专项基金(B210201036)

Received: 2022-04-1 Revised: 2022-07-21 Accepted: 2022-07-27

作者简介 About authors

毛佳(1991-),副教授,从事滑坡涌浪和流固耦合方面的研究.

摘要

根据流固耦合问题的普遍性,提出了基于浸入边界法的高解析度计算流体力学-离散元法(CFD-DEM)流固耦合方法.新方法采用基于欧拉框架的计算流体力学方法描述流体的运动,采用基于拉格朗日框架的离散元法描述固体的运动及碰撞,在离散单元的表面布置浸入边界点,解决固体运动过程中与流体间的移动且未知的边界问题.为验证方法的准确性,模拟了圆柱绕流涡激振动、方块驰振两个经典算例,计算结果与数值解吻合度高,说明新方法能够准确描述流固耦合作用.最终,将该方法应用于多块体沉降的模拟,结果表明新方法能够反映流场的复杂变化,有效处理包含大量任意形状离散块体碰撞的流固耦合问题.

关键词： 流固耦合; 浸入边界; 计算流体力学; 离散元法; 高解析度

Abstract

Based on the immersed boundary method, a resolved CFD-DEM algorithm is proposed to tackle fluid-solid interaction problems which widely exist. In the proposed method, the fluid filed is described by the computational fluid dynamics in the Eulerian framework, while the movement and collision of the solids are simulated by the discrete element method in the Lagrangian framework. In order to deal with the moving interfaces between the fluid and the solids, several immersed boundary points are allocated on the boundaries of the solids. Two classic test cases are calculated to verify the accuracy of the proposed method, including the vortex-induced vibration of a cylinder and the rotational galloping of a rectangular rigid body. Good agreements are achieved between the current results and those in previous references and the reliability of the present method in modelling the fluid-solid interaction problems are proved. Finally, the sedimentation of multiple solids is simulated and the ability of the proposed CFD-DEM method in solving the complex fluid field and the collision among the solids with arbitrary shapes are verified.

Keywords： fluid-solid interaction; immersed boundary method; computational fluid dynamics (CFD); discrete element method (DEM); high resolution

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

本文引用格式

毛佳, 肖景文, 赵兰浩, 底瑛棠. 基于浸入边界法的高解析度CFD-DEM流固耦合方法[J]. 上海交通大学学报, 2023, 57(8): 988-995 doi:10.16183/j.cnki.jsjtu.2022.095

MAO Jia, XIAO Jingwen, ZHAO Lanhao, DI Yingtang. A Resolved CFD-DEM Approach Based on Immersed Boundary Method[J]. Journal of Shanghai Jiaotong University, 2023, 57(8): 988-995 doi:10.16183/j.cnki.jsjtu.2022.095

流固耦合问题常见于不同工程领域,根据固体的物理性质,可划分为流体与连续介质、流体与非连续介质间的流固耦合问题.目前,处理流体与连续介质间的流固耦合问题的数值模拟技术已较为成熟,广泛应用于水利工程^[1]、船舶工程^[2]等领域;而流体与非连续介质间的流固耦合,例如滑坡涌浪^[3]、冰区航行^[4]等问题的数值模拟仍具有较大的挑战性,主要难点在于固体运动过程中流固耦合界面的准确描述.

有研究者采用特殊的网格处理技术模拟流固耦合运动界面,如动网格^[5-6]、嵌套网格^[7]等.但由于固体运动的无规律性以及运动过程中相互间的碰撞,采用动网格方法时将会产生网格扭曲,甚至网格拓扑结构发生变化的问题.嵌套网格方法虽然适应性较好,但是此方法在运动固体周围设置的部件网格限制了其应用范围,当计算大量块体碰撞的问题时,各个部件网格将发生重叠,导致计算难以顺利进行.因此,这些网格处理技术只适用于模拟少量块体的非碰撞流固耦合问题,而对于涉及大量离散块体的流固耦合问题则存在明显的局限性.

考虑到流固耦合运动界面的模拟难度较高,有研究者提出简化的欧拉-欧拉多相流模型^[8],采用流体本构模型模拟离散块体,使得固体和流体间的界面成为内部界面,避免了采用特殊的措施处理流体域的移动且未知的边界.然而,此类方法不能体现固体的非连续性,较适用于解决固体呈现出显著流体化特征的问题.实际上,流体与固体具有截然不同的物理力学性质,需要采用不同的方程进行描述^[9],因此,有研究者提出了属于欧拉-拉格朗日方法的离散相模型(DPM)^[10].此方法能够考虑固体的非连续性,但由于无法模拟固体间的相互碰撞,只能用于模拟稀疏两相流问题.

近年来,国内外研究者普遍采用耦合数值方法模拟流固耦合问题,即将基于欧拉框架的计算流体力学(CFD)方法与基于拉格朗日框架的离散元法(DEM)耦合,建立CFDEM流固耦合数值模拟方法.目前的耦合方法大多基于经验公式计算流固耦合作用力,实际上可以将这些耦合数值模拟方法看作DPM的一种改进,通常将此类方法称为非解析的(Unresolved)CFDEM耦合方法^[11-12].非解析的CFDEM耦合方法是将计算流体力学方法与非连续介质方法相结合进行流固耦合问题研究的一次有益尝试,以简化的方式处理固体与流体的相互作用力,回避了两种介质之间边界移动且未知的问题,能够有效描述非连续介质在流场中的宏观运动规律,具有计算量小的优点.但与实际物理过程相比,该方法仍存在一些明显的不足:① 目前的经验公式仅限于圆形颗粒;② 流固耦合边界条件无法准确满足;③ 多颗粒流固耦合问题中,尾流对颗粒运动的影响难以充分反映.

为了解决非解析流固耦合方法的上述缺点,本文提出高解析度(Resolved)CFD-DEM流固耦合方法.采用Navier-Stokes方程组作为控制方程模拟流体的运动,并在固定网格上采用投影法求解流场运动变量;采用离散单元法模拟非连续固体的任意运动;引入浸入边界法(IBM)表示固体和流体间的移动边界,并计算两者间的相互作用力.通过强耦合迭代求解算法,提高计算的收敛性和稳定性,获得高精度的流场特性,准确反映固体运动过程中的流固耦合界面.

1 CFD-DEM流固耦合方法

相比于基于经验公式的非解析流固耦合方法,本文提出的高解析度方法能够准确反映流固耦合作用力,同时,在流固耦合界面上严格满足速度边界条件,能够在块体周围获得高解析度的流场.本文提出的基于浸入边界法的CFD-DEM流固耦合方法示意图如图1所示.

图1

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

图1 基于浸入边界法的CFD-DEM流固耦合方法示意图

Fig.1 Sketch of resolved CFD-DEM algorithm based on immersed boundary method

采用基于欧拉框架的CFD描述流体的运动,采用基于拉格朗日框架的DEM描述固体的运动及碰撞,在离散单元的表面布置浸入边界点,解决固体运动过程中与流体间的移动且未知的边界问题.为了实现流体域与固体域间的强耦合,提出流固耦合系统迭代求解方案.

1.1 流体域数值模拟

在黏性不可压缩Navier-Stokes方程组右端项中引入附加体力项,得到如下控制方程组:

(1)

\begin{matrix} \begin{matrix} Δ \cdot u = 0 \end{matrix} \end{matrix}

(2)

\frac{\partial u}{\partial t} + Δ \cdot (u u) = - \frac{1}{ρ} Δ p + \frac{1}{ρ} Δ \cdot τ + f_{b} + f

式中: 为流场速度; 为时间;p为压力;ρ为密度;$f_b$为体力项;τ为黏性应力张量;f为附加体力项.

在直接力法中,附加体力项由动量方程结合浸入边界点上的速度相容条件推导而来.速度相容条件是指在流体固体交界面上需要满足速度相同条件,即

(3)

\begin{matrix} U^{n + 1} = V^{n + 1} \end{matrix}

式中:Uⁿ⁺¹表示将笛卡尔网格点上的流体速度插值到浸入边界点上得到的速度;Vⁿ⁺¹表示同一浸入边界点上的固体速度;n表示计算步.

引入插值函数 $I (ϕ, X_{i}), 将定义在笛卡尔网格上的物理量插值到浸入边界点上, 同时, 引入分配函数 D (Φ, x),$ 将定义在浸入边界点上的物理量分配到笛卡尔网格上,插值函数与分配函数分别满足如下条件:

(4)

\begin{matrix} Φ (X_{i}) = I (ϕ, X_{i}) = \sum_{x \in g_{b}} ϕ (x) δ (x - X_{i}) \end{matrix}

(5)

\begin{matrix} ϕ (x) = D (Φ, x) = \overset{N_{I B P}}{\sum_{i = 1}} Φ (X_{i}) δ (x - X_{i}) Δ S_{i} \end{matrix}

式中: 表示定义在笛卡尔网格节点x上的物理量,包括u,p,f等; 表示定义在浸入边界点$X_i $上的物理量; $g_b $表示笛卡尔网格点的集合; 为插值计算某笛卡尔网格节点上物理量时所涉及到的浸入边界点数目; 表示浸入边界点的控制域;δ表示狄拉克函数,在二维情况下满足,

(6)

d (r) = \{\begin{array}{l} \frac{1}{8 h} [3 - 2 \frac{|r|}{h} + \\ \sqrt{1 + 4 \frac{|r|}{h} - 4 {(\frac{|r|}{h})}^{2}}], \\ |r| \leq h \\ 0, |r| \geq 2 h \\ \frac{1}{8 h} [5 - 2 \frac{|r|}{h} - \\ \sqrt{- 7 + 12 \frac{|r|}{h} - 4 {(\frac{|r|}{h})}^{2}}], \\ h < |r| < 2 h \end{array}

h取为网格尺寸的倍数;r表示有限元网格节点与浸入边界点的距离.结合式(4)~(6),可实现物理量的插值与分配.

基于特征线法,式(2)的时间离散格式如下:

(7)

Δ u = K^{n} + R^{n + 1} + f^{n + 1} Δ t

(8)

K^{n} = - Δ t [Δ \cdot (u u) - Δ \cdot \frac{τ}{ρ}] + \frac{Δ t^{2}}{2} u \cdot Δ [Δ \cdot (u u)]

(9)

\begin{array}{l} R^{n + 1} = - Δ t (\frac{1}{ρ} Δ p^{n + 1} - f_{b}) + \frac{Δ t^{2}}{2} u \cdot Δ (\frac{1}{ρ} Δ p^{n + 1} - f_{b}) \end{array}

在此基础上,将第n个时间步的速度增量分为两个部分:

(10)

\begin{matrix} Δ u^{*} = K^{n} + f^{n + 1} Δ t \end{matrix}

(11)

\begin{matrix} Δ u^{* *} = R^{n + 1} \end{matrix}

将中间速度场增量式(10)与速度场增量的修正量式(11)代入速度边界条件 $V^{n + 1} = U^{n + 1} 中,$ 整理可得:

(12)

\begin{matrix} f^{n + 1} Δ t = D [I (f^{n + 1} Δ t)] = D [V^{n + 1} - I (u^{n} + K^{n} + R^{n + 1})] \end{matrix}

由式(12)可见,附加体力项与速度场、压力场耦合,因此需要通过迭代计算求解上述三类未知量,具体步骤如下.

步骤1 忽略压力及附加体力项,计算中间速度场 $u^{*, K}, u^{*, K} = u^{n} + K^{n} .$

步骤2 计算压力场增量 $Δ p^{n, k}$ 及附加体力项 $f^{n + 1, k}$ .令 $p^{n + 1,0} = p^{n}$ , $f^{n + 1, 0} = 0$ , $k = 1$ .

步骤2.1 根据附加体力项的计算公式计算fⁿ^+1,^k.

值得注意的是,当浸入边界点与网格点重合,式(12)中的fⁿ⁺¹Δt=D[I(fⁿ⁺¹Δt)]成立,附加体力项的计算无需迭代.然而,通常情况下,浸入边界点与网格点不重合,通过单次计算附加体力项并不能满足速度边界条件,因此需要迭代求解体力项.令fⁿ^+1,^k^,0=fⁿ^+1,^k^-1,i=1,根据fⁿ^+1,^k^,ⁱΔt=fⁿ^+1,^k^,ⁱ^-1×Δt+D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹+fⁿ^+1,^k^,ⁱ^-1Δt)]计算fⁿ^+1,^k^,ⁱ.在此基础上,计算‖Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹+fⁿ^+1,^k^,ⁱ^-1Δt)‖,进行收敛性判断,其中‖·‖表示任一范数,如果范数小于容差ε,则结束步骤2.1,否则令i=i+1,迭代计算直到满足收敛条件.最终,fⁿ^+1,^k=fⁿ^+1,^k^,ⁱ.

步骤2.2 更新中间速度场u^*,计算式为u^*=u^*^,^K+fⁿ^+1,^kΔt.

步骤2.3 由 $Δ \cdot R^{n + 1, k} = - Δ \cdot (u^{n} + K^{n} + f^{n + 1, k} Δ t)$ 计算压力场增量 $Δ p^{n, k}$ .

步骤2.4 收敛性检查.计算‖I(uⁿ+Kⁿ+Rⁿ^+1,^k)-I(uⁿ+Kⁿ+Rⁿ^+1,^k^-1)‖,如果小于容差ε,则满足收敛性检查,结束步骤2的迭代计算,否则,令k=k+1,重复步骤2.1~2.4,直到满足收敛条件.

步骤3 更新压力场 $p^{n + 1} = p^{n} + Δ p^{n, k}$ ,校正速度场 $u^{n + 1} = u^{*} + R^{n + 1, k}$ .

1.2 固体域的数值模拟

根据牛顿第二定律,固体运动的控制方程如下:

(13)

m \ddot{δ} = F

(14)

I_{M} \ddot{θ} = M

式中:m为固体质量; $\ddot{δ}$ 为平动加速度;F为作用力在固体上的总力,表达式满足F=F_e+F_c+F_f-c $\overset{\cdot}{δ}$ ,F_e为外力,F_c为固体间接触力,F_f表示流体作用于固体的力,c为阻尼系数, $\overset{\cdot}{δ}$ 为平动速度;I_M为转动惯量张量; $\ddot{θ}$ 为角加速度;M为作用力在固体上的力矩.

F_c可分为法向接触力F_cn和切向接触力F_ct,F_cn的计算式如下:

(15)

\begin{matrix} F_{c n} = p_{f} \oint_{S_{β_{c} ⋂ β_{t}}} n_{S_{β_{c} ⋂ β_{t}}} (φ_{c} - φ_{t}) d S \end{matrix}

式中:p_f为罚函数;β_c为接触单元;β_t为目标单元;n为重叠域边界的单位外法向向量;φ_c、φ_t分别为定义在接触单元、目标单元上任一点的势函数^[13];S为重叠区域.

根据力-位移法则计算切向接触力,具体地,通过计算单元所受法向接触力合力及其作用点,并依据法向接触力合力的方向确定当前时间步切向力方向,以位移增量法计算切向接触力大小^[14],详细步骤如下:

步骤1 计算切向接触力增量 $Δ F_{c t}^{i} = k_{t} Δ δ_{t}^{i}$ , $k_{t}$ 为切向刚度系数, $Δ δ_{t}^{i}$ 为第i时间步的切向位移增量.

步骤2 将上一时间步的切向接触力转化至当前时间步切向接触力增量的方向,即 ${\bar{F}}_{c t}^{i - 1}$ =R $F_{c t}^{i - 1}$ ,R为旋转矩阵.

步骤3 计算当前时间步切向接触力为 $F_{c t}^{i}$ =max{ ${\bar{F}}_{c t}^{i - 1}$ +Δ $F_{c t}^{i}$ , $F_{c t}^{i}$ tan φ_μ},φ_μ表示最大静摩擦角.

1.3 流固耦合系统迭代求解方案

为了实现流体域与固体域间的强耦合,提高计算的收敛性和稳定性,获得高精度的流场特性,准确反映固体运动过程中的流固耦合界面,提出了流固耦合系统迭代求解方案.为简化起见,流体和固体控制方程迭代格式记为

(16)

\begin{matrix} (Δ u^{n, k}, Δ p^{n, k}) = R_{f} \end{matrix}

(17)

\begin{matrix} Δ δ^{n, k} = R_{s} \end{matrix}

式中:R_f表示流体方程右端项,与uⁿ、pⁿ、δⁿ^+1,^k^-1相关;R_s表示固体方程右端项,与δⁿ、uⁿ^+1,^k、pⁿ^+1,^k相关.

假设第n步计算已经完成,现在进入第n+1步的第k次迭代,执行过程如下:

步骤1 求解流体方程,更新速度场uⁿ^+1,^k与压力场pⁿ^+1,^k.

步骤2 求解固体方程,更新固体的位置δⁿ^+1,^k.

步骤3 收敛性检查.计算 $‖ Δ ζ^{n, k} - Δ ζ^{n, k - 1} ‖$ ,如果小于容差 $ε_{ζ}$ ,则满足收敛性检查,结束当前时间步的流固耦合系统迭代计算,否则,令 $k = k + 1$ ,重复步骤1~3,直到满足收敛条件.其中, $ζ$ 可取为变量 $Δ u$ 、 $Δ p$ 以及 $Δ δ$ , $ε_{ζ}$ 为对应于各变量ζ的给定的小值.

2 数值算例分析

2.1 圆柱绕流涡激振动

圆柱绕流涡激振动^[15-16]被广泛用于验证流固耦合算法解决固体平动耦合问题的准确性.圆柱具有水平向x及竖直向y两个方向的自由度,运动方程如下:

(18)

\begin{matrix} \ddot{X} + 2 ξ \frac{2 π}{U^{*}} \overset{\cdot}{X} + {(\frac{2 π}{U^{*}})}^{2} X = \frac{2}{π m^{*}} C_{D} \end{matrix}

(19)

\begin{matrix} \ddot{Y} + 2 ξ \frac{2 π}{U^{*}} \overset{\cdot}{Y} + {(\frac{2 π}{U^{*}})}^{2} Y = \frac{2}{π m^{*}} C_{L} \end{matrix}

式中: $ξ$ 为阻尼比; $U^{*}$ 为折合流速;D为圆柱直径; $m^{*} = m_{s} / m_{f}$ 为质量比, $m_{s}$ 、 $m_{f}$ 分别为圆柱质量与同等体积下的流体质量; $C_{L}$ 、 $C_{D}$ 分别为升力系数、阻力系数;X、Y为圆柱无量纲位移,分别满足 $X = x / D$ 、 $Y = y / D$ .

计算域大小为[-15D,45D]×[-20D,20D],初始时刻圆柱圆心坐标为(0,0).雷诺数Re=200,ξ=0.01,$U^*$=5.0, $m^*$=4/π.计算域左侧边界为均匀入流边界,右侧边界为自由出流边界,上、下边界为滑移固壁.为分析网格敏感性,采用了3种网格尺寸,分别记为L1、L2及L3,对应地,计算域离散为204×112、408×224及816×448个非结构化四边形单元,圆柱周围的网格尺寸分别为D/20、D/40及D/80,在圆柱表面布置753个浸入边界点.

固定圆柱直至下游产生稳定的卡门涡街,此后释放圆柱,可在圆柱下游观察到如图2所示的涡(图中Ω为涡量).圆柱振动频率为0.187,与漩涡脱落频率接近.

图2

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

图2 圆柱下游涡量场(L3)

Fig.2 Vortexes field downstream of cylinder (L3)

由于流场拖曳力的作用,圆柱的振动中心偏离(0, 0),本文方法计算得到的圆柱振动中心为(0.632D, 0),与参考文献[15]中的(0.62D, 0)接近.为方便对比,将圆柱振动中心平移至(0, 0),在此基础上绘制了圆柱“8字形”振动轨迹对比图,如图3所示.由于L1网格较粗,无法精确捕捉下游脱落的漩涡,所以使得x方向与y方向的位移较小,当采用较密网格L2、L3时,本文计算结果与前人的贴体网格计算结果^[15]、基于浸入边界法的强耦合数值模拟结果^[16]极为接近,一方面说明了加密网格可显著提高数值计算的精度,另一方面说明了本文方法模拟固体平动耦合问题的准确性.

图3

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

图3 圆柱中心振动轨迹

Fig.3 Centerline displacement phase plot

为说明流固耦合迭代对计算效率的影响,比较了计算时间t=3 s时,迭代次数取1和3的耗时.迭代1次的工况耗时为 7 899.38 s,迭代3次的工况耗时为 23 150.77 s,约为迭代1次工况耗时的2.9倍.可见,增加流固耦合迭代次数对计算效率有显著影响.然而,对于流场的变化或流固耦合作用剧烈的强耦合问题,若采用弱耦合求解格式,即用上一个时间步的力来求解当前步的运动,即使流体、固体的计算收敛,仍不能代表流固耦合系统一定收敛.因此,对于此类问题,即使会损失计算效率,也必须采用强耦合求解格式来保证计算精度.

2.2 方块驰振

为说明本文方法计算固体转动耦合问题的准确性,计算了经典的方块驰振问题^[16⇓-18],参考文献[17]中采用贴体网格计算,参考文献[16,18]中采用基于浸入边界法的流固耦合数值模拟方法计算,提供了强耦合模拟结果.方块宽度为L,厚度为D_s,运动方程如下:

(20)

\begin{matrix} \frac{d^{2} θ}{d t^{2}} + 2 ξ_{θ} \frac{2 π}{U_{θ}^{*}} \frac{d θ}{d t} + {(\frac{2 π}{U_{θ}^{*}})}^{2} θ = \frac{C_{T}}{2 I^{*}} \end{matrix}

式中:θ为方块转角;$ξ_θ$为方块阻尼比; $U_θ^*$为折合流速;C_T为转动力矩系数;I^*=I_s/(ρ_f)为转动惯量比,I_s为方块单位长度转动惯量,ρ_f为流体密度.

计算域尺寸为[-11D_s, 50D_s]×[-21D_s, 21D_s],初始时刻方块中心坐标为(0, 0).方块宽厚比满足L/D_s=4,雷诺数为Re=250,阻尼比ξ_θ=0.25,折合流速为 $U_{θ}^{*}$ =40,方块单位长度转动惯量为I_s=400.计算域左侧边界为均匀入流边界,右侧边界为自由出流边界,上、下边界为滑移固壁.计算域离散为347 236个非结构化四边形,方块周围网格尺寸为0.2D_s,在方块表面布置1 538个浸入边界点.

固定方块直至下游产生稳定脱落的涡,此后释放方块.方块转动过程中的流场涡量示意图如图4所示.当方块转动角度为正最大值及负最大值时,涡量场呈反对称状态,如图4(a)、图4(c)所示.

图4

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

图4 方块转动过程中的瞬时涡量场

Fig.4 Instantaneous vorticity surrounding rectangular rigid body

表1给出了本文计算得到的方块最大转角θ_max及转动频率f_rot,当前方法得出的最大转角为15.4°,与文献[16]中的结果误差不超过4.4%;当前方法得出的转动频率为0.019 2,与文献[18]中的结果误差不超过3.1%.

表1 方块最大转角 $θ_{m a x}$ 及转动频率 $f_{r o t}$

Tab.1 Maximum rotational angle $θ_{m a x}$ and rotational frequency $f_{r o t}$

来源	θ_max/(°)	f_rot
文献[16]	16.1	0.019 7
文献[17]	15.0	0.019 1
文献[18]	15.7	0.019 8
本文方法	15.4	0.019 2

新窗口打开| 下载CSV

图5中绘制了方块转角变化过程对比图.图中:U为流体自由流速.由于已有参考文献中方块释放时间不一致,初始阶段的方块运动有一定的差异,但是当转动逐渐发展稳定之后,各结果间的吻合度较高,说明了本文方法模拟转动耦合问题的准确性.

图5

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

图5 方块转角变化过程

Fig.5 Evolution of rotational angle

2.3 多块体沉降

为说明本文方法模拟流体与多固体间的流固耦合作用及固体间碰撞作用的能力,模拟了多块体沉降问题.

计算域尺寸为[10 m]×[10 m],采用两类正方形块体,共计28个,边长分别为0.698 m、0.548 m,在重力作用下自由下落.初始时刻位置如图6所示.A点坐标为(0.802 m, 7.945 m),B点坐标为(1.427 m, 6.750 m),同层块体间的x向距离均为1.100 m,相邻两层块体间y向距离为1.120 m.水体密度为1 000 kg/m³,动力黏性为0.01 kg/(m·s),块体密度为 2 600 kg/m³.计算域离散为480×480个结构化四边形单元,在边长为 0.698 m、0.548 m 的块体表面分别布置404、320个浸入边界点.本算例中所有边界均设为无滑固壁.

图6

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

图6 块体初始时刻位置

Fig.6 Initial positions of solids

图7中给出了块体运动速度(v)及流体域速度矢量场变化过程.初始阶段(见图7(a)),水体受块体运动扰动的同时,向块体施加了作用力.随着块体的进一步沉降,水体和块体间的相互作用逐渐增强(见图7(b)~图7(d)),此阶段在块体周围可以观察到大量的涡,受涡的影响,块体在沉降过程中伴随着转动.部分块体在2.8 s左右开始接触计算域底部,此后流场的扰动显著降低(见图7(e)~图7(f)).此算例说明了本文方法不限于模拟流体和圆形颗粒的相互作用,能够描述流场的复杂变化及固体间的碰撞,可实现复杂流固耦合问题中固体运动界面的准确描述.

图7

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

图7 块体运动速度及流体域速度矢量场变化过程

Fig.7 Instantaneous velocity vector of fluid field and position of rigid bodies

3 结论

提出了基于浸入边界法的高解析度CFD-DEM流固耦合方法,主要结论如下:

(1) 区别于非解析CFDEM耦合方法依赖于经验公式计算流固耦合作用力,新方法采用浸入边界法表示固体和流体间的移动边界,并计算两者间的相互作用力.

(2) 计算了圆柱绕流涡激振动、方块驰振两个经典算例,并与已有数值计算结果进行比较,验证了新方法计算平动耦合问题、转动耦合问题的准确性.

(3) 模拟了多块体沉降问题,说明了新方法不局限于圆形颗粒,可以有效描述复杂流固耦合问题中的运动界面,获得高解析度的流场.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李师尧, 程永光, 张春泽.

IB-LB耦合格式模拟贯流式水轮机三维瞬变流

[J]. 华中科技大学学报(自然科学版), 2016, 44(1): 122-127.

[本文引用: 1]

Shiyao

, CHENG

Yongguang

, ZHANG

Chunze

3-demensional transient flow simulation of tubluar turbine based on IB-LB

[J]. Journal of Huazhong University of Science and Technology (Nature Science Edition), 2016, 44(1): 122-127.

[本文引用: 1]

[2]

杨兆臣, 张怀新.

基于双向流固耦合的柔性表面覆盖层减阻性能

[J]. 上海交通大学学报, 2016, 50(8): 1165-1172.

[本文引用: 1]

YANG

Zhaochen

, ZHANG

Huaixin

Investigation of drag reduction of flexible surface based on bi-directional fluid-structure interaction

[J]. Journal of Shanghai Jiao Tong University, 2016, 50(8): 1165-1172.

[本文引用: 1]

[3]

徐文杰.

滑坡涌浪流-固耦合分析方法与应用

[J]. 岩石力学与工程学报, 2020, 39(7): 1-14.

[本文引用: 1]

Wenjie

Fluid-solid coupling method of landslide tsunamis and its application

[J]. Chinese Journal of Rock Mechanics and Engineering, 2020, 39(7): 1-14.

[本文引用: 1]

[4]

郭春雨, 李夏炎, 王帅, 等.

冰区航行船舶碎冰阻力预报数值模拟方法

[J]. 哈尔滨工程大学学报, 2016, 37(2): 145-150.

[本文引用: 1]

GUO

Chunyu

, LI

Xiayan

, WANG

Shuai

, et al.

A numerical simulation method for resistance prediction of ship in pack ice

[J]. Journal of Harbin Engineering University, 2016, 37(2): 145-150.

[本文引用: 1]

[5]

, BAYLY

A E

, HASSANPOUR

Coupling CFD-DEM with dynamic meshing: A new approach for fluid-structure interaction in particle-fluid flows

[J]. Powder Technology, 2018, 325: 620-631.

DOI:10.1016/j.powtec.2017.11.045 URL [本文引用: 1]

[6]

施红辉, 周栋, 温俊生, 等.

基于ALE方法的弹性圆柱壳入水时的流固耦合模拟

[J]. 弹道学报, 2020, 32(1): 9-14.

DOI:10.12115/j.issn.1004-499X(2020)01-002 [本文引用: 1]

为了获得头型及材料弹性对圆柱壳入水过程中头部变形、压力分布、入水空泡形态及入水运动状态的影响,采用ALE方法,基于有限元软件ANSYS\LS-DYNA,对平头、120°锥角、90°锥角弹性圆柱壳的入水过程进行了数值模拟,对平头弹性圆柱实体和平头刚性圆柱壳进行计算,并对模拟结果进行对比。结果表明:入水过程中空泡直径随着锥角的增大而增大,锥角越大,圆柱壳上表面的变形越大,最大变形出现的时间越早:平头弹性圆柱壳在0.1 ms时出现最大变形0.84 mm,120°锥角弹性圆柱壳在0.3 ms时出现最大变形0.54 mm,90°锥角弹性圆柱壳在0.5 ms时出现最大变形0.43 mm; 下表面受到的压力越大,圆柱壳的速度衰减越快; 平头圆柱壳下表面的变形与振动频率大于上表面; 上表面的变形是由惯性引起的,下表面的变形是由流体冲击力引起的。

SHI

Honghui

, ZHOU

Dong

, WEN

Junsheng

, et al.

Fluid-solid interaction simulation of elastic cylindrical shell penetrating water based on the ALE method

[J]. Journal of Ballistics, 2020, 32(1): 9-14.

DOI:10.12115/j.issn.1004-499X(2020)01-002 [本文引用: 1]

In order to obtain the influence of the head shape and material elasticity on the deformation,pressure distribution,cavity shape and motion state of cylindrical shell in the process of water entry,Arabic Lagrange Euler(ALE)method was used to simulate and analyze the water entry process of elastic cylindrical shell with flat head,120° cone angle and 90° cone angle based on the finite element software ANSYS\LS-DYNA. The simulation result of rigid cylindrical shell was compared with that of flat head projectile. The results show that the diameter of cavity increases with the increase of the cone angle. The larger the cone angle is,the larger the deformation of the top surface of cylindrical shell is,and the earlier the maximum deformation occur. The maximum deformation of flat elastic cylindrical shell is 0.84 mm at 0.1 ms,and the maximum deformation of elastic cylindrical shell with 120° cone angle is 0.54 mm at 0.3 ms,and the maximum deformation of elastic cylindrical shell with 90° cone angle is 0.43 mm at 0.5 ms. The greater the pressure on the bottom surface,the faster the velocity of cylindrical shell decays. The deformation and vibration frequency of the bottom surface of flat head cylindrical shell is higher than that of the top surface. The deformation of the top surface is caused by inertia,and the deformation of the bottom surface is caused by the impact force of the fluid.

[7]

Z H

, QIAN

, MARTINEZ-FERRER

P J

, et al.

An overset mesh based multiphase flow solver for water entry problems

[J]. Computers & Fluids, 2018, 172: 689-705.

DOI:10.1016/j.compfluid.2018.01.025 URL [本文引用: 1]

[8]

郭剑, 沈伟, 李同录, 等.

一种流动性滑坡涌浪动力学模型

[J]. 水科学进展, 2019, 30(2): 273-281.

[本文引用: 1]

GUO

Jian

, SHEN

Wei

, LI

Tonglu

, et al.

Establishment of dynamic model of a flow-like landslide-induced surge

[J]. Advances in Water Science, 2019, 30(2): 273-281.

[本文引用: 1]

[9]

黄波林, 殷跃平, 陈小婷, 等.

地质灾害涌浪研究综述

[J]. 水文地质工程地质, 2014, 41(1): 112-118.

[本文引用: 1]

HUANG

Bolin

, YIN

Yueping

, CHEN

Xiaoting

, et al.

A review on impulse wave generated by geohazards

[J]. Hydrogeology & Engineering Geology, 2014, 41(1): 112-118.

[本文引用: 1]

[10]

ORONA

J D

, ZORRILLA

S E

, PERALTA

J M

Computational fluid dynamics combined with discrete element method and discrete phase model for studying a food hydrofluidization system

[J]. Food and Bioproducts Processing, 2017, 102: 278-288.

DOI:10.1016/j.fbp.2017.01.005 URL [本文引用: 1]

[11]

王胤, 艾军, 杨庆.

考虑粒间滚动阻力的CFD-DEM流-固耦合数值模拟方法

[J]. 岩土力学, 2017, 38(6): 1771-1780.

[本文引用: 1]

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.

[本文引用: 1]

[12]

ZHAO

, XU

, GE

, et al.

A CFD-DEM-IBM method for Cartesian grid simulation of gas-solid flow in complex geometries

[J]. Chemical Engineering Journal, 2020, 389: 124343-1-14.

DOI:10.1016/j.cej.2020.124343 URL [本文引用: 1]

[13]

MUNJIZA

The combined finite-discrete element method[M]. Chichester, UK: Wiley, 2004.

[本文引用: 1]

[14]

刘勋楠, 赵兰浩, 毛佳, 等.

距离势离散单元法

[J]. 岩石力学与工程学报, 2017, 36(6): 1353-1362.

[本文引用: 1]

LIU

Xunnan

, ZHAO

Lanhao

, MAO

Jia

, et al.

Discrete element method based on distance potential

[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(6): 1353-1362.

[本文引用: 1]

[15]

BLACKBURN

H M

, KARNIADAKTS

G E

Two-and three-dimensional simulations of vortex-induced vibration of a circular cylinder

[C]//Proceedings of the Third (1993) International Offshore and Polar Engineering Conference. Golden, Colorado, USA: ISOPE, 1993: 715-720.

[本文引用: 3]

[16]

YANG

J M

, STERN

A non-iterative direct forcing immersed boundary method for strongly-coupled fluid-solid interactions

[J]. Journal of Computational Physics, 2015, 295: 779-804.

[本文引用: 5]

[17]

ROBERTSON

, LI

, SHERWIN

S J

, et al.

A numerical study of rotational and transverse galloping rectangular bodies

[J]. Journal of Fluids and Structures, 2003, 17(5): 681-699.

DOI:10.1016/S0889-9746(03)00008-2 URL [本文引用: 2]

[18]

YANG

J M

, STERN

A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions

[J]. Journal of Computational Physics, 2012, 231(15): 5029-5061.

DOI:10.1016/j.jcp.2012.04.012 URL [本文引用: 3]

IB-LB耦合格式模拟贯流式水轮机三维瞬变流

2016

... 流固耦合问题常见于不同工程领域,根据固体的物理性质,可划分为流体与连续介质、流体与非连续介质间的流固耦合问题.目前,处理流体与连续介质间的流固耦合问题的数值模拟技术已较为成熟,广泛应用于水利工程^[1]、船舶工程^[2]等领域;而流体与非连续介质间的流固耦合,例如滑坡涌浪^[3]、冰区航行^[4]等问题的数值模拟仍具有较大的挑战性,主要难点在于固体运动过程中流固耦合界面的准确描述. ...

IB-LB耦合格式模拟贯流式水轮机三维瞬变流

2016

基于双向流固耦合的柔性表面覆盖层减阻性能

2016

基于双向流固耦合的柔性表面覆盖层减阻性能

2016

滑坡涌浪流-固耦合分析方法与应用

2020

滑坡涌浪流-固耦合分析方法与应用

2020

冰区航行船舶碎冰阻力预报数值模拟方法

2016

冰区航行船舶碎冰阻力预报数值模拟方法

2016

Coupling CFD-DEM with dynamic meshing: A new approach for fluid-structure interaction in particle-fluid flows

2018

... 有研究者采用特殊的网格处理技术模拟流固耦合运动界面,如动网格^[5-6]、嵌套网格^[7]等.但由于固体运动的无规律性以及运动过程中相互间的碰撞,采用动网格方法时将会产生网格扭曲,甚至网格拓扑结构发生变化的问题.嵌套网格方法虽然适应性较好,但是此方法在运动固体周围设置的部件网格限制了其应用范围,当计算大量块体碰撞的问题时,各个部件网格将发生重叠,导致计算难以顺利进行.因此,这些网格处理技术只适用于模拟少量块体的非碰撞流固耦合问题,而对于涉及大量离散块体的流固耦合问题则存在明显的局限性. ...

基于ALE方法的弹性圆柱壳入水时的流固耦合模拟

2020

基于ALE方法的弹性圆柱壳入水时的流固耦合模拟

2020

An overset mesh based multiphase flow solver for water entry problems

2018

一种流动性滑坡涌浪动力学模型

2019

... 考虑到流固耦合运动界面的模拟难度较高,有研究者提出简化的欧拉-欧拉多相流模型^[8],采用流体本构模型模拟离散块体,使得固体和流体间的界面成为内部界面,避免了采用特殊的措施处理流体域的移动且未知的边界.然而,此类方法不能体现固体的非连续性,较适用于解决固体呈现出显著流体化特征的问题.实际上,流体与固体具有截然不同的物理力学性质,需要采用不同的方程进行描述^[9],因此,有研究者提出了属于欧拉-拉格朗日方法的离散相模型(DPM)^[10].此方法能够考虑固体的非连续性,但由于无法模拟固体间的相互碰撞,只能用于模拟稀疏两相流问题. ...

一种流动性滑坡涌浪动力学模型

2019

地质灾害涌浪研究综述

2014

地质灾害涌浪研究综述

2014

Computational fluid dynamics combined with discrete element method and discrete phase model for studying a food hydrofluidization system

2017

考虑粒间滚动阻力的CFD-DEM流-固耦合数值模拟方法

2017

... 近年来,国内外研究者普遍采用耦合数值方法模拟流固耦合问题,即将基于欧拉框架的计算流体力学(CFD)方法与基于拉格朗日框架的离散元法(DEM)耦合,建立CFDEM流固耦合数值模拟方法.目前的耦合方法大多基于经验公式计算流固耦合作用力,实际上可以将这些耦合数值模拟方法看作DPM的一种改进,通常将此类方法称为非解析的(Unresolved)CFDEM耦合方法^[11-12].非解析的CFDEM耦合方法是将计算流体力学方法与非连续介质方法相结合进行流固耦合问题研究的一次有益尝试,以简化的方式处理固体与流体的相互作用力,回避了两种介质之间边界移动且未知的问题,能够有效描述非连续介质在流场中的宏观运动规律,具有计算量小的优点.但与实际物理过程相比,该方法仍存在一些明显的不足:① 目前的经验公式仅限于圆形颗粒;② 流固耦合边界条件无法准确满足;③ 多颗粒流固耦合问题中,尾流对颗粒运动的影响难以充分反映. ...

考虑粒间滚动阻力的CFD-DEM流-固耦合数值模拟方法

2017

A CFD-DEM-IBM method for Cartesian grid simulation of gas-solid flow in complex geometries

2020

2004

... 式中:p_f为罚函数;β_c为接触单元;β_t为目标单元;n为重叠域边界的单位外法向向量;φ_c、φ_t分别为定义在接触单元、目标单元上任一点的势函数^[13];S为重叠区域. ...

距离势离散单元法

2017

... 根据力-位移法则计算切向接触力,具体地,通过计算单元所受法向接触力合力及其作用点,并依据法向接触力合力的方向确定当前时间步切向力方向,以位移增量法计算切向接触力大小^[14],详细步骤如下: ...

距离势离散单元法

2017

Two-and three-dimensional simulations of vortex-induced vibration of a circular cylinder

1993

... 圆柱绕流涡激振动^[15-16]被广泛用于验证流固耦合算法解决固体平动耦合问题的准确性.圆柱具有水平向x及竖直向y两个方向的自由度,运动方程如下: ...

... 由于流场拖曳力的作用,圆柱的振动中心偏离(0, 0),本文方法计算得到的圆柱振动中心为(0.632D, 0),与参考文献[15]中的(0.62D, 0)接近.为方便对比,将圆柱振动中心平移至(0, 0),在此基础上绘制了圆柱“8字形”振动轨迹对比图,如图3所示.由于L1网格较粗,无法精确捕捉下游脱落的漩涡,所以使得x方向与y方向的位移较小,当采用较密网格L2、L3时,本文计算结果与前人的贴体网格计算结果^[15]、基于浸入边界法的强耦合数值模拟结果^[16]极为接近,一方面说明了加密网格可显著提高数值计算的精度,另一方面说明了本文方法模拟固体平动耦合问题的准确性. ...

... [15]、基于浸入边界法的强耦合数值模拟结果^[16]极为接近,一方面说明了加密网格可显著提高数值计算的精度,另一方面说明了本文方法模拟固体平动耦合问题的准确性. ...

A non-iterative direct forcing immersed boundary method for strongly-coupled fluid-solid interactions

2015

... 圆柱绕流涡激振动^[15-16]被广泛用于验证流固耦合算法解决固体平动耦合问题的准确性.圆柱具有水平向x及竖直向y两个方向的自由度,运动方程如下: ...

... 为说明本文方法计算固体转动耦合问题的准确性,计算了经典的方块驰振问题^[16⇓-18],参考文献[17]中采用贴体网格计算,参考文献[16,18]中采用基于浸入边界法的流固耦合数值模拟方法计算,提供了强耦合模拟结果.方块宽度为L,厚度为D_s,运动方程如下: ...

... ]中采用贴体网格计算,参考文献[16,18]中采用基于浸入边界法的流固耦合数值模拟方法计算,提供了强耦合模拟结果.方块宽度为L,厚度为D_s,运动方程如下: ...

... 表1给出了本文计算得到的方块最大转角θ_max及转动频率f_rot,当前方法得出的最大转角为15.4°,与文献[16]中的结果误差不超过4.4%;当前方法得出的转动频率为0.019 2,与文献[18]中的结果误差不超过3.1%. ...

A numerical study of rotational and transverse galloping rectangular bodies

2003

... ,参考文献[17]中采用贴体网格计算,参考文献[16,18]中采用基于浸入边界法的流固耦合数值模拟方法计算,提供了强耦合模拟结果.方块宽度为L,厚度为D_s,运动方程如下: ...

A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions

2012

... ,18]中采用基于浸入边界法的流固耦合数值模拟方法计算,提供了强耦合模拟结果.方块宽度为L,厚度为D_s,运动方程如下: ...

〈

〉