AN AXISYMMETRIC IMMERSED BOUNDARY METHOD BASED ON 2D MESH
-
摘要: 浸没边界法是处理颗粒两相流中运动边界问题的一种常用数值模拟方法. 当研究的物理问题的无量纲参数满足一定要求时, 该流场结构呈现轴对称状态. 为此本文提出了一种基于2D笛卡尔网格和柱坐标系的轴对称浸没边界法. 该算法采用有限体积法(FVM)对动量方程进行空间离散, 并通过阶梯状锐利界面替代真实的固体浸没边界来封闭控制方程. 为了提高计算效率, 本文采用自适应网格加密技术提高浸没边界附近网格分辨率. 由于柱坐标系的使用, 使得动量方程中的黏性项产生多余的源项, 我们对其作隐式处理. 此外, 在对小球匀速近壁运动进行直接数值模拟时, 由于球壁间隙很小, 间隙内的压力变化比较剧烈. 因此想要精确地解析流场需要很高的网格分辨率. 此时, 需要在一个时间步内多次实施投影步来保证计算的稳定性. 而在小球自由碰壁运动中, 我们通过引入一个润滑力模型使得低网格分辨率下也能模拟小球近壁处的运动. 最后通过小球和圆盘绕流、Stokes流小球近壁运动以及小球自由下落碰壁弹跳算例验证本算法对于轴对称流的静边界和动边界问题均是适用和准确的.Abstract: Immersed boundary methods is a common numerical simulation method to deal with moving boundary problems in particle two-phase flows. When the dimensionless parameters of studied physical problems meet certain requirements, its corresponding flow structure becomes axisymmetric. Hence, an axisymmetric immersed boundary method based on a 2D mesh and cylindrical coordinates is proposed in the present paper. A finite volume method is used as the spatial discretization in the present algorithm. And the governing equation is closed by a sharp stepped interface, which is used to replace the real solid immersed boundary. In order to improve the efficiency of computation, an adaptive mesh refinement technology is used to improve the mesh resolution near the immersed boundary. The using of cylindrical coordinates will produce a redundant source item from the viscous term in the momentum equation. The additional source term will be handled by using an implicit scheme. Moreover, in the direct numerical simulation of a sphere approaching a wall with a constant velocity, the pressure of the fluid in the gap dramatically changes because of the small gap. So, in order to accurately analyze the flow field, the required grid resolution is very high in the gap. Multiple projection-step calculations are carried out in one time step for maintaining the stability of the simulation. In the movement of a sphere free impacting a wall, a lubrication force model will be introduced to simulate the movement of the sphere near the wall even with a low grid resolution. Finally, simulation results on the flow past a fixed sphere, the flow past a circular disk, the stokes flow by a sphere approaching a wall and the flow caused by the sphere-wall collision prove that the present axisymmetric IBM algorithm is applicable and accurate for dealing with stationary and moving boundary problems in an axisymmetric flow.
-
Keywords:
- axisymmetric /
- 2D mesh /
- immersed boundary method /
- sphere-wall collision
-
引 言
流体领域中的颗粒两相流[1-2]一直是工程和学术研究上的重要内容. 除了理论和实验研究之外, 数值模拟也是一种重要的研究手段. 而在颗粒两相流的相关研究中, 复杂几何和动边界问题是十分普遍的. 当我们采用数值模拟手段研究上述问题时, 使用贴体网格是有一定挑战的. 例如固体几何比较复杂, 那么壁面边界层附近的网格处理将会是麻烦的. 同样当固体边界发生运动时, 贴体网格会发生畸变, 使得网格的正交性变差. 但是浸没边界法(immersed boundary method, IBM)为相关问题提供了非常好的技术手段.
浸没边界法首先是由Peskin[3-4]为了模拟心脏力学以及相关的血液流动而发展出来一种方法. 该方法对流体和固体采用两套不同网格, 流体一般采用正交的欧拉网格而固体一般采用拉格朗日网格. 固体边界上的信息通过直接在Navier-Stokes方程中添加力源项反映到流场中, 并引入狄拉克δ函数来联系固体域和流体域. 随后此方法也被后续的科研工作者推广到了其他众多领域, 如文献[5]将IBM方法推广到了湍流领域、自然和强迫对流传热领域; 文献[6]将IBM方法分别应用到颗粒液滴运动问题以及可压缩黏性流动问题; Grigoriadis等[7]则将其推广到磁流体动力学(magnetohydrodynamics, MHD)领域. 另外, 为了对钠冷快堆中的停堆组件落棒过程进行数值模拟研究, 秦如冰等[8]开发了2D的浸没边界法. 同时, 近年来由于格子玻尔兹曼方法兴起, IBM方法与格子法的结合使用也是热门的领域. 周睿等利用格子玻尔兹曼浸没边界法分别对双层刚性植被明渠水流[9]以及变动水面水库[10]进行了数值模拟和验证. 佟莹等[11]发展了一种直接力格式的浸没边界格子模型用于动边界流动数值计算. 李桥忠等[12]提出了一种基于浸没边界−简化热格子玻尔兹曼方法的耦合模型, 该模型简化了边界条件的实现形式, 减小了计算成本.
通常根据浸没固体边界条件引入流体计算域的差异, 后续发展的IBM方法主要分为两类: 一类是连续力法如Peskin[3-4], 另一类是离散力法如文献[13]. 相比于连续力法, 离散力法并不通过对控制方程改动来施加浸没边界条件, 而是在浸没边界附近单元直接进行插值引入边界条件. 考虑到后者在满足物理量守恒特性上要优于前者, 本文所采用的IBM方法也是属于离散力法. 但是使用离散力方式处理浸没边界时, 固体界面将背景的流体网格切割成不规则形状. 如何处理切割后的网格单元的守恒性是一个难点. 文献[14]通过耦合IBM以及切割网格技术有效地处理了由于动边界引发的压力震荡, 但其方法实现是比较困难的. 随后Pan等[15]发展了一种相容守恒格式的浸没边界法用以研究MHD流动以及动边界问题. 该方法既能够抑制动边界产生的压力震荡又相对前者更易实现.
同时经过本文调研发现, 目前大部分关于IBM方法的研究[1-15]主要针对的是较为普遍的二维或三维问题, 而对于轴对称这种特定物理问题(比如雷诺数Re < 210时的小球绕流[16], 伽利略数G < 156时的小球自由运动[17])的相关研究是比较少的. Lai等[18]基于柱坐标发展了一种轴对称IBM用以研究伴随不溶性表面活性物质的轴对称界面流. Li等[19]在研究细胞分裂增长的轴对称过程中发展了一种轴对称浸没边界模型. Li[20]采用涡量−流函数耦合控制方程发展了一种快速轴对称浸没边界投影法用以处理小球颗粒与壁面正碰问题. 上述三种模型均是基于连续力法来实现的, 而目前基于离散力格式的轴对称IBM方法还未发展.
因此本文发展了一种基于有限体积法(FVM)以及离散力格式的2D轴对称浸没边界法. 此算法是以Pan等[15]所发展的2D和3D相容守恒浸没边界法扩展而来. 因此对于浸没边界条件的引入以及锐利阶梯界面的划分, 本文沿袭其处理方式. 但是Pan等[15]所发展的原始IBM方法并未考虑碰撞情况, 当需考虑颗粒壁面碰撞, 在小球离壁面只有一层网格时, 由于插值数据信息不足导致数值计算崩溃. 所以在其原始算法基础之上, 本文进一步完善了浸没边界离壁面边界只有一层网格时的情况.
如图1所示为例, 所有背景网格单元被分为三种: 第一种是棕色表示的固体网格单元, 其网格中心在浸没边界法内部; 第二种是绿色表示的IB网格单元, 其为离浸没边界最近的第一层网格单元. 第三种是剩下的蓝色表示的流体单元(直接参与离散矩阵计算). 除了上述网格单元之外, 黄色的壁面面元代表了流体域的壁面边界, 面心的数据由边界条件给出. 图中的红色实线则是表示浸没边界, 阶梯状虚线表示锐利界面(实际计算的边界).
![]()
图 1 浸没边界和网格单元划分示意图Figure 1. Schematic diagram of immersed boundary and dividing cells type由于IB单元一般是不规则的切割单元(如图1中所示, 绿色IB单元被红色浸没边界切割), 本文采用文献[21]提出的加权最小二乘法(LSM)插值得到. 而插值所用到的信息一般为浸没固体边界条件以及周围流体网格中心的物理量. 当浸没边界开始接触壁面边界时(相距只有几层网格), 流体侧壁面边界上的面心数据也会作为插值数据点. 得到IB网格单元中心数据之后, 再通过引入阶梯状锐利界面作为近似边界来封闭离散方程, 从而实现固体边界条件的引入.
一般根据物体运动速度的确定方式, 本文将其分为两种基本类型. 第一种是运动物体的强迫运动, 人为给出物体的速度随时间的变化. 二是通过耦合运动学方程和Navier-Stokes方程, 由流体−颗粒的相互作用来确定物体的速度. 通常, 运动学方程中的力系数在时间上是显示离散的, 以便将运动学方程从Navier-Stokes方程中解耦[22]. 本文通过2D轴对称IBM方法模拟Stokes流小球近壁运动以及小球自由下落碰壁弹跳分别对上述两种运动方式进行验证.
1. 数值模型
1.1 控制方程
本文算法考虑的是一种不可压缩的、运动黏性为
$ \nu $ 、密度为$ {\rho }_{f} $ 的流体所产生的无环流轴对称流动. 在轴对称假设之下, 本文可以利用柱坐标系将三维流动的控制方程简化成两维形式. 本文使用$\boldsymbol{u}=(u_{{r}}, u_{{z}})$ 以及$ p $ 来分别表示子午面上的速度和压力, 其中$u_r$ 和$u_z$ 分别表示径向和轴向的速度分量. 因此无量纲的流体动力学控制方程可以写成math-ident="5em" $ \dfrac{1}{r}\frac{\partial }{{\partial r}}\left(r{u_r}\right) + \dfrac{{\partial {u_z}}}{{\partial {\rm{z}}}} = 0 $ (1) math-ident="5em" $ \dfrac{{\partial {u_r}}}{{\partial t}} + {u_r}\dfrac{{\partial {u_r}}}{{\partial r}} + {u_z}\dfrac{{\partial {u_r}}}{{\partial {\rm{z}}}} = - \dfrac{{\partial p}}{{\partial r}} + \dfrac{1}{{{Re} }}\left(\Delta {u_r} - \dfrac{{u_r}}{{{r^2}}}\right) $ (2) math-ident="5em" $ \frac{{\partial {u_z}}}{{\partial t}} + {u_r}\frac{{\partial {u_z}}}{{\partial r}} + {u_z}\frac{{\partial {u_z}}}{{\partial {\rm{z}}}} = - \frac{{\partial p}}{{\partial {\rm{z}}}} + \frac{1}{{{Re} }}\Delta {u_z} $ (3) 这里无量纲参数雷诺数
$ \mathit{Re}=UL/\nu $ , 其中$ U $ 和$ L $ 分别表示轴对称流动的特征速度和特征长度. 对于小球算例特征长度$ L=D $ ,$ D $ 为小球直径. 而对于圆盘特征长度$ L=D $ ,$ D $ 为圆盘径向直径. 另外需要注意的是轴对称柱坐标系下的拉普拉斯算符$ \mathrm{\Delta } $ 以及梯度算符$ \nabla $ 与直角坐标是不一样的, 分别为math-ident="5em" $ \Delta = \frac{1}{r}\frac{\partial }{{\partial r}}\left(r\frac{\partial }{{\partial r}}\right) + \frac{{{\partial ^2}}}{{\partial {{\textit{z}}^2}}} $ (4) math-ident="5em" $ \nabla = \left(\frac{\partial }{{\partial r}},\frac{\partial }{{\partial {\textit{z}}}}\right) $ (5) 另外, 对比式(2)和式(3), 我们能够清楚地发现在径向动量方程中由于坐标系的变化导致多出来的源项
$ {u_r}/{{r}^{2}} $ . 由于这一项在靠近对称轴位置变得很大, 因此对于该项的处理是十分重要的.因为本文的数值模拟计算也涉及到小球颗粒自由下落的运动过程. 所以除了上述流体的控制方程之外, 颗粒运动的控制方程也需要计算. 一般来说我们可以将颗粒的运动速度表示为
${\boldsymbol{U}}_{{\it{\Gamma}} }={\boldsymbol{U}}_{p} + {\boldsymbol{\omega }}_{p}\times {\boldsymbol{r}}_{{\it{\Gamma}} }$ , 其中${\boldsymbol{U}}_{p}$ 表示颗粒质心的平移速度、${\boldsymbol{\omega }}_{p}$ 表示颗粒的角速度、${\boldsymbol{r}}_{{\it{\Gamma}} }$ 表示颗粒表面到质心的距离矢量. 颗粒运动方程则分别是由牛顿动量方程(6)以及角动量方程(7)决定, 即math-ident="5em" $ \dfrac{{\partial \left({m_p}{{\boldsymbol{U}}_p}\right)}}{{\partial t}} = {{\boldsymbol{f}}_p} = {\rho _f}{{\boldsymbol{\oint }}{{\boldsymbol{\tau }} \cdot {\boldsymbol{n}}} _p}{\rm{d}}S + {V_p}\left({\rho _p} - {\rho _f}\right){\boldsymbol{g}} $ (6) math-ident="5em" $ \dfrac{{\partial \left({{\boldsymbol{I}}_p}{{\boldsymbol{\omega}} _p}\right)}}{{\partial t}} = {{\boldsymbol{T}}_p} = {\boldsymbol{\oint }}{{{\boldsymbol{r}}_p} \times {{\boldsymbol{f}}_p}} {\rm{d}}S $ (7) 这里
$\boldsymbol{\tau }=-\boldsymbol{I}p/{\rho }_{f} + \nu \left(\nabla \boldsymbol{u} + \nabla {\boldsymbol{u}}^{{\rm{T}}}\right)$ 表示应力张量, 其中$ {\boldsymbol{I}} $ 和$ p $ 分别为单位张量和流体动压(不含静水压力). 而下标$ p $ 和$ f $ 分别对应颗粒和流体.$ {V}_{p} $ ,$ {m}_{p} $ 和${\boldsymbol{n}}_{p}$ 则分别表示固体颗粒的体积、颗粒质量以及颗粒表面的单位外法向矢量.${\boldsymbol{f}}_{p}$ 为颗粒所受外力之和, 包含流体作用力、重力以及浮力.${\boldsymbol{I}}_{p}$ 和${\boldsymbol{r}}_{p}$ 分别表示固体颗粒的惯性张量以及颗粒所受到的外力作用点到颗粒重心的距离, 而${\boldsymbol{T}}_{p}$ 则表示颗粒所受到的合力矩. 因为本文只考虑无环流轴对称运动, 颗粒在自由下落运动过程中是没有转动发生的即固体颗粒的角速度${\boldsymbol{\omega }}_{p}={\boldsymbol{0}}$ . 这也意味着我们只需要单独求解颗粒的牛顿运动方程 (6) 即可.1.2 数值离散
本文的算法是基于柱坐标系以及有限体积法来开展的. 因此对于无环流轴对称流动, 子午面上生成两维的平面网格后, 将网格沿轴线(
${{z}}$ 轴)分别向正反环向方向旋转${\rm{d}}\theta/2$ 角度后便形成了一层楔形三维网格(夹角为$ {\rm{d}}\theta $ ). 如图2所示: 其中坐标$ r $ 轴位于网格中界面,${{S}}_{f}$ 表示网格面面积, nf表示网格面单位法向量.网格体积微元表示为
${{\it{\Omega}} }_{c}={\rm{d}}r{\rm{d}}{{z}}{r}_{c}{{\rm{d}}}{\theta }$ ,${{z}}$ 方向网格表面积计算表示为${S}_{{\textit{z}}}={r}_{f}{\rm{d}}\theta {\rm{d}}r$ ,$ r $ 方向网格表面积计算表示为${S}_{r}={r}_{f}{\rm{d}}\theta {\rm{d}}{\textit{z}}$ ,$ \theta $ 方向网格表面积计算表示为${S}_{r}={\rm{d}}\theta{\rm{d}}{\textit{z}}$ . 这里$ {r}_{c} $ 和$ {r}_{f} $ 分别为网格微元中心和网格面心径向坐标. 虽然在进行空间离散时, 由于有限体积法的使用引入了一层三维网格单元, 但在计算过程中本质上仍然是2D计算. 然后对流体控制方程进行的完整离散, 本文采用一个经典的具有二阶时间和空间精度的分步式投影法[23]. 主要计算过程可以分为以下四个步骤.第一步分别计算预测步径向速度
$ u_r $ 和轴向速度$u_z$ . 对于对流项以及黏性项中的拉普拉斯项均采用二阶中心差分格式离散, 而对于径向方程中多出来的${u_r}/{{r}^{2}}$ 项采用隐式处理. 则上述径向和轴向动量方程的离散形式为$$ \begin{split} \dfrac{{u_{rc}^* - u_{rc}^k}}{{\Delta t}}{\varOmega _c} + \dfrac{1}{2}\sum\limits_{}^{} {\left[m_f^k\left(u_{rf}^*\right) + m_f^k\left(u_{rf}^k\right)\right]} = - \left(\nabla {p_r}\right)_c^k{\varOmega _c} {\text{ }} + \\ \qquad\dfrac{1}{{2{Re} }}\sum\limits_{}^{} {\left[\left(\dfrac{{\partial u_r}}{{\partial n}}\right)_f^* + \left(\dfrac{{\partial u_r}}{{\partial n}}\right)_f^k\right]} {S_f} - \dfrac{{u_{rc}^*}}{{{Re} r_c^2}}{\varOmega _c} \end{split}$$ (8) $$\begin{split} \dfrac{{u_{zc}^* - u_{zc}^k}}{{\Delta t}}{\varOmega _c} &+ \dfrac{1}{2}\sum\limits_{}^{} {\left[m_f^k\left(u_{zf}^*\right) + m_f^k\left(u_{zf}^k\right)\right]} = - \left(\nabla {p_{\textit{z}}}\right)_c^k{\varOmega _c} {\text{ }} + \\ & \dfrac{1}{{2{Re} }}\sum\limits_{}^{} {\left[\left(\dfrac{{\partial u_z}}{{\partial n}}\right)_f^* + \left(\dfrac{{\partial u_z}}{{\partial n}}\right)_f^k\right]} {S_f} \end{split}$$ (9) 这里需要注意的是式(8)和式(9)中的上标
$ k $ 和$ * $ 分别表示当前时层(已知)和预测步时层(未知), 下标$ c $ 和$ f $ 则分别表示网格单元中心和网格单元面.${m}_{f}^{k}={{\boldsymbol{u}}}_{f}^{k}\cdot {\boldsymbol{n}}_{f}{{{S}}}_{f}$ 表示网格单元面上的质量通量.$ \nabla {p}_{r} $ 和$\nabla {p}_{{\it{{\textit{z}}}}}$ 则分别表示压力梯度的径向和轴向分量.${{\it{\Omega}} }_{c}$ 表示网格单元的体积(在环向$ \theta $ 方向单位长度取为1, 只有一层网格).第二步计算中间速度场
${\widetilde {\boldsymbol{u}}_c}=\left({\widetilde u_{rc}},{\widetilde u_{zc}}\right)$ . 这里我们利用已知$ k $ 时层的压力场$ {p}_{c}^{k} $ 以及第一步计算得到的预测步速度场${\boldsymbol{u}}_{c}^{*}=\left({u}_{rc}^{*},{u}_{zc}^{*}\right)$ 来计算中间速度$$ {\widetilde {\boldsymbol{u}}_c} = {\boldsymbol{u}}_c^* + \Delta t\left(\nabla p\right)_c^k$$ (10) 第三步计算
$ k + 1 $ 时层压力场$ {p}_{c}^{k + 1} $ . 首先根据上一步计算得到的网格中心的中间步速度场, 我们通过插值可以得到网格单元面上中间步质量通量$$ {\widetilde {\boldsymbol{m}}_f} = {\widetilde {\boldsymbol{u}}_{c \to f}} \cdot {{\boldsymbol{n}}_f}{S_f} $$ (11) 这里为了防止速度压力求解失耦, 发生棋盘压力震荡现象, 我们需要将网格单元中心的速度插值到网格单元面心上(
${\widetilde {\boldsymbol{u}}_{c \to f}}$ ). 其次通过上述的网格单元面上的中间步质量通量来求解压力泊松方程得到新时层的压力场. 这里压力泊松方程为$$ \sum\limits_{}^{} {\left(\dfrac{{\partial p}}{{\partial n}}\right)} _f^{k + 1}{S_f} = \dfrac{1}{{\Delta t}}\sum\limits_{}^{} {{{\widetilde m}_f}} $$ (12) 第四步根据
$ k + 1 $ 时层的压力场, 修正中间步速度场, 得到$ k + 1 $ 时层的网格单元中心和面上的速度场以及面上的质量通量分别为math-ident="5em" $ {\boldsymbol{u}}_c^{k + 1} = {\widetilde {\boldsymbol{u}}_c} - \Delta t\left(\nabla p\right)_c^{k + 1} $ (13) math-ident="5em" $ m_f^{k + 1} = {\widetilde m_f} - \Delta t\left(\frac{{\partial p}}{{\partial n}}\right)_f^{k + 1}{S_f}$ (14) 至此完整的一次投影步迭代计算结束. 需要说明的是对于网格中心处的梯度计算是采用高斯定理方式
$\nabla {p}_{c}=\dfrac{\sum\limits {p}_{f}{{{\boldsymbol{n}}}}_{f}{S}_{f}}{{\mathrm{\varOmega }}_{c}}$ , 而对于网格单元界面上的法向压力梯度是采用二阶中心差分形式${\left(\dfrac{\partial p}{\partial n}\right)}_{f}=\dfrac{{p}_{nb}-{p}_{A}}{{d}_{nb-A}}$ , 下标A和nb分别表示网格A与网格A的相邻网格nb,$ {d}_{nb-A} $ 表示网格A与相邻网格nb中心距离.另外对于在不可压缩黏性流体中颗粒小球垂直靠近壁面运动时, 由于运动颗粒和静止壁面之间的间隙很小, 数值模拟所需的网格分辨率非常高. 因此为了保证计算的稳定性, 需要在一个时间步长内重复多次执行投影步计算, 如图3所示.
根据本文测试的经验, 这里的迭代次数一般设置为3~5次即可. 另外在前两次迭代计算时, 一般会对压力和速度做松弛处理来保证计算的稳定性.
除了求解流体的Navier-Stokes方程之外, 本文也要求解浸没固体颗粒的牛顿运动方程. 这里我们直接采用Euler的显式离散格式处理, 如下式所示
$$ \dfrac{{{m_p}\left({\boldsymbol{U}}_p^{k + 1} - {\boldsymbol{U}}_p^k\right)}}{{\Delta t}} = {\rho _f}{{\boldsymbol{\oint}} {{{\boldsymbol{\tau}} ^{k + 1}} \cdot {\boldsymbol{n}}} _p}{\rm{d}}S + {V_p}\left({\rho _p} - {\rho _f}\right){\boldsymbol{g}} $$ (15) 这里的应力张量
${\boldsymbol{\tau }}^{k + 1} = -\boldsymbol{I}{p}_{ib}^{k + 1} + \mu \left[\nabla {\boldsymbol{u}}_{ib}^{k + 1} + {\left(\nabla {\boldsymbol{u}}_{ib}^{k + 1}\right)}^{{\rm{T}}}\right]$ , 其中浸没边界上的压力以及速度梯度分别为$ {p}_{ib}^{k + 1}={p}_{IB}^{k + 1} $ 和$\nabla {\boldsymbol{u}}_{ib}^{k + 1}=\dfrac{\left({\boldsymbol{u}}_{IB}^{k + 1}-{\boldsymbol{u}}_{ib}^{k}\right)}{{d}_{IB-ib}}{\boldsymbol{n}}_{ib}$ , 而$ {d}_{IB-ib} $ 表示IB网格单元到浸没边界上的映射ib点的距离. 需要注意的是${\boldsymbol{u}}_{ib}^{k}$ 可由当前浸没颗粒速度即${\boldsymbol{U}}_{p}^{k}$ 给出, 而$ {p}_{IB}^{k + 1} $ 和${\boldsymbol{u}}_{IB}^{k + 1}$ 需要通过由求解流体离散方程(8)~(14)得到的流体单元数据经过最小二乘法插值得到. 关于IB单元和映射ib点的位置可以参考图1. 在下一小节中我们会说明如何插值得到IB单元数据. 因此通过使用式(15), 我们很容易求解得到$ k + 1 $ 时层的颗粒速度. 进而更新$ k + 1 $ 时层的浸没边界条件, 继续求解下一时层流场.1.3 IB边界条件引入
因为本文采用离散力法, 所以固体浸没边界的引入主要分为三个步骤. 第一步通过最小二乘法插值得到IB网格单元的值, 第二步利用阶梯状锐利界面替代真实的固体浸没边界封闭控制方程求解流场信息, 第三步利用新求解的流场信息修正IB网格单元信息. 接下来我们将分别说明这三个步骤的实现.
第一步: 我们对IB单元中心物理量
$\varphi$ 在浸没边界投影点(图1$ ib $ 点)在子午面上做二维平面泰勒展开, 即$$ \begin{split} \varphi \left(x,y\right) \approx \varPhi \left( {x,y} \right) = {\varphi _{ib}} + \dfrac{{\partial \varphi }}{{\partial x}}{\Big|_{ib}}\Delta x + \dfrac{{\partial \varphi }}{{\partial y}}{\Big|_{ib}}\Delta y +\\\qquad \dfrac{1}{2}\dfrac{{{\partial ^2}\varphi }}{{\partial {x^2}}}{\Big|_{ib}}{\left(\Delta x\right)^2} + \dfrac{1}{2}\dfrac{{{\partial ^2}\varphi }}{{\partial {y^2}}}{\Big|_{ib}}{\left(\Delta y\right)^2} + \cdots \end{split} $$ (16) 其中投影点坐标用
$ ({x}_{ib},{y}_{ib}) $ ,$ \mathrm{\Delta }x=x-{x}_{ib} $ 和$ \mathrm{\Delta }y=y-{y}_{ib} $ 则分别表示IB单元中心到投影点的各个方向上的距离. 这里需要说明一下, 尽管本文流场的方程离散计算是在轴对称柱坐标系上进行的, 但是IB网格单元中心的插值还是在二维直角坐标基础上完成的. 那么只要确定了式(16)中的泰勒展开系数${\left. {\dfrac{{\partial \varphi }}{{\partial x}}} \right|_{ib}}$ ,${\left. {\dfrac{{\partial \varphi }}{{\partial y}}} \right|_{ib}}$ ,${\left. {\dfrac{{{\partial ^2}\varphi }}{{\partial {x^2}}}} \right|_{ib}}$ ,${\left. {\dfrac{{{\partial ^2}\varphi }}{{\partial {y^2}}}} \right|_{ib}}$ 等, 就能得到IB单元中心的值. 针对这些未知系数, 本文采用IB单元周围的已知流体网格单元以及壁面边界面元信息构建一个代数方程组(未知系数作为待求解变量). 通过最小二乘法求解方程系数, 以此确定式(16)中的泰勒展开系数. 当插值搜索的数据点多于插值方程未知系数的个数时, 我们通过带权重的最小二乘法[21]求解:$\gamma = {\displaystyle\sum\limits }_{m=1}^{M}{\omega }_{m}^{2}{\left[\varPhi \left({{\Delta}} x,{{\Delta}} y\right)-\varphi \left({{\Delta}} x,{{\Delta}} y\right)\right]}^{2}$ . 这里M和m分别表示数据点的总数和第m个数据点, 而${\omega }_{m} = \dfrac{1}{2}\left[1 + \mathit{\cos}\left(\dfrac{{\text{π}} {d}_{m}}{R}\right)\right]$ 表示距离权函数. 其中${d_m} = \sqrt {{{\Delta}} x_m^2 + {{\Delta}} y_m^2}$ 表示数据点与投影点之间的间距. R表示这些间距中的最大值.针对纽曼边界条件也是类似处理, 这里不再赘述, 详情可以参考文献[15]. 考虑到插值的精度和效率, 本文的插值精度采用三阶形式.
第二步: 阶梯近似界面引入, 如图1所示. 首先忽略浸没边界的存在, 则离散方程的最终离散形式在
$ k + 1 $ 时层的一般表达为$${A_p}\varphi _p^{k + 1} + \sum\limits_{fl = 1}^{nfl/(nIB)} {{{\left({A_{nb}}\varphi _{nb}^{k + 1}\right)}_{fl}}} + \sum\limits_{IB = 1}^{nIB} {{{\left({A_{nb}}\varphi _{nb}^{k + 1}\right)}_{IB}}} = S_p^k$$ (17) 其中下标
$ nb $ 表示网格单元$ p $ 周围的网格单元,$nfl/(nIB)$ 和$ nIB $ 则分别表示相邻网格单元中不包括IB单元的总数以及IB网格单元的数量. 这里方程(17)因为缺少边界条件是无法求解的. 但是我们可以利用第一步得到的IB单元中心的插值信息($ k $ 时层)以及阶梯近似边界来封闭流体离散方程. 此时新的离散方程最终离散形式为$$ {A_p}\varphi _p^{k + 1} + \sum\limits_{fl = 1}^{nfl/(nIB)} {{{\left({A_{nb}}\varphi _{nb}^{k + 1}\right)}_{fl}}} = S_p^k - \sum\limits_{IB = 1}^{nIB} {{{\left({A_{nb}}\varphi _{nb}^k\right)}_{IB}}} $$ (18) 这里IB单元被当做源项处理放在公式右侧.
第三步: 修正IB单元的信息. 通过第二步计算, 我们得到了
$ k + 1 $ 时层的流场信息. 但是IB单元还是$ k $ 时层的, 我们需要对其再做一次最小二乘法插值, 将其更新到$ k + 1 $ 时层. 至此IB边界条件引入完成.1.4 2D自适应网格技术实现
对于动边界问题, 由于颗粒在靠近壁面或与壁面发生碰撞弹跳之前需要运动较长距离, 而颗粒附近的黏性边界层和颗粒碰壁间隙处相对于远处流场需要较密的网格才能解析. 因此如果在颗粒运动路径都采用较密网格显然会造成很高的计算成本, 这显然是不能接受的. 而自适应网格技术很好地替我们解决这一困难, 因此本文开发了一套基于浸没边界法的2D自适应网格技术. 对于浸没边界分层加密策略, 主要通过设置三个加密参数来执行. 首先是最大加密次数k1, 其决定最小网格尺寸为
$ \mathrm{\Delta }h/{2}^{{k}_{1}} $ , 其中$ \mathrm{\Delta }h $ 为背景基础网格尺寸. 另外为了保证浸没边界附近的网格加密层过度平滑(有助于计算稳定), 我们针对在浸没边界附近的加密层, 设置前k2个加密层中每一个加密层具有N层网格. 一个$ {k}_{1}=3, {k}_{2}=2, N=8 $ 的加密结果如图4所示.图4给出的流场子午面上的加密网格. 粗网格和细网格加密界面上的法向和网格中心连线存在一定倾角. 针对这一现象, 本文采用倾斜校正处理使得计算更加稳定.
2. 数值算法验证
为了验证本文所提出的轴对称IBM算法的有效性与可靠性. 本章节将分别给出小球绕流、圆盘绕流、Stokes流小球近壁运动以及小球自由下落碰壁弹跳等四种轴对称流动的数值模拟验证结果.
2.1 静止小球绕流验证
根据Wu[24]的试验研究以及Johnson[16]和Tomboulides[25]数值模拟结果, 当雷诺数
$\mathit{Re}={U}_{\infty }D/ {\nu }_{f}$ 小于210时, 均匀来流通过一个固定小球所形成的流动是保持轴对称的. 这里$ {U}_{\infty } $ ,$ D $ ,$ {\nu }_{f} $ 分别表示远场均匀来流速度、静止小球直径以及流体的运动黏度. 因此本文分别测试了雷诺数为1, 10, 50, 100, 150, 200等六种不同工况, 并将本文的数值模拟结果与前人结果进行比较. 首先我们给出计算模型如图5 所示. 模型流向总长为$ 40 D $ , 小球中心距离入口和出口截面距离均为$ 20 D $ , 小球中心距离上侧壁为$ 20 D $ . 其次本文对于边界条件设置: 进口采用均匀来流边界条件$ {U}_{\infty } $ , 出口设置为对流出口边界条件, 对称轴设置为对称边界条件, 侧壁设置为远场边界条件, 浸没小球表面设置为无滑移边界条件. 模拟过程中本文采用自适应时间步长, 保证最大的库朗数(CFL)不超过0.5.这里考虑到本文的工况最高雷诺数为200以及Pan[26]的数值经验, 一个网格无关性验证如表1所示.
表 1 雷诺数Re = 200网格无关性验证Table 1. Grid independence test at Re = 200Case $D/\Delta x$ $ {C_d} $ $ {L_{re}} $ coarse mesh 160 0.768 1.435 fine mesh 320 0.770 1.430 表中
$ \mathrm{\Delta } $ x为网格最小尺寸,${C}_{d}={F}_{{\textit{z}}}/\left(\dfrac{1}{8}\rho {U}_{\infty }^{2}{\text{π}} {D}^{2}\right)$ 为阻力系数($ {F}_{z} $ 为流向阻力),$ {L}_{re} $ 表示小球尾流回流区长度. 通过粗细网格阻力系数和回流区长度结果对比, 我们发现粗网格对于本问题是足够解析的. 因此后续工况的数值结果均是基于粗网格实现. 在图6中, 从上到下依次给出雷诺数$ {Re} $ 为50, 100, 150, 200时采用轴对称IBM计算得到的流线结果, 白色半圆表示浸没小球. 从图中我们可以看到四种雷诺数下, 随着雷诺数的增加回流区长度也是增加的.图7(a)和图7(b)中, 将采用本文算法计算得到的回流区长度
$ {L}_{re} $ 以及阻力系数$ {C}_{d} $ 与Johnson[16]和Magnaudet[27]的结果进行定量对比. 根据图7(a)可以发现: 我们计算的回流区长度在雷诺数为50和100时与前人结果均符合较好, 在雷诺数为150和200时与Johnson[16]结果也符合很好. 而根据图7(b)可以发现: 我们通过轴对称IBM算法计算的阻力系数在低雷诺数和高雷诺数与Johnson[16]和Magnaudet[27]的均符合很好. 综上所述, 本文的轴对称浸没边界算法对于求解小球绕流中轴对称流动是可靠的、有效的.![]()
图 7 小球回流区长度和阻力系数对比Figure 7. Comparisons of the length of recirculation zone$ {L}_{re} $ and drag coefficient$ {C}_{d} $ of sphere2.2 静止圆盘绕流验证
除了小球这种特殊几何结构之外, 圆盘也是一种典型的轴对称结构. 为了更好地说明本文的算法对于固体几何的适应性, 我们对于低雷诺数下的轴对称圆盘绕流也进行了相应的测试. 根据文献[28-29]的数值研究: 对于纵横比
${\chi} ={D}/{H}_{{\textit{z}}}=10$ 的圆盘, 当雷诺数$ \mathit{Re}={U}_{\infty }D/{\nu }_{f} $ 小于135时, 远场均匀来流通过固定圆盘后形成稳态轴对称流动. 这里$ {U}_{\infty } $ ,$ D $ ,$ {H}_{z} $ ,$ {\nu }_{f} $ 分别表示远场均匀来流速度、静止圆盘径向直径、静止圆盘轴向厚度以及流体的运动黏度. 本文分别计算了雷诺数为10, 20 , 40, 60, 80, 100等六种工况. 为了更好地与文献[28]研究对比, 我们在计算模型的设置上采取与之相一致的参数. 计算模型如图8所示: 进口距离圆盘$ 2.5 D $ , 出口距离圆盘$ 15 D $ , 侧壁距离对称轴$ 6 D $ .这里关于边界条件的设置与上一小节的小球绕流中的设置是一样, 便不再复述. 我们同样采用两种不同分辨率的网格进行计算, 用以验证网格无关性, 如表2所示.
表 2 雷诺数Re = 100网格无关性验证Table 2. Grid independence test at Re = 100Case $D/\Delta x$ $ {C_d} $ $ {L_{re}} $ coarse mesh 160 1.327 1.960 fine mesh 320 1.329 1.965 根据表2 中的阻力系数以及回流区长度对比, 我们发现粗网格的数值计算结果已经足够求解该问题, 所以后续的低雷诺数下的工况我们均采用粗网格进行计算.
图9 中, 我们分别给出了雷诺数为10和100两种工况下的流场流线图, 白色矩形表示浸没圆盘. 我们可以直观地看出回流区长度也是随着雷诺数增加而增加的.
除此之外, 我们同样对比了圆盘尾迹回流区长度
$ {L}_{re} $ 以及阻力系数$ {C}_{d} $ , 对比结果如图10(a)和图10(b)所示. 我们发现二者在不同雷诺数时均与文献[28]符合较好, 并且阻力系数$ {C}_{d} $ 随着雷诺数的增加是减小的, 这与小球绕流结果是相似的. 综上所述, 本文的轴对称IBM算法对于圆盘几何的轴对称绕流问题求解结果也是有效和可靠的, 同时也说明本文算法对复杂的轴对称几何也是适用的.![]()
图 10 圆盘回流区长度和阻力系数对比Figure 10. Comparisons of the length of recirculation zone$ {L}_{re} $ and drag coefficient$ {C}_{d} $ of circular disk2.3 Stokes流小球近壁运动验证
在前面两个小节中, 本文对于静止浸没边界问题已经做了详细的验证和讨论, 充分说明了本文的轴对称IBM算法对于各种复杂几何静边界轴对称流动问题是能够正确求解的. 但是工业工程更多的是动边界问题, 本小节将通过数值模拟的方法计算匀速小球靠近壁面运动所受阻力, 并与经典的Brenner[30]小球碰壁润滑力解析解式(19)和式(20)进行对比
math-ident="5em" $ {F_{{lub} }} = 3{\text{π}} \rho \nu DU\lambda \left(\varepsilon \right) $ (19) math-ident="5em" $ \begin{split} \lambda \left(\varepsilon \right) = \dfrac{4}{3} & {\rm{sinh}} \alpha \sum\limits_1^\infty {\dfrac{{n\left(n + 1\right)}}{{\left(2n - 1\right)\left(2n + 1\right)}}} \cdot \\& \left\{\dfrac{{2{\rm{sinh}}\left[\left(2n + 1\right)\alpha \right]+ \left(2n + 1\right){\rm{sinh}}(2\alpha) }}{{4{\rm{sin}}{{\rm{h}}^2}\left[\left(n + 1/2\right)\alpha \right] - {{\left(2n + 1\right)}^2}{\rm{sin}}{{\rm{h}}^2}\alpha }} - 1\right\} \end{split} $ (20) 其中
$\varepsilon =\dfrac{H-1/(2D)}{1/(2D)}$ 表示无量纲球壁之间的间隙,$ H $ 为小球中心到壁面的距离,$ D $ 为小球直径. 该润滑力实质上为一个光滑刚体小球在静止不可压缩黏性流体中以恒定速度垂直向静止壁面靠近时所受阻力. Brenner[30]进行理论求解时假设流动是一种蠕动流动, 即不考虑流体的惯性. 但在数值模拟中我们还是要考虑惯性项, 不过可以通过将雷诺数设置很小达到理论近似下的蠕动流动状态. 这里我们采用$\mathit{Re}=\dfrac{UD}{{\nu }_{f}}=1$ 的设置. 相应的计算模型如图11 所示, 小球初始$ H $ 设置为$ 2 D $ . 对于上壁面和侧壁我们使用远场边界条件, 而对于下壁面和浸没边界均采用无滑移边界条件, 对称轴采用对称边界条件.![]()
图 11 小球匀速靠近壁面几何参数Figure 11. Configuration for a sphere with a constant velocity approaching a wall图12给出了不同间隙时的流场压力云图. 从图中我们可以看出由于小球不断靠近壁面, 间隙流体被挤压向径向外侧流动的同时, 间隙处的压力会急剧升高导致小球阻力上升. 当间隙特别小时, 压力会变得非常大, 相应所需要的网格分辨率也就会很高. 本文阻力数值模拟结果与Brenner[30]理论结果对比如图13 所示: 采用了800, 1600和3200三种不同分辨率的网格进行计算, 很明显三种网格的数值模拟结果在
$\varepsilon > 0.1$ 时和理论结果都能符合很好, 但是在小球无限接近壁面时, 即$\varepsilon \to 0$ 时结果开始出现偏差, 而且网格分辨率越小偏差越大.![]()
图 12 压力云图间隙为(a)$ \varepsilon = 1.1 $ 和(b)$ \varepsilon = 0.1 $ Figure 12. Contour of pressure around sphere with (a)$\varepsilon =1.1$ and (b)$\varepsilon =0.1$ 这是因为随着间隙越来越小, 压力场梯度剧烈变化, 想要精确地解析流场, 所需要的网格分辨率是非常高的. 虽然这种极限小间隙情况很难解析正确, 但是在网格分辨率足够的情况下, 我们的结果依然是有效的. 这可以说明本文的轴对称算法对于给定速度的强迫运动产生的动边界问题也是可以正确求解的.
2.4 小球自由下落碰壁弹跳验证
上一小节中, 本文验证了强迫运动的动边界问题, 本节我们通过耦合颗粒运动学方程以及流体控制方程来数值模拟小球颗粒在不可压缩黏性流体中自由下落碰撞壁面过程. 针对这样的一个球壁碰撞运动, Gondret等[31]通过实验研究, 给出了小球碰撞后的运动轨迹以及相应的恢复系数. 文献[32-33]通过数值模拟的手段再现了Gondret等[31]的实验结果. 为了验证本文轴对称IBM算法的正确性, 我们选取与之相同的工况参数设置并进行结果对比. 由于本节问题和上一小节类似, 上一小节的计算模型在本节仍然继续使用. 在进行碰撞数值模拟之前, 我们需要引入两个理论模型: 一个是润滑力模型, 另一个是碰撞模型. 首先对于润滑力模型, 结合上一小节, 我们发现间隙越来越小所需网格分辨率越来越高. 直接模拟碰撞将会耗费巨大的资源, 而且还不能求解准确, 因此在小间隙时一个润滑力修正模型[34-35]的引入是有必要的. 再加入润滑力模型后, 颗粒运动方程(6)修正为式(21), 其中润滑力由式(22)[35]和式(23)[36]给出(式(23)为式(20)的近似形式). 其中
$ {U}_{n} $ 为颗粒垂直壁面法向(z轴方向)速度大小,${\varepsilon }_{al}$ 为激活润滑力模型间隙条件. 由于本文的模拟参数$ \mathit{Re} < 210 $ , 流场保持为轴对称状态, 因此小球速度${\boldsymbol{U}}_{p}$ 只有垂直壁面法向分量$ {U}_{n} $ math-ident="5em" $ \begin{split} \dfrac{{\partial ({m_p}{{\boldsymbol{U}}_p})}}{{\partial t}} =& {\rho _f}{{\boldsymbol{\oint}} {{\boldsymbol{\tau}} \cdot {\boldsymbol{n}}} _p}{\rm{d}}S + {V_p}\left({\rho _p} - {\rho _f}\right){\boldsymbol{g}} {\text{ }} + \\& {F_{{lub} correct}}{{ {\boldsymbol{e}}}_{\textit{z}}} {\text{ = }}{{\boldsymbol{F}}_h} + {F_{{lub} correct}}{{ {\boldsymbol{e}}}_{\textit{z}}} \end{split} $ (21) math-ident="5em" $ {F_{{lub} correct}} = - 3{\text{π}} \rho \nu D{U_n}\left[ {\lambda \left(\varepsilon \right) - \lambda \left({\varepsilon _{al}}\right)} \right] $ (22) math-ident="5em" $ \lambda \left(\varepsilon \right) = \dfrac{1}{\varepsilon } - \dfrac{1}{5}{\rm{ln}}\varepsilon - \dfrac{1}{{21}}\varepsilon {\rm{ln}}\varepsilon + 0.971\;3 $ (23) 这里因为
${1}/{\varepsilon }$ 的存在, 导致当间隙趋近于零时校正润滑力趋近于无穷大. 这显然不符合物理实际, 事实上由于颗粒表面粗糙度的存在, 一般当$\varepsilon$ 接近于表面粗糙度时碰撞已经发生($\varepsilon ={\varepsilon }_{w}$ ). 此时润滑力不会继续增大, 颗粒进入固体碰撞阶段. 这一阶段流体对于颗粒小球的作用力相比于固体的碰撞力可以忽略. 为了解析这样的碰撞过程, 使用一个经典的软球模型[34], 即$$ \frac{{\partial ({m_p}{U_n})}}{{\partial t}} = - {k_n}{\delta _{s - w}} - {\eta _n}{U_n} $$ (24) 其中
$ {k}_{n} $ 和$ {\eta }_{n} $ 分别为刚度系数和阻尼系数且由式(25)给出,$ {e}_{n,d} $ 表示空气中的干碰撞恢复系数与颗粒材料有关, 一般直接参考实验数据. 根据Gondret等[31]文献, 本文碰撞算例均采用$ {e}_{n,d}=0.98 $ 的设置.$ {\delta }_{s-w} $ 表示颗粒与壁面重叠位移. 碰撞时间步长${{\Delta}} {t}_{c}=\dfrac{1}{10}{{\Delta}} {t}_{f}$ [33],${{\Delta}} {t}_{f}$ 表示流体时间步长.$N{{\Delta}}{t}_{f}$ 为人为设定的碰撞时间, 本文中$ N=10 $ [32-34]. 实际上碰撞发生的时间是比较短的, 这里人为地放大了碰撞接触时间, 避免速度改变过于剧烈导致计算误差太大$$ {k_n} = \frac{{{m_p}({{\text{π}} ^2} + {\rm{l}}{{\rm{n}}^2}{e_{n,d}})}}{{{{(N\Delta {t_f})}^2}}},\;{\eta _n} = - \frac{{2{m_p}{\rm{ln}}{e_{n,d}}}}{{N\Delta {t_f}}} $$ (25) 完整小球与壁面正碰撞数值模拟过程如图14所示: 当间隙
$\varepsilon > {\varepsilon }_{al}$ 时, 颗粒只受到来自流体的作用力与重力; 当${\varepsilon }_{w} < \varepsilon < {\varepsilon }_{al}$ 时, 颗粒除了受到来自流体的作用力及重力之外, 还有一项润滑修正力; 当$\varepsilon < {\varepsilon }_{w}$ 时, 颗粒和壁面已经发生碰撞接触, 颗粒运动由方程(24)控制.这里我们对比了冲击壁面Stokes数为27和152两种工况结果, 具体的物理参数如表3所示. 其中Stokes数
$ S{t}_{c}={\rho }_{p}{U}_{in}D/\left(9{\rho }_{f}\nu \right) $ , 雷诺数$ \mathit{Re}={U}_{in}D/{\nu }_{f} $ ,$ {U}_{in} $ 表示小球自由下落后在壁面影响前达到的稳定速度.表 3 物理计算参数Table 3. Physical and computational parametersCase A B $ S{t_c} $ 27 152 $ {Re} $ 30 165 $ {U_{in}} $/(m·s−1) 0.518 0.585 $ D $/m 0.006 0.003 $ {\rho _p} $/(kg·m−3) 7800 7800 $ {\rho _f} $/(kg·m−3) 965 935 $ \nu $/(s·m−2) $1.036\;3 \times {10^{ - 4} }$ $1.069\;5 \times {10^{ - 5} }$ $ D/\Delta x $ 32 32 $ {\varepsilon _{al}} $ 0.5 0.5 $ {\varepsilon _w} $ 0.001 0.001 图15 (a)和15(b) 给出了两种工况下的数值模拟结果与Gondret等[31]实验结果对比. 图中
${H}^{*}=\dfrac{H-1/(2 D)}{1/(2 D)}$ 和${T}^{*}=\dfrac{(t-{t}_{c}){U}_{in}}{1/(2 D)}$ 表示无量纲高度和时间, 其中$ {t}_{c} $ 为首次碰撞发生时间. 我们可以发现小球的轨迹均符合较好. 因此本文的轴对称算法对于耦合颗粒运动学方程的动边界问题也是准确的.3. 结 论
本文针对颗粒两相流中的无环流轴对称流动问题, 发展了基于2D笛卡尔网格的轴对称IBM算法. 该IBM算法通过一个阶梯锐利状界面实现固体浸没边界引入, 对于界面上的IB网格单元通过最小二乘法插值得到. 本文通过柱坐标系建立流体控制方程, 并采用一个有限体积格式的经典投影算法进行方程离散. 同时对方程中由于柱坐标系的使用而产生的源项
${u_r}/{{r}^{2}}$ 采用隐式格式处理. 针对动边界, 本文也发展了一个2D自适应网格技术以提升计算效率. 通过小球绕流、圆盘绕流、Stokes流小球近壁运动以及小球自由下落碰壁弹跳等数值模拟算例, 验证该轴对称IBM算法对于复杂几何边界以及动边界轴对称流动的求解是高效的和准确的. -
![]()
图 1 浸没边界和网格单元划分示意图
Figure 1. Schematic diagram of immersed boundary and dividing cells type
![]()
图 7 小球回流区长度和阻力系数对比
Figure 7. Comparisons of the length of recirculation zone
$ {L}_{re} $ and drag coefficient$ {C}_{d} $ of sphere![]()
图 10 圆盘回流区长度和阻力系数对比
Figure 10. Comparisons of the length of recirculation zone
$ {L}_{re} $ and drag coefficient$ {C}_{d} $ of circular disk![]()
图 11 小球匀速靠近壁面几何参数
Figure 11. Configuration for a sphere with a constant velocity approaching a wall
![]()
图 12 压力云图间隙为(a)
$ \varepsilon = 1.1 $ 和(b)$ \varepsilon = 0.1 $ Figure 12. Contour of pressure around sphere with (a)
$\varepsilon =1.1$ and (b)$\varepsilon =0.1$ 表 1 雷诺数Re = 200网格无关性验证
Table 1 Grid independence test at Re = 200
Case $D/\Delta x$ $ {C_d} $ $ {L_{re}} $ coarse mesh 160 0.768 1.435 fine mesh 320 0.770 1.430 
下载: 导出CSV
表 2 雷诺数Re = 100网格无关性验证
Table 2 Grid independence test at Re = 100
Case $D/\Delta x$ $ {C_d} $ $ {L_{re}} $ coarse mesh 160 1.327 1.960 fine mesh 320 1.329 1.965
下载: 导出CSV
表 3 物理计算参数
Table 3 Physical and computational parameters
Case A B $ S{t_c} $ 27 152 $ {Re} $ 30 165 $ {U_{in}} $/(m·s−1) 0.518 0.585 $ D $/m 0.006 0.003 $ {\rho _p} $/(kg·m−3) 7800 7800 $ {\rho _f} $/(kg·m−3) 965 935 $ \nu $/(s·m−2) $1.036\;3 \times {10^{ - 4} }$ $1.069\;5 \times {10^{ - 5} }$ $ D/\Delta x $ 32 32 $ {\varepsilon _{al}} $ 0.5 0.5 $ {\varepsilon _w} $ 0.001 0.001
下载: 导出CSV
-
[1] 丛成华, 邓小刚, 毛枚良. 绕椭球的低速流动研究. 力学进展, 2021, 51(3): 467-619 (Cong Chenghua, Deng Xiaogang, Mao Meiliang. Advances in complex low speed flow around a prolate spheroid. Advances in Mechanics, 2021, 51(3): 467-619 (in Chinese) doi: 10.6052/1000-0992-20-036 Cong Chenghua, Deng Xiaogang, Mao Meiliang. Advances in complex low speed flow around a prolate spheroid. Advances in Mechanics, 2021, 51(3): 467-619(in Chinese) doi: 10.6052/1000-0992-20-036
[2] 吴介之, 刘罗勤, 刘天舒. 定常升阻力普适理论的特色和升力的物理来源. 力学进展, 2021, 51(1): 106-129 (Wu Jiezhi, Liu Luoqin, Liu Tianshu. The universal steady lift and drag theory and the physical origin of lift. Advances in Mechanics, 2021, 51(1): 106-129 (in Chinese) doi: 10.6052/1000-0992-20-014 Wu Jiezhi, Liu Luoqin, Liu Tianshu. The universal steady lift and drag theory and the physical origin of lift. Advances in Mechanics, 2021, 51(1): 106-129(in Chinese) doi: 10.6052/1000-0992-20-014
[3] Peskin CS. Flow patterns around heart valves: a digital computer method for solving the equations of motion. [PhD Thesis]. New York: Yeshiva University, 1972
[4] Peskin CS. The fluid dynamics of heart valves: Experimental, theoretical, and computational methods. Annual Review of Fluid Mechanics, 1982, 14(1): 235-259 doi: 10.1146/annurev.fl.14.010182.001315
[5] Iaccarino G, Verzicco R. Immersed boundary technique for turbulent flow simulations. Applied Mechanics Reviews, 2003, 56(3): 331-347 doi: 10.1115/1.1563627
[6] Yang J, Stern F. A simple and efficient direct forcing immersed boundary framework for fluid –structure interactions. Journal of Computational Physics, 2012, 231(15): 5029-5061 doi: 10.1016/j.jcp.2012.04.012
[7] Grigoriadis D, Kassinos S, Votyakov E, et al. Immersed boundary method for the MHD flows of liquid metals. Journal of Computational Physics, 2009, 228(3): 903-920 doi: 10.1016/j.jcp.2008.10.017
[8] 秦如冰, 柴翔, 程旭. 基于浸没边界法的流固耦合模拟分析. 核科学与工程, 2020, 40(5): 763-770 (Qin Rubing, Chai Xiang, Cheng Xu. Fluid-solid coupling simulation analysis based on immersed boundary method. Nuclear Science and Engineering, 2020, 40(5): 763-770 (in Chinese) Qin Rubing, Chai Xiang, Cheng Xu. Fluid-solid Coupling Simulation Analysis Based on Immersed Boundary Method. Nuclear Science and Engineering. 2020, 40(5): 763-770(in Chinese)
[9] 周睿, 程永光, 吴家阳. 用浸没边界-格子 Boltzmann 方法模拟双层刚性植被明渠水流特性. 水动力学研究与进展(A辑), 2019, 34(4): 503-511 (Zhou Rui, Cheng Yongguang, Wu Jiayang. Simulation of the open channel flow over double-layered rigid submerged vegetation by the immersed boundary-lattice Boltzmann method. Hydro-Science and Engineering (A) , 2019, 34(4): 503-511 (in Chinese) Zhou Rui, Cheng Yongguang, Wu Jiayang. Simulation of the open channel flow over double-layered rigid submerged vegetation by the immersed boundary-lattice Boltzmann method. Hydro-Science and Engineering (A), 2019, 34(04): 503-511(in Chinese)
[10] 周睿, 程永光, 吴家阳. 基于浸没边界法的水库变动水面模拟及验证. 水利水运工程学报, 2020, 1: 66-73 (Zhou Rui, Cheng Yongguang, Wu Jiayang. Simulation and verification of reservoir fluctuating water level based on immersed boundary method. Hydro-Science and Engineering, 2020, 1: 66-73 (in Chinese) doi: 10.12170/20180821001 Zhou Rui, Cheng Yongguang, Wu Jiayang. Simulation and verification of reservoir fluctuating water level based on immersed boundary method. Hydro-Science and Engineering, 2020,1: 66-73(in Chinese) doi: 10.12170/20180821001
[11] 佟莹, 夏健, 陈龙等. 基于隐式扩散的直接力格式浸没边界格子Boltzmann 方法. 力学学报, 2022, 54(1): 94-105 (Tong Ying, Xia Jian, Chen Long, et al. An immersed boundary lattice Boltzmann method based on implicit diffuse directforcing scheme. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 94-105 (in Chinese) doi: 10.6052/0459-1879-21-315 Tong Ying, Xia Jian, Chen Long, Xue Haotian. An Immersed boundary lattice Boltzmann method based on implicit diffuse directforcing scheme. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 94-105(in Chinese) doi: 10.6052/0459-1879-21-315
[12] 李桥忠, 陈木凤, 李游等. 浸没边界–简化热格子 Boltzmann 方法研究及其应用. 力学学报, 2019, 51(2): 392-404 (Li Qiaozhong, Chen Mufeng, Li You, et al. Immersed boundary-simplified thermal lattice Boltzmann method for fluid-structure interaction problem with heat transfer and its application. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(2): 392-404 (in Chinese) Li Qiaozhong, Chen Mufeng, Li You, Niu Xiaodong, Adnan Khan. Immersed boundary-simplified thermal lattice Boltzmann method for fluid-structure interaction problem with heat transfer and its application. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(2): 392-404 (in Chinese))
[13] Mittal R, Iaccarino G. Immersed boundary methods. Annual Review of Fluid Mechanics, 2005, 37: 239-261 doi: 10.1146/annurev.fluid.37.061903.175743
[14] Seo JH, Mittal R. A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations. Journal of Computational Physics, 2011, 230(19): 7347-7363 doi: 10.1016/j.jcp.2011.06.003
[15] Pan JH, Ni MJ, Zhang NM. A consistent and conservative immersed boundary method for MHD flows and moving boundary problems. Journal of Computational Physics, 2018, 373: 425-445 doi: 10.1016/j.jcp.2017.12.034
[16] Johnson TA, Patel VC. Flow past a sphere up to a Reynolds number of 300. Journal of Fluid Mechanics, 2000, 378: 19-70
[17] Jenny M, Dušek J, Bouchet G. Instabilities and transition of a sphere falling or ascending freely in a Newtonian fluid. Journal of Fluid Mechanics, 2004, 508: 201-239 doi: 10.1017/S0022112004009164
[18] Lai MC, Huang C, Huang YM. Simulating the axisymmetric interfacial flows with insoluble surfactant by immersed boundary method. International Journal of Numerical Analysis & Modeling, 2011, 8(1): 105-117
[19] Li Y, Yun A, Kim J. An immersed boundary method for simulating a single axisymmetric cell growth and division. Journal of Mathematical Biology, 2012, 65(4): 653-675 doi: 10.1007/s00285-011-0476-7
[20] Li X. An experimental and numerical study of normal particle collisions in a viscous liquid. [PhD Thesis]. California: California Institute of Technology, 2010
[21] Seo JH, Mittal R. A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries. Journal of Computational Physics, 2011, 230(4): 1000-1019 doi: 10.1016/j.jcp.2010.10.017
[22] Robertson I, Li L, Sherwin SJ, et al. A numerical study of rotational and transverse galloping rectangular bodies. Journal of Fluids and Structures, 2003, 17(5): 681-699 doi: 10.1016/S0889-9746(03)00008-2
[23] Ni MJ, Munipalli R, Huang P, et al. A current density conservative scheme for incompressible MHD flows at a low magnetic Reynolds number. Part II: On an arbitrary collocated mesh. Journal of Computational Physics, 2007, 227(1): 205-228
[24] Wu JS, Faeth GM. Sphere wakes in still surroundings at intermediate Reynolds numbers. AIAA Journal, 2012, 31(8): 1448-1455
[25] Tomboulides AG, Orszag SA. Numerical investigation of transitional and weak turbulent flow past a sphere. Journal of Fluid Mechanics, 2000, 416: 45-73 doi: 10.1017/S0022112000008880
[26] Pan JH, Zhang NM, Ni MJ. The wake structure and transition process of a flow past a sphere affected by a streamwise magnetic field. Journal of Fluid Mechanics, 2018, 842: 248-272 doi: 10.1017/jfm.2018.133
[27] Magnaudet J, Rivero M, Fabre J. Accelerated flows past a rigid sphere or a spherical bubble. Part 1. Steady straining flow. Journal of Fluid Mechanics, 1995, 284: 97-135
[28] Shenoy AR, Kleinstreuer C. Flow over a thin circular disk at low to moderate Reynolds numbers. Journal of Fluid Mechanics, 2008, 605: 253-262 doi: 10.1017/S0022112008001626
[29] Fernandes PC, Risso F, Ern P, et al. Oscillatory motion and wake instability of freely rising axisymmetric bodies. Journal of Fluid Mechanics, 2007, 573: 479-502 doi: 10.1017/S0022112006003685
[30] Brenner H. The slow motion of a sphere through a viscous fluid towards a plane surface. Chemical Engineering Science, 1961, 16(3-4): 242-251 doi: 10.1016/0009-2509(61)80035-3
[31] Gondret P, Lance M, Petit L. Bouncing motion of spherical particles in fluids. Physics of Fluids, 2002, 14(2): 643-652 doi: 10.1063/1.1427920
[32] Kempe T, Fröhlich J. Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids. Journal of Fluid Mechanics, 2012, 709: 445-489 doi: 10.1017/jfm.2012.343
[33] Xia Y, Xiong H, Yu Z, et al. Effects of the collision model in interface-resolved simulations of particle-laden turbulent channel flows. Physics of Fluids, 2020, 32(10): 103303 doi: 10.1063/5.0020995
[34] Costa P, Boersma BJ, Westerweel J, et al. Collision model for fully resolved simulations of flows laden with finite-size particles. Physical Review E, 2015, 92(5): 053012 doi: 10.1103/PhysRevE.92.053012
[35] Brändle JC, Breugem WP, Gazanion B, et al. Numerical modelling of finite-size particle collisions in a viscous fluid. Physics of Fluids, 2013, 25(8): 083302 doi: 10.1063/1.4817382
[36] Jeffrey DJ. Low-Reynolds-number flow between converging spheres. Mathematika, 1982, 29(1): 58-66 doi: 10.1112/S002557930001216X
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 694
- HTML全文浏览量: 254
- PDF下载量: 104
- 被引次数: 1
