基于2D 网格的轴对称浸没边界法1)

蔡政刚 潘君华 倪明玖

(中国科学院大学工程科学学院,北京 100049)

引言

流体领域中的颗粒两相流[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 浸没边界和网格单元划分示意图Fig.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 控制方程

本文算法考虑的是一种不可压缩的、运动黏性为 ν、密度为 ρf的流体所产生的无环流轴对称流动.在轴对称假设之下,本文可以利用柱坐标系将三维流动的控制方程简化成两维形式.本文使用u=(ur,uz)以及p来分别表示子午面上的速度和压力,其中ur和uz分别表示径向和轴向的速度分量.因此无量纲的流体动力学控制方程可以写成

这里无量纲参数雷诺数Re=UL/ν,其中U和L分别表示轴对称流动的特征速度和特征长度.对于小球算例特征长度L=D,D为小球直径.而对于圆盘特征长度L=D,D为圆盘径向直径.另外需要注意的是轴对称柱坐标系下的拉普拉斯算符 Δ 以及梯度算符 ∇ 与直角坐标是不一样的,分别为

另外,对比式(2)和式(3),我们能够清楚地发现在径向动量方程中由于坐标系的变化导致多出来的源项ur/r2.由于这一项在靠近对称轴位置变得很大,因此对于该项的处理是十分重要的.

因为本文的数值模拟计算也涉及到小球颗粒自由下落的运动过程.所以除了上述流体的控制方程之外,颗粒运动的控制方程也需要计算.一般来说我们可以将颗粒的运动速度表示为UΓ=Up+ωp×rΓ,其中Up表示颗粒质心的平移速度、ωp表示颗粒的角速度、rΓ表示颗粒表面到质心的距离矢量.颗粒运动方程则分别是由牛顿动量方程(6)以及角动量方程(7)决定,即

这里 τ=-Ip/ρf+ν(∇u+∇uT)表示应力张量,其中I和p分别为单位张量和流体动压(不含静水压力).而下标p和f分别对应颗粒和流体.Vp,mp和np则分别表示固体颗粒的体积、颗粒质量以及颗粒表面的单位外法向矢量.fp为颗粒所受外力之和,包含流体作用力、重力以及浮力.Ip和rp分别表示固体颗粒的惯性张量以及颗粒所受到的外力作用点到颗粒重心的距离,而Tp则表示颗粒所受到的合力矩.因为本文只考虑无环流轴对称运动,颗粒在自由下落运动过程中是没有转动发生的即固体颗粒的角速度 ωp=0 .这也意味着我们只需要单独求解颗粒的牛顿运动方程 (6)即可.

1.2 数值离散

本文的算法是基于柱坐标系以及有限体积法来开展的.因此对于无环流轴对称流动,子午面上生成两维的平面网格后,将网格沿轴线(z轴)分别向正反环向方向旋转 dθ/2 角度后便形成了一层楔形三维网格(夹角为 dθ).如图2 所示:其中坐标r轴位于网格中界面,S f表示网格面面积,nf表示网格面单位法向量.

图2 柱坐标系下的体积微元Fig.2 Volume element in cylindrical coordinates

网格体积微元表示为 Ωc=drdzrcdθ,z方向网格表面积计算表示为Sz=rfdθdr,r方向网格表面积计算表示为Sr=rfdθdz,θ 方向网格表面积计算表示为Sr=dθdz.这里rc和rf分别为网格微元中心和网格面心径向坐标.虽然在进行空间离散时,由于有限体积法的使用引入了一层三维网格单元,但在计算过程中本质上仍然是2D 计算.然后对流体控制方程进行的完整离散,本文采用一个经典的具有二阶时间和空间精度的分步式投影法[23].主要计算过程可以分为以下四个步骤.

第一步分别计算预测步径向速度ur和轴向速度uz.对于对流项以及黏性项中的拉普拉斯项均采用二阶中心差分格式离散,而对于径向方程中多出来的ur/r2项采用隐式处理.则上述径向和轴向动量方程的离散形式为

这里需要注意的是式(8)和式(9)中的上标k和 * 分别表示当前时层(已知)和预测步时层(未知),下标c和f则分别表示网格单元中心和网格单元面.表示网格单元面上的质量通量.∇pr和∇pz则分别表示压力梯度的径向和轴向分量.Ωc表示网格单元的体积(在环向 θ 方向单位长度取为1,只有一层网格).

第三步计算k+1时层压力场.首先根据上一步计算得到的网格中心的中间步速度场,我们通过插值可以得到网格单元面上中间步质量通量

这里为了防止速度压力求解失耦,发生棋盘压力震荡现象,我们需要将网格单元中心的速度插值到网格单元面心上().其次通过上述的网格单元面上的中间步质量通量来求解压力泊松方程得到新时层的压力场.这里压力泊松方程为

第四步根据k+1 时层的压力场,修正中间步速度场,得到k+1 时层的网格单元中心和面上的速度场以及面上的质量通量分别为

至此完整的一次投影步迭代计算结束.需要说明的是对于网格中心处的梯度计算是采用高斯定理方式 ∇pc=,而对于网格单元界面上的法向压力梯度是采用二阶中心差分形式,下标A和nb分别表示网格A与网格A的相邻网格nb,dnb-A表示网格A与相邻网格nb中心距离.

另外对于在不可压缩黏性流体中颗粒小球垂直靠近壁面运动时,由于运动颗粒和静止壁面之间的间隙很小,数值模拟所需的网格分辨率非常高.因此为了保证计算的稳定性,需要在一个时间步长内重复多次执行投影步计算,如图3 所示.

图3 流体方程计算流程图Fig.3 Flow chart of fluid equation calculation

根据本文测试的经验,这里的迭代次数一般设置为3～5 次即可.另外在前两次迭代计算时,一般会对压力和速度做松弛处理来保证计算的稳定性.

除了求解流体的Navier-Stokes 方程之外,本文也要求解浸没固体颗粒的牛顿运动方程.这里我们直接采用Euler 的显式离散格式处理,如下式所示

这里的应力张量 τk+1=,其中浸没边界上的压力以及速度梯度分别为,而dIB-ib表示IB 网格单元到浸没边界上的映射ib点的距离.需要注意的是可由当前浸没颗粒速度即给出,而需要通过由求解流体离散方程(8)～(14)得到的流体单元数据经过最小二乘法插值得到.关于IB 单元和映射ib点的位置可以参考图1.在下一小节中我们会说明如何插值得到IB 单元数据.因此通过使用式(15),我们很容易求解得到k+1 时层的颗粒速度.进而更新k+1 时层的浸没边界条件,继续求解下一时层流场.

1.3 IB 边界条件引入

因为本文采用离散力法,所以固体浸没边界的引入主要分为三个步骤.第一步通过最小二乘法插值得到IB 网格单元的值,第二步利用阶梯状锐利界面替代真实的固体浸没边界封闭控制方程求解流场信息,第三步利用新求解的流场信息修正IB 网格单元信息.接下来我们将分别说明这三个步骤的实现.

第一步:我们对IB 单元中心物理量 φ 在浸没边界投影点(图1ib点)在子午面上做二维平面泰勒展开,即

其中投影点坐标用 (xib,yib),Δx=x-xib和Δy=y-yib则分别表示IB 单元中心到投影点的各个方向上的距离.这里需要说明一下,尽管本文流场的方程离散计算是在轴对称柱坐标系上进行的,但是IB 网格单元中心的插值还是在二维直角坐标基础上完成的.那么只要确定了式(16)中的泰勒展开系数,等,就能得到IB 单元中心的值.针对这些未知系数,本文采用IB 单元周围的已知流体网格单元以及壁面边界面元信息构建一个代数方程组(未知系数作为待求解变量).通过最小二乘法求解方程系数,以此确定式(16)中的泰勒展开系数.当插值搜索的数据点多于插值方程未知系数的个数时,我们通过带权重的最小二乘法[21]求解:γ=.这里M和m分别表示数据点的总数和第m个数据点,而表示距离权函数.其中表示数据点与投影点之间的间距.R表示这些间距中的最大值.

针对纽曼边界条件也是类似处理,这里不再赘述,详情可以参考文献[15].考虑到插值的精度和效率,本文的插值精度采用三阶形式.

第二步:阶梯近似界面引入,如图1 所示.首先忽略浸没边界的存在,则离散方程的最终离散形式在k+1 时层的一般表达为

其中下标nb表示网格单元p周围的网格单元,nfl/(nIB)和nIB则分别表示相邻网格单元中不包括IB 单元的总数以及IB 网格单元的数量.这里方程(17)因为缺少边界条件是无法求解的.但是我们可以利用第一步得到的IB 单元中心的插值信息(k时层)以及阶梯近似边界来封闭流体离散方程.此时新的离散方程最终离散形式为

这里IB 单元被当做源项处理放在公式右侧.

第三步:修正IB 单元的信息.通过第二步计算,我们得到了k+1 时层的流场信息.但是IB 单元还是k时层的,我们需要对其再做一次最小二乘法插值,将其更新到k+1 时层.至此IB 边界条件引入完成.

1.4 2D 自适应网格技术实现

对于动边界问题,由于颗粒在靠近壁面或与壁面发生碰撞弹跳之前需要运动较长距离,而颗粒附近的黏性边界层和颗粒碰壁间隙处相对于远处流场需要较密的网格才能解析.因此如果在颗粒运动路径都采用较密网格显然会造成很高的计算成本,这显然是不能接受的.而自适应网格技术很好地替我们解决这一困难,因此本文开发了一套基于浸没边界法的2D 自适应网格技术.对于浸没边界分层加密策略,主要通过设置三个加密参数来执行.首先是最大加密次数k1,其决定最小网格尺寸为,其中Δh为背景基础网格尺寸.另外为了保证浸没边界附近的网格加密层过度平滑(有助于计算稳定),我们针对在浸没边界附近的加密层,设置前k2个加密层中每一个加密层具有N层网格.一个k1=3,k2=2,N=8的加密结果如图4 所示.

图4 给出的流场子午面上的加密网格.粗网格和细网格加密界面上的法向和网格中心连线存在一定倾角.针对这一现象,本文采用倾斜校正处理使得计算更加稳定.

图4 2D 自适应网格Fig.4 An 2D AMR grid around immersed boundary

2 数值算法验证

为了验证本文所提出的轴对称IBM 算法的有效性与可靠性.本章节将分别给出小球绕流、圆盘绕流、Stokes 流小球近壁运动以及小球自由下落碰壁弹跳等四种轴对称流动的数值模拟验证结果.

2.1 静止小球绕流验证

根据Wu[24]的试验研究以及Johnson[16]和Tomboulides[25]数值模拟结果,当雷诺数Re=U∞D/νf小于210 时,均匀来流通过一个固定小球所形成的流动是保持轴对称的.这里U∞,D,νf分别表示远场均匀来流速度、静止小球直径以及流体的运动黏度.因此本文分别测试了雷诺数为1,10,50,100,150,200 等六种不同工况,并将本文的数值模拟结果与前人结果进行比较.首先我们给出计算模型如图5所示.模型流向总长为 40D,小球中心距离入口和出口截面距离均为 20D,小球中心距离上侧壁为 20D.其次本文对于边界条件设置:进口采用均匀来流边界条件U∞,出口设置为对流出口边界条件,对称轴设置为对称边界条件,侧壁设置为远场边界条件,浸没小球表面设置为无滑移边界条件.模拟过程中本文采用自适应时间步长,保证最大的库朗数(CFL)不超过0.5.

图5 小球绕流几何参数Fig.5 Configuration for flow past a fixed sphere

这里考虑到本文的工况最高雷诺数为200 以及Pan[26]的数值经验,一个网格无关性验证如表1所示.

表1 雷诺数Re=200 网格无关性验证Table 1 Grid independence test at Re=200

表中 Δx为网格最小尺寸,Cd=为阻力系数(Fz为流向阻力),Lre表示小球尾流回流区长度.通过粗细网格阻力系数和回流区长度结果对比,我们发现粗网格对于本问题是足够解析的.因此后续工况的数值结果均是基于粗网格实现.在图6中,从上到下依次给出雷诺数Re为50,100,150,200 时采用轴对称IBM 计算得到的流线结果,白色半圆表示浸没小球.从图中我们可以看到四种雷诺数下,随着雷诺数的增加回流区长度也是增加的.

图6 不同雷诺数流线对比Fig.6 Streamline at different Re with IBM

图7(a)和图7(b)中,将采用本文算法计算得到的回流区长度Lre以及阻力系数Cd与Johnson[16]和Magnaudet[27]的结果进行定量对比.根据图7(a)可以发现:我们计算的回流区长度在雷诺数为50 和100 时与前人结果均符合较好,在雷诺数为150 和200 时与Johnson[16]结果也符合很好.而根据图7(b)可以发现:我们通过轴对称IBM 算法计算的阻力系数在低雷诺数和高雷诺数与Johnson[16]和Magnaudet[27]的均符合很好.综上所述,本文的轴对称浸没边界算法对于求解小球绕流中轴对称流动是可靠的、有效的.

图7 小球回流区长度和阻力系数对比Fig.7 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of sphere

2.2 静止圆盘绕流验证

除了小球这种特殊几何结构之外,圆盘也是一种典型的轴对称结构.为了更好地说明本文的算法对于固体几何的适应性,我们对于低雷诺数下的轴对称圆盘绕流也进行了相应的测试.根据文献[28-29]的数值研究:对于纵横比χ=D/Hz=10 的圆盘,当雷诺数Re=U∞D/νf小于135 时,远场均匀来流通过固定圆盘后形成稳态轴对称流动.这里U∞,D,Hz,νf分别表示远场均匀来流速度、静止圆盘径向直径、静止圆盘轴向厚度以及流体的运动黏度.本文分别计算了雷诺数为10,20,40,60,80,100 等六种工况.为了更好地与文献[28]研究对比,我们在计算模型的设置上采取与之相一致的参数.计算模型如图8 所示:进口距离圆盘 2.5D,出口距离圆盘 15D,侧壁距离对称轴 6D.

图8 圆盘绕流几何参数Fig.8 Configuration for flow past a circular disk

这里关于边界条件的设置与上一小节的小球绕流中的设置是一样,便不再复述.我们同样采用两种不同分辨率的网格进行计算,用以验证网格无关性,如表2 所示.

表2 雷诺数Re=100 网格无关性验证Table 2 Grid independence test at Re=100

根据表2 中的阻力系数以及回流区长度对比,我们发现粗网格的数值计算结果已经足够求解该问题,所以后续的低雷诺数下的工况我们均采用粗网格进行计算.

图9 中,我们分别给出了雷诺数为10 和100 两种工况下的流场流线图,白色矩形表示浸没圆盘.我们可以直观地看出回流区长度也是随着雷诺数增加而增加的.

图9 不同雷诺数流线对比Fig.9 Comparison of streamline at different Re

除此之外,我们同样对比了圆盘尾迹回流区长度Lre以及阻力系数Cd,对比结果如图10(a)和图10(b)所示.我们发现二者在不同雷诺数时均与文献[28]符合较好,并且阻力系数Cd随着雷诺数的增加是减小的,这与小球绕流结果是相似的.综上所述,本文的轴对称IBM 算法对于圆盘几何的轴对称绕流问题求解结果也是有效和可靠的,同时也说明本文算法对复杂的轴对称几何也是适用的.

图10 圆盘回流区长度和阻力系数对比Fig.10 Comparisons of the length of recirculation zone Lre and drag coefficient Cd of circular disk

2.3 Stokes 流小球近壁运动验证

在前面两个小节中,本文对于静止浸没边界问题已经做了详细的验证和讨论,充分说明了本文的轴对称IBM 算法对于各种复杂几何静边界轴对称流动问题是能够正确求解的.但是工业工程更多的是动边界问题,本小节将通过数值模拟的方法计算匀速小球靠近壁面运动所受阻力,并与经典的Brenner[30]小球碰壁润滑力解析解式(19)和式(20)进行对比

图11 小球匀速靠近壁面几何参数Fig.11 Configuration for a sphere with a constant velocity approaching a wall

图12 给出了不同间隙时的流场压力云图.从图中我们可以看出由于小球不断靠近壁面,间隙流体被挤压向径向外侧流动的同时,间隙处的压力会急剧升高导致小球阻力上升.当间隙特别小时,压力会变得非常大,相应所需要的网格分辨率也就会很高.本文阻力数值模拟结果与Brenner[30]理论结果对比如图13 所示:采用了800,1600 和3200 三种不同分辨率的网格进行计算,很明显三种网格的数值模拟结果在 ε＞0.1 时和理论结果都能符合很好,但是在小球无限接近壁面时,即 ε→0 时结果开始出现偏差,而且网格分辨率越小偏差越大.

图12 压力云图间隙为(a)ε=1.1和(b)ε=0.1Fig.12 Contour of pressure around sphere with (a)ε=1.1 and(b)ε=0.1

图13 数值模拟阻力和理论结果对比Fig.13 Comparison of numerical and theoretical results

这是因为随着间隙越来越小,压力场梯度剧烈变化,想要精确地解析流场,所需要的网格分辨率是非常高的.虽然这种极限小间隙情况很难解析正确,但是在网格分辨率足够的情况下,我们的结果依然是有效的.这可以说明本文的轴对称算法对于给定速度的强迫运动产生的动边界问题也是可以正确求解的.

2.4 小球自由下落碰壁弹跳验证

上一小节中,本文验证了强迫运动的动边界问题,本节我们通过耦合颗粒运动学方程以及流体控制方程来数值模拟小球颗粒在不可压缩黏性流体中自由下落碰撞壁面过程.针对这样的一个球壁碰撞运动,Gondret 等[31]通过实验研究,给出了小球碰撞后的运动轨迹以及相应的恢复系数.文献[32-33]通过数值模拟的手段再现了Gondret 等[31]的实验结果.为了验证本文轴对称IBM 算法的正确性,我们选取与之相同的工况参数设置并进行结果对比.由于本节问题和上一小节类似,上一小节的计算模型在本节仍然继续使用.在进行碰撞数值模拟之前,我们需要引入两个理论模型:一个是润滑力模型,另一个是碰撞模型.首先对于润滑力模型,结合上一小节,我们发现间隙越来越小所需网格分辨率越来越高.直接模拟碰撞将会耗费巨大的资源,而且还不能求解准确,因此在小间隙时一个润滑力修正模型[34-35]的引入是有必要的.再加入润滑力模型后,颗粒运动方程(6)修正为式(21),其中润滑力由式(22)[35]和式(23)[36]给出(式(23)为式(20)的近似形式).其中Un为颗粒垂直壁面法向(z轴方向)速度大小,εal为激活润滑力模型间隙条件.由于本文的模拟参数Re＜210,流场保持为轴对称状态,因此小球速度Up只有垂直壁面法向分量Un

这里因为 1/ε 的存在,导致当间隙趋近于零时校正润滑力趋近于无穷大.这显然不符合物理实际,事实上由于颗粒表面粗糙度的存在,一般当 ε 接近于表面粗糙度时碰撞已经发生(ε=εw).此时润滑力不会继续增大,颗粒进入固体碰撞阶段.这一阶段流体对于颗粒小球的作用力相比于固体的碰撞力可以忽略.为了解析这样的碰撞过程,使用一个经典的软球模型[34],即

其中kn和 ηn分别为刚度系数和阻尼系数且由式(25)给出,en,d表示空气中的干碰撞恢复系数与颗粒材料有关,一般直接参考实验数据.根据Gondret 等[31]文献,本文碰撞算例均采用en,d=0.98 的设置.δs-w表示颗粒与壁面重叠位移.碰撞时间步长 Δtc=[33],Δtf表示流体时间步长.NΔtf为人为设定的碰撞时间,本文中N=10[32-34].实际上碰撞发生的时间是比较短的,这里人为地放大了碰撞接触时间,避免速度改变过于剧烈导致计算误差太大

完整小球与壁面正碰撞数值模拟过程如图14所示:当间隙 ε＞εal时,颗粒只受到来自流体的作用力与重力;当 εw＜ε＜εal时,颗粒除了受到来自流体的作用力及重力之外,还有一项润滑修正力;当ε＜εw时,颗粒和壁面已经发生碰撞接触,颗粒运动由方程(24)控制.

图14 小球与壁面正碰过程Fig.14 Process of a sphere impacting normally on a wall

这里我们对比了冲击壁面Stokes 数为27 和152两种工况结果,具体的物理参数如表3 所示.其中Stokes 数Stc=ρpUinD/(9ρfν),雷诺数Re=UinD/νf,Uin表示小球自由下落后在壁面影响前达到的稳定速度.

表3 物理计算参数Table 3 Physical and computational parameters

图15 (a)和15(b)给出了两种工况下的数值模拟结果与Gondret 等[31]实验结果对比.图中和T*=表示无量纲高度和时间,其中tc为首次碰撞发生时间.我们可以发现小球的轨迹均符合较好.因此本文的轴对称算法对于耦合颗粒运动学方程的动边界问题也是准确的.

图15 颗粒碰撞轨迹Fig.15 Trajectory of the sphere after colliding with a wall

3 结论

本文针对颗粒两相流中的无环流轴对称流动问题,发展了基于2D 笛卡尔网格的轴对称IBM 算法.该IBM 算法通过一个阶梯锐利状界面实现固体浸没边界引入,对于界面上的IB 网格单元通过最小二乘法插值得到.本文通过柱坐标系建立流体控制方程,并采用一个有限体积格式的经典投影算法进行方程离散.同时对方程中由于柱坐标系的使用而产生的源项ur/r2采用隐式格式处理.针对动边界,本文也发展了一个2D 自适应网格技术以提升计算效率.通过小球绕流、圆盘绕流、Stokes 流小球近壁运动以及小球自由下落碰壁弹跳等数值模拟算例,验证该轴对称IBM 算法对于复杂几何边界以及动边界轴对称流动的求解是高效的和准确的.