一种模拟动边界绕流的锐利界面浸入边界法*

2022-09-17张晋铭张纹惠王文全

爆炸与冲击 2022年8期

郭涛，张晋铭，张纹惠,3，王文全

（1. 昆明理工大学建筑工程学院，云南昆明 650500；2. 四川大学水力学与山区河流开发保护国家重点实验室，四川成都 610065；3. 中国能源建设集团云南省电力设计院有限公司，云南昆明 650051）

绕流是自然界以及实际工程运用中一个十分普遍的物理现象，也是一个经典的流体力学问题。例如：机翼在流动空气中的振动问题；导弹无侧滑大攻角飞行时，背风面分离涡的不对称运动导致的附加偏航力问题；水轮机导叶、叶片在水流作用下的摆动和振动问题；鱼在水中通过鱼鳍的摆尾提供推力；昆虫利用翅膀的拍动使其在空中飞行等等。这些运动及工程都有涉及到翼型、椭圆在流场中的绕流问题。当流体绕过不具有流线型的钝体结构时，在钝体尾部会发生周期性的漩涡脱落，产生周期性的上下拖曳力，也就是升、阻力。随着来流情况的变化，当涡脱频率接近结构的自振频率时，漩涡脱落和结构振动相互锁定，形成共振，也就是涡激振动，其在实际生活中产生的危害十分常见并且影响也十分广泛。例如：1940 年刚通车数月的塔科马大桥发生颤振风毁事件，震惊世界；以及近期的虎门大桥涡激振动事件等。因此，流固耦合问题一直是计算流体力学领域较为关注的一类重要问题。

一般地，研究流固耦合问题时，通常采用贴体网格，这类网格质量良好且生成速度快，但是在实际计算中，为了适应不断变化的边界条件，在每一个时间步中都要重新划分网格，这样就大大增加了计算成本。鉴于这些原因，Peskin提出了浸入边界法，并应用到模拟人类心脏血液流动中的瓣膜运动中。浸入边界法（immersed boundary method，IBM）在整个物理计算区域中采用一套笛卡尔网格，避免了在每一个时间步更新网格的繁琐，使用力场来代替固体边界，使流固耦合问题中的动边界处理起来更加简单和精确。因此，力源项的获得是浸入边界法求解流体问题的关键。Peskin 方法，早期是将固体边界离散为一组由弹簧连接的单元，当边界发生变形时，就会产生一个意图使其回到原位置的回复力，通过边界上的回复力结合delta 函数求得力源项。但该方法涉及到固体的弹性变形，主要适用于柔性体，而且该变形不易通过流场物理量求得，不具有一般性。Goldstein 等在此基础上发展了虚拟边界法来求解刚性边界问题，主要使用反馈力法，配合谱方法来渐进地施加所需的边界条件，该法也称反馈力法。在浸入边界法的发展过程中，Jamaludin 等提出了直接力法，即通过直接构造受力区的速度场来推导力源项，从而确保满足边界条件，基本实现了对边界的准确描述，成为学者们的主要研究方向之一。如Fadlum 等将其应用于流固耦合研究，获得了与贴体网格下计算结果一致性较好的结果；Le 等提出了一种隐式直接力浸入边界法，Wu 等、王文全等提出了一种对速度进行二次修正的隐式直接力浸入边界法；Uhlmann等提出了显式直接力法。直接力法相比于Peskin 的经典IBM 法，无须通过繁琐的delta 函数完成拉格朗日点与欧拉点上信息交换，通过直接施加浸入边界处的速度边界，由此推导出力源项，不依赖于固体的本构关系，也不进行反馈调整，故名直接力法。

按照是否使用delta 函数或能否对界面进行清晰描述，可将IBM 法分为扩散界面（diffused-interface）浸入边界法和锐利界面（sharp-interface）浸入边界法。理论上来讲，经典浸入边界法、反馈力法和直接力法均属于扩散界面浸入边界法。锐利界面法不直接计算力源项来施加边界条件，主要是在界面处通过局部插值，重构边界网格点附近的流动参数（速度），以此作为边界条件将其直接施加到计算域进行求解。理论上可以具有高阶精度，在可压缩流体以及激波冲击问题中均得到了利用。本文将亦采用一种改进的锐利界面浸入边界法，将整个物理区域划分成纯流体区域以及包含固体的次流体区域。在流体次区域边界的法向，构造“虚拟点—受力点—垂足点”的计算结构。受力点速度由流体域内某一虚拟点和边界上垂足点之间的速度线性插值确定,实现复杂动边界的流动数值模拟。本文利用C++编程，通过圆柱绕流经典算例与其他方法的文献结果及试验结果进行对比，验证该方法的有效性和准确性。基于该方法，进一步探究不同轴长比、不同攻角 θ 下的椭圆柱绕流的涡结构分布特征、流线、压力等流场现象和升阻力系数、涡脱频率等水力不稳定现象。

1 锐利界面浸入边界法求解过程

1.1 控制方程

将整个物理区域（包括流体、固体）看作不可压缩牛顿流体的粘性流动。在整个计算域内，流体采用Euler 描述，固体采用Lagrange 描述。其连续性方程和动量方程可表示为

1.2 时间推进

采用二阶时间分布投影法对时间进行离散，来求解控制方程（1）、（2），具体步骤如下。

（1）忽略力源项，计算考虑初始压力作用下的中间速度，即：

1.3 空间离散

2 验证算例

如图2 所示，计算域为一矩形区域，长×宽=15×10，流体次区域为边长1.5的正方形，次区域左侧距流场入口的距离为4.25，次区域中央浸没一个刚性圆柱，其约束方式为全固定边界，圆柱直径为，圆心坐标为(0,0)。计算域左侧为均匀来流入口边界，采用Dirichlet 边界条件，即=,=0；上、下两侧均为无穿透边界；右侧为自由出流边界，采用Neumenn 边界，即∂/∂=0，∂/∂=0 。整个流场采用一套间距为Δ=Δ=0.025的均匀四边形网格，固体边界采取与流体网格尺度相等的等间距离散，时间步长为Δ=0.002。

图1 固体边界处理Fig. 1 Treatment of the solid boundary

图2 整体计算域及边界条件Fig. 2 Computational domain and boundary conditions

图3 为=160 s、雷诺数=300 时圆柱绕流的流场速度分布（范围： - 1≤/≤1 ）。从图中可看出，速度等值线表现光滑，说明计算过程中对速度的更新、修正并没有造成流场的震荡。由此可知，采用设置中间速度并不断更新速度的方法是适用的。同时，固体边界周围的速度分布也符合真实的圆柱绕流流场分布规律。说明利用这种锐利界面法描述固体对流体的耦合作用和运用双线性插值方法计算虚拟点的速度是可行的。

图3 流向速度分布Fig. 3 Isolines of velocity

图4 为一个涡脱落周期内圆柱尾迹涡随时间的演化过程（范围： - 8≤wD/≤8 ，为涡量；实线为正值，虚线为负值）。从中可看出柱体表面的边界层分布、分离剪切层的卷起和圆柱后方的分离区域清晰可见。在一个周期内圆柱尾部出现明显的正、负漩涡交替脱落，向下游发展，呈现出典型的卡门涡街现象，与试验得到的卡门涡街现象一致。说明本文对边界附近的流动参数插值的计算方法及程序是可行的。

图4 一个周期内尾迹涡的演化Fig. 4 Isolines of vorticity against time

图5 为=300 时尾涡脱落稳定后的升力系数及阻力系数的时程曲线，可以看出，随着时间的推进，流场逐渐趋于稳定，在流场稳定后升阻力系数都随时间的发展而呈现出稳定的周期摆动。其阻力系数、Strouhal 数()的全时域统计结果见表1 所示。其中升力系数=2F/ρ，阻力系数=2F/ρ，Strouhal 数=/,为涡的脱落频率。对比表中结果可知，相比于文献[16,29]的数值解，本文的计算结果与文献[30]的实验值更为接近。验证了本方法及程序的准确性。

图5 升力、阻力系数时程曲线Fig. 5 Variations of the lift and drag coefficients with time

表1 本文结果与其他文献结果的对比（ R e=300 ）Table 1 The results comparison of average drag coefficients and Strouhal number at Re=300

3 椭圆柱主动运动算例

3.1 计算模型

目前，以圆柱为对象的钝体绕流的研究已较为充分，对椭圆柱的关注度相对要少一些，但椭圆柱作为介于圆柱和平板之间更具有代表性的钝体，其形状更具流线型，流动阻力低，具有较好的对抗流体的作用，而且剪切边界层的卷起、尾迹涡的脱落、转捩和耗散等也很复杂，其研究对工程实际具有重要的现实意义。取计算域尺寸、边界条件、网格尺寸、时间步长等均与验证算例一致，雷诺数=800。不同之处在于：流体次区域中央浸没的固体为一系列不同轴长比=/的椭圆柱（见图6），其绕枢轴点（原点）旋转摆动时攻角 θ （摆角）的变化规律由下式控制：

图6 椭圆柱绕流计算模型Fig. 6 Computational model of flow around an elliptical cylinder

式中： θ为最大攻角值，取80°； ω =2π为振荡角频率，为运动频率， β 为曲线形状参数。该椭圆柱的摆动形式由水翼在正弦振荡的基础上演变而来（当 β=1 时为正弦波的半个波长）。 β=2.5 代表非正弦振荡，前后0.2 s 为匀速开、关时间段，中间持续0.6 s，其1 个开、关周期内的方波曲线如图7 所示。

图7 一个周期内椭圆柱的瞬时攻角变化曲线Fig. 7 Changes of the instantaneous angle of attack of the elliptical cylinder within a period

3.2 不同轴长比对椭圆柱绕流的影响

由图8 可以看出，随着的不断增大，平均阻力系数以及最大升力系数 ||均逐渐增大。由于椭圆柱垂直于水流方向的投影面积不断增大，当流体流过椭圆柱时，对椭圆柱的冲击力增大，因此椭圆柱对流体的反作用力也会增大，椭圆柱所受到来自流体的阻力必然会增大。但是，相比阻力而言，升力的增大幅度比较小。由图9 可以看出，涡脱频率随着轴长比的变化趋势可以分成3 个阶段。初期：随着轴长比的增大而降低，当=0.2 时，椭圆柱的涡脱频率最大，约为0.6 Hz；第二阶段：当进入某一区域（=0.5～0.7）后，开始稳定下来，约为0.38 Hz；后期：当超出这一区间后，涡脱频率将再次随着的增大而递减，“稳定现象”消失，在=0.9 时，达到最小值。

图8 平均阻力系数、最大升力系数随轴长比的变化Fig. 8 Variation of lift and drag coefficients with axis ratio

图9 涡脱频率随轴长比的变化Fig. 9 Variation of vortex shedding frequency with axis ratio

由图10 可知，当=0.2 时，椭圆柱更具流线型，流动阻力低，尾部并没有明显的分离涡，流场比较平稳；当=0.5 时，尾部流线波动开始逐渐明显；随着轴长比的增大，流态波动更加剧烈；当=0.7 时，在后方开始产生一个尺度较小的分离涡，其位置向上偏离椭圆柱的中心。

图10 不同轴长比下的流线图Fig. 10 The instantaneous streamlines with different axis ratio

图11 为=800 时不同轴长比工况下，同一时刻下的瞬时涡量图（ - 8≤wD/≤8 ，w为涡量；实线为正值，虚线为负值）。可以看出：当=0.2 时，尾迹涡以S 型模态脱落，呈现出上侧负漩涡与下侧正漩涡交替脱落向下游演化发展的卡门涡街现象。脱落的正负漩涡之间的纵向间距、周期较小，形成明显的2 列，尾涡强度（尺寸）也较小，漩涡的分离点基本位于轴上；随着轴长比的增大，卡门涡街规律依然明显，但漩涡强度（尺寸）和周期逐渐增大，漩涡分离点偏离轴，且脱落位置增大，正、负漩涡中心间距减小逐渐趋于一列。

图11 不同轴长比工况下的瞬时涡量Fig. 11 Variation of instantaneous vorticity with axis ratio

3.3 不同攻角对椭圆柱绕流的影响

图12 给出了=0.2，=800 时，不同攻角 θ 下，升、阻力系数的变化规律。

图12 不同攻角对升、阻力系数的影响Fig. 12 Variations of lift and drag coefficients with angle of attack

从图13 中也可看出，当 θ =0时椭圆柱上下表面的压力对称，压差接近零，升力系数很小；随着攻角增大，椭圆柱的上表面的压力，不管从数值和承压面积上来看都明显比下表面大，即上下表面的压差增大，故升力系数随攻角的增大而增大。当攻角发展到临界攻角 θ =45时，压差达到最大，之后升力系数的值随着攻角的增大而逐渐减小。

图13 不同攻角下的瞬时压力场Fig. 13 Instantaneous pressure fields at different angles of attack

由图14 漩涡脱频率曲线可以看到，涡脱频率随着攻角 θ 的不断增大，整体呈现先增加再逐渐递减的趋势。其变化趋势可以分成3 个阶段：初期， θ =2°～5°阶段，随着攻角 θ 的增大而增加， θ =2时，涡脱频率达到最大，约为=0.61 Hz；第二阶段，当攻角进入 θ =2°～10°区域后，有所稳定，递减速率较小，在0.6 Hz 左右浮动；后期，当攻角 θ ＞10后，涡脱频率将随着攻角 θ 增大而逐渐减小，后期曲线逐渐平缓，减小的速率逐渐减慢。

图14 不同攻角对涡脱频率的影响Fig. 14 Variation of vortex shedding frequency with angle of attack

图15 为=800 时不同攻角工况下椭圆柱尾迹涡的对比，由图可见，固体表面分布的边界层、分离剪切层的卷起和椭圆柱后侧的分离区，以及涡的脱落等现象都清晰可见。与圆柱绕流一样，椭圆柱绕流尾迹涡也是呈现出明显的一对正负漩涡从上下两侧交替脱落，向下游发展的卡门涡街现象。但是，椭圆柱摆动的攻角幅值对流态起决定性作用。随着攻角的增加，涡的脱落周期、尺寸、形状和演化形式等均有较大差异，涡脱落的反对称性越发明显。类似于水轮机活动导叶大攻角运行或导弹无侧滑大攻角飞行时，背流（风）面分离涡的不对称运动。攻角越大越容易产生不希望的附加偏航力，或加大水轮机活动导叶转动轴折断的风险。

图15 不同攻角下的瞬时涡量图Fig. 15 Variation of the instantaneous vorticity with the angle of attack

随着椭圆柱摆动攻角的增大，漩涡脱落频率降低，周期加长，漩涡尺寸加大。当攻角 θ ≤30时，在每个周期内，有两个符号相反的涡从椭圆柱上下侧脱落，在下游形成交替出现的涡对，是典型的反对称S 型脱落模态，如图16(a)所示，类似于静止圆柱的涡脱落方式。当 θ =45°～60°时，上下侧脱落的漩涡尾部均拖着一个方向相同的子涡（2 个涡核）；在向下游演化发展的过程中，上侧脱落的负漩涡的子涡与母涡融合，而下侧脱落的正漩涡的子涡则产生分裂，分裂后的次正漩涡位于最上侧，靠近负漩涡，最终尾迹涡表现为3 排，在下游形成了一对涡和一个单涡交替脱落的“P+S”Ⅰ型模态，如图16(b) 所示。当θ≥75时，则情况相反，上侧脱落的负漩涡的子涡发生分裂,分裂后的次负漩涡位于最下侧，靠近正漩涡，而椭圆柱下侧脱落的正漩涡的子涡则与母涡则融合，在下游形成了新的“P+S”Ⅱ型模态，如图16(c)所示。