涉及移动边界的流动问题一直以来都是计算流体力学(CFD)的研究热点．目前，嵌套网格方法已广泛应用于模拟涉及移动边界的流动问题[1-3]．嵌套网格方法的关键在于如何处理重叠区域，跨越不同网格间的内部边界传递流场信息以耦合分离的背景域和子域．一种比较直接的方式是删除重叠区域的部分网格，然后在该区域生成网格连接分离的背景域和子域[4]．这种方式消除了重叠区域，使原网格变为贴体网格，但网格的生成十分繁琐，而且会降低重叠区域的网格质量．另一种方式是根据流场信息在内部边界上施加约束，然后进行迭代求解，耦合不同的计算域[5]．更加高效可行的方式是通过插值过程传递流场信息．近年来，关于重叠网格的插值方法研究成果较多，文献[6]借助于supermesh方法，通过二阶泰勒插值过程传递信息；文献[7]提出了一种四阶精度的拉格朗日插值格式；文献[8]借助于浸入边界法在重叠区域人为设置了一层无滑固壁边界，修正了质量不守恒问题．关于重叠网格方法的研究进展可以参考文献[9]，不再详细叙述．尽管不同的插值格式之间差异较大，但根据信息能否同步传递，可以分为显式插值格式和隐式插值格式．其中，显式插值格式即逐子域迭代法，根据相邻网格的流场信息，建立狄利克雷、纽曼及罗宾传输条件，作为边界条件施加在计算域上，然后在每个计算域内独立求解控制方程[10]．流场信息不能同步传递，同样须构建迭代过程求解，但计算精度和收敛性较差．隐式插值格式能够跨越内部边界同步传递流场信息，改善计算精度和收敛性．文献[11]针对压力修正求解格式，将压力修正量作为插值变量分配到插值点上求解压力修正方程，数值算例的结果表明这种做法是有效的，但难以应用于其他离散格式．另一种方式是将传输条件嵌入到代数方程组中，在代数层面进行耦合[9]，但需要传递方程残差来施加纽曼条件，只能够采用迭代求解器求解．本研究提出了一种简单、普适的隐式插值格式．首先基于速度和压力保持连续性的特征建立了插值条件，然后将插值条件嵌入到代数方程组中，通过数学运算重构代数方程组．虽然该方法是通过数学运算建立的，但具有明确的物理意义．插值条件可视为狄利克雷条件，系数矩阵和右端项的重构则施加了纽曼条件．通过点对点元素替换的方式可以简单高效地重构代数方程组，对于原始CFD代码改动较小，易于实现．1 隐式插值方法1.1　基本思想插值条件可表示为X=CX'，(1)式中：X为插值点处的插值变量；C为插值系数矩阵；X'为贡献点处的插值变量．根据不同的离散格式离散控制方程，可建立代数方程组Ma=b，(2)式中：M为系数矩阵；a为未知量，一般为速度或压力；b为右端项．建立插值格式的关键就在于如何耦合求解式(1)和(2)．考虑到插值点的插值变量可由贡献点处变量及插值系数矩阵确定，因此只须求解与贡献点相关的方程即可．若选择压力和速度为插值变量，则可将式(1)直接嵌入到下式MCX'=b．(3)上述方程组不涉及插值点处未知量的求解，但系数矩阵为MC，在方程两端同乘矩阵CT，CTMCX'=CTb．(4)重构后的方程组只涉及贡献点处未知量的求解，而且通过嵌入插值条件实现了不同计算域的耦合．重构方程组可表示为M'X'=b'，(5)式中：M'=CTMC为重构系数矩阵，仍然保持正定对阵的性质；b'=CTb为重构右端项．通过数学运算将插值条件嵌入到代数方程组完成了代数方程组的重构，实现了不同计算域的耦合．由于插值条件被直接嵌入到代数方程组内，因此相邻计算域能够同步传递流场信息，相较于基于贴体网格的求解程序，只须替换系数矩阵和右端项的部分元素，对原始代码改动较小，而且没有额外增加工作量．从代数层面出发，能够适用于任意离散格式，不强制采用迭代求解器求解．1.2　挖洞挖洞是嵌套网格的前处理过程，主要目的是筛除部分单元及建立插值条件．在求解过程中嵌套网格方法一般对背景域施加边界条件，子域的所有流场信息都是未知的，所以子域外边界必须通过插值过程获得流场信息．采用基于距离的线性插值方法，插值点所在的背景单元即为该点的贡献单元．如图1所示，子域外边界的贡献单元已用黄色标注．贡献单元的所有节点均为该插值点的贡献点，相应的插值系数可根据相对距离求解，但须保证插值系数的总和为1．10.13245/j.hust.220108.F001图1挖洞过程在确定子域外边界的贡献单元之后，可设置重叠区的厚度，然后向内搜索重叠区单元，图1中设置了两层重叠区单元．较大的重叠区能够确保流场信息跨越内部边界时光滑过渡，略微提高计算精度，但单元的增加会降低计算效率，在计算中可根据具体问题进行设置[12-13]．重叠区的内边界定义为洞边界，洞边界所包围的单元即为洞单元，洞边界是不许参与计算的．但为反映实体边界的存在，洞边界须通过插值过程由子域提供流场信息，贡献单元及插值系数的选择与子域外边界节点类似．1.3　耦合方法耦合方法的求解过程：构建代数方程组MX=b；通过挖洞建立插值条件X=CX'；重构系数矩阵M'=CTMC；重构右端项b=CTb；求解方程组M'X'=b'；求解插值点变量X=CX'．求解涉及移动边界的流动问题时，上述步骤在每个时间步内重复执行．插值条件可视为狄利克雷边界条件，对于某插值节点i处对应的插值系数满足∑j=1ncij=1，(6)式中：cij为插值点i关于贡献点j的插值系数；n为贡献点总数．压力和速度跨越内部边界时能够保持连续性，贡献域对插值域施加了狄利克雷条件．插值域反过来对贡献域施加了纽曼条件．一般而言，纽曼条件是通过传递方程组残差施加的，因此必须采用迭代求解器求解．本文方法通过重构系数矩阵和右端项对贡献域施加了纽曼条件，所以可以采用任意求解器求解．本文方法通过数学运算建立的耦合格式，然后给出明确的物理意义．通过点对点元素替换的方式重构代数方程组．根据节点属性，可以将所有节点分为插值点、贡献点和普通点三类．重构右端项过程中，当节点i为普通点时，右端项元素保持不变，仍为bi；当节点i为贡献点时，右端项元素由该节点处ni个插值点插值得到，即∑i=1nicijbi．重构系数矩阵的过程类似，当节点i和节点j均为普通点时，保持不变；当节点i为贡献点时，系数矩阵元素为∑k=1ncikmik，mik为M中的元素；当节点i和节点j均为贡献点时，稀疏矩阵元素为∑k=1ni∑l=1njcikcjlmkl．采用点对点元素替换的方式重构系数矩阵和右端项，只须根据节点属性将相应的元素乘上插值系数替换即可，对于原始代码改动较小，易于编程实现，而且并没有增加额外的工作量，计算高效．1.4　动态耦合格式在CFD的数值求解中，通常根据离散格式建立增量形式的代数方程组求解控制方程．对于静态嵌套网格方法，上述格式能够直接应用，因为插值条件恒定不变，速度和压力的增量同样满足插值条件．求解涉及移动边界的流动问题时，插值条件会随子域的移动而改变，速度和压力的增量不满足插值条件，须重新推导耦合格式．增量方程为M(Xt+1-Xt)=Δbt+1，(7)式中：Xt+1为t+1时刻物理量；Xt为t时刻物理量；Δbt+1为增量项，以Xt+1为未知量，可表示为MXt+1=Δbt+1-MXt，(8)考虑到全量仍然满足插值条件，即Xt+1=Ct+1Xt+1'，(9)将式(9)嵌入到式(8)，Ct+1TMCt+1Xt+1'=CTt+1(Δbt+1-MXt)，(10)以增量为未知量，可改写为Ct+1TMCt+1ΔXt+1'=CTt+1[Δbt+1-M(Xt-Ct+1Xt')]. (11)与式(4)相比，式(11)的未知量为速度和压力增量，能直接采用增量格式求解．重构系数矩阵与式(4)一致，仅右端项发生变化，通过点对点元素替换的方式，右端项仍然能简单高效地进行重构．2 算例验证2.1　圆柱横向振荡直径为D=1的圆柱位于[-15,15]×[-10,10]的计算域坐标原点处，圆柱运动速度为u(t)=-u∞cos(2πt/T);                v(t)=0,式中：u(t)和v(t)为圆柱横纵向运动速度；u∞为最大运动速度；T为运动周期．计算参数取值与文献[14-15]一致，雷诺数取为Re=100，最大运动速度为1 m/s，运动周期为5 s．计算所采用的网格模型如图2所示，圆柱附近的背景网格单元尺寸为0.08，子域网格单元尺寸为0.03．所有边界均为滑移固壁．10.13245/j.hust.220108.F002图2圆柱横向振荡计算网格模型为定量地分析计算结果，验证本文方法的精确性，选取x=-0.6，0.6两个截面对比了流场速度．图3和4给出了t=nT+T/2时刻x=-0.6 m截面的水平速度和竖向速度，并与其他文献结果及实验结果并行对比．图5给出了阻力系数Cx=Fx/(0.5ρu2∞D)对比，其中：Fx为阻力；ρ为流体密度．由图3~5可知：当前结果与文献结果符合较好，说明本文方法模拟涉及已知运动的移动边界流动问题具有较高的精度．10.13245/j.hust.220108.F003图3t = nT+T/2时刻x=-0.6 m截面水平速度对比10.13245/j.hust.220108.F004图4t = nT+T/2时刻x=-0.6 m截面竖向速度对比10.13245/j.hust.220108.F005图5阻力系数对比2.2　涡激振动圆柱具有水平向和竖直向两个方向的自由度，运动方程如下X¨+2ξ(2π/U*)X˙+(2π/U*)2X=2CD/(πm*);Y¨+2ξ(2π/U*)Y˙+(2π/U*)2Y=2CL/(πm*),式中：ξ为阻尼比；U*为折合流速；m*为质量比；CL和CD分别为升力系数和阻力系数；X=x/D和Y=y/D为圆柱无量纲位移，D为圆柱直径．计算域大小为[-15D，45D]×[-20D，20D]，初始时刻圆心位于坐标原点处．计算参数取值与文献[16]相同，雷诺数Re=200，ξ=0.01，U*=5.0，m*=4/π．网格模型如图6所示，圆柱附近的背景网格单元尺寸为D/50，子域单元尺寸为D/125．左侧为速度进口边界，右侧为自由出流边界，上下两侧为滑移边界．10.13245/j.hust.220108.F006图6涡激振动计算网格模型计算开始时先将圆柱固定，下游产生稳定的卡门涡街后释放圆柱，可在圆柱尾部观察到如图7所示的双涡，圆柱振动频率为0.187，与旋涡脱落频率一致．由于流畅拖拽力的原因，圆柱振动中心并不是坐标原点，如图8所示，当前方法计算得到的圆柱振动中心为(0.63D，0)，与文献中的(0.62D，0)接近．为方便对比，将圆柱振动中心平移至坐标原点，在此基础上绘制了“8”字形振动轨迹，如图9所示，模拟结果与文献[16-18]结果基本一致，验证了本文方法的准确性．10.13245/j.hust.220108.F007图7圆柱周围涡量云图10.13245/j.hust.220108.F008图8圆柱运动轨迹10.13245/j.hust.220108.F009图9圆柱运动轨迹对比2.3　颗粒沉降计算参数取值与文献[19]保持一致．直径为2.5 mm的圆形颗粒在长0.02 m、高0.06 m的区域内进行自由沉降，令计算域左下角的坐标为(-0.01，0)m．初始时刻，颗粒圆心位于点(0，0.04)m处，整个计算域内流体以及圆形颗粒运动速度为0 m/s．流体密度取为1 000 kg/m3，动力黏性为0.01 kg/(m⋅s)，颗粒密度为1 250 kg/m3．网格模型如图10所示，其中背景网格的正方形单元尺寸为D/8，颗粒附近的子域网格单元尺寸为D/32．整个计算域的外边界均为无滑固壁．10.13245/j.hust.220108.F010图10颗粒沉降计算网格模型图11给出了不同时刻的流速分布云图．计算开始时，颗粒在重力作用下自由沉降，颗粒运动速度不断增大，并在颗粒两侧形成两个反向旋转的漩涡．随着时间的推进，颗粒运动对整个计算域流场的扰动范围不断扩大．t=0.78 s时颗粒已经沉降到计算区域底部．10.13245/j.hust.220108.F011图11不同时刻的流速分布云图(色标单位：m·s-2)为了定量、直观地分析计算结果，图12绘制了颗粒y向速度及y向坐标随时间变化曲线，由图12可知：圆柱在开始沉降后很快达到最大沉降速度，然后匀速沉降，在t=0.73 s之后在底面效应的作用下速度迅速降为0 m/s．此外，将计算结果与文献[16，19-21]结果进行比较，符合度较高．10.13245/j.hust.220108.F012图12颗粒y向速度及y向坐标随时间变化曲线3 结语提出了一种用于动态嵌套网格的隐式插值方法．基于插值条件，通过数学运算将代数方程组重构以耦合分离的计算域．重构后的代数方程组与原方程组的格式相似，系数矩阵仍保持对称正定的性质．本方法是基于狄利克雷和纽曼型传输条件建立的，其中，基于变量连续性建立的插值条件可视为狄利克雷条件，而系统矩阵和右端项的重构则成功施加了纽曼条件．此外，针对CFD模拟中通常采用的增量求解格式重新推导了右端项的重构，使该方法能够直接模拟设计移动边界的流动问题．采用点对点元素替换的方式法能够简单高效地重构代数方程组，与基于贴体网格的CFD方法相比，只须对原代码进行少量修改，并不增加额外的计算量．最后，通过三个基准数值算例验证了方法的可靠性和准确性．该方法能较准确地模拟具有移动边界的流动．