适合复杂边界及运动边界的CFD方法：浸入式边界方法的研究及应用 | 前沿研究

自然界中的鱼类、鸟类等生物均依靠自身部位 (鱼鳍和鸟的翅膀) 在流体介质中的摆动来获得前进的动力 (Tong et al. 2021, Zhang et al. 2022a), 这一过程通常表现出较高的推进效率以及较低的噪声水平。因此人类始终希望能够深入理解自然界生物流动中的复杂机理, 以此来设计具有更高性能的推进装置, CFD方法是目前广泛采用的一种研究手段。

以鱼鳍为例, 真实的鱼鳍通常具有复杂的外形, 单个鱼鳍在摆动过程中的推力状态可能与鱼的身体或上下游的其他鳍的摆动相互耦合 (Park & Sung 2018, Han et al. 2020, Mendelson &Techet 2021)。与此同时, 由于某些鱼类经常组成庞大的群体, 这种群体行为同样会对单个个体的推力状态以及流场结构产生重要的影响(Weihs 1973, Gazzola et al. 2016, Bao et al. 2017)。这些基本特征使得传统的贴体网格方法在处理这一类问题上面临网格生成困难、多体运动问题难以处理的挑战, 而这些挑战正是IB方法所擅长之处, 因此这一方法也在鱼类游泳等生物流动领域获得了广泛的应用。

为降低预测难度同时把握关键物理因素, 早期对于生物流动的数值模拟研究多采用简化固体模型, 如用振动翼型来表示鱼鳍或鸟翅 (Wang et al. 2014, 2019, Deng & Caulfield 2015, Denget al. 2015, 2016). 而近年来的部分研究工作已经开始将生物体真实几何外形以及摆动规律包含于数值模拟中。图3(a)给出了一种外形类似蝙蝠的鱼类示意图, 对于蝙蝠鱼而言, 其推力产生机制难以通过简化的固体模型进行阐述, 因此Huang等(2020) 通过对真实蝙蝠鱼身体摆动规律的细致观察建立了如图3(b)所示的固体模型, 其中鱼身的摆动通过指定固体边界在不同时刻的空间分布来实现。Zhang等(2022a)采用IB方法对图3(b)中的蝙蝠鱼模型的推力产生机制展开了系统性的参数研究, 其结果表明鱼身沿弦向(或流向)的摆动波长对于推力的产生、推进效率以及能量消耗等参数有着至关重要的影响, 合适的弦向摆动波长有助于获得更高的推进性能, 而摆动频率同样能够对于前缘涡、翼尖涡等涡系结构产生明显的影响进而改变推进状态。上述研究结果对于深入理解包含复杂三维变形鱼类的推力产生机理具有重要意义, 而这一问题对于采用贴体网格的CFD方法而言具有极大的挑战性。

 

正如前文所提到, 鱼类在游泳的过程中多个鳍之间的流动结构会产生相互干涉进而影响推力状态, 贴体网格在处理多鳍干涉问题时往往面临不同程度的网格重构困难, 通常来讲运动边界的数量越多重构过程越复杂, 重构过程消耗的计算资源也会越多。而对于IB方法而言, 由于采用固定的笛卡尔网格, 运动边界的数量并不会对数值模拟的过程产生根本性的影响, 因此在处理此类多运动边界问题上IB方法同样具有明显的优势(Pan et al. 2016)。金枪鱼、鳟鱼等是典型的依靠多鳍干涉产生推力的鱼类, 其中主要包含腹鳍、臀鳍以及尾鳍等之间的相互干涉, Zhang等(2022b) 采用IB方法对金枪鱼中腹鳍与尾鳍之间的相互干涉及推力产生机理进行了数值模拟研究, 结果表明在腹鳍尾迹的激励下, 尾鳍所产生的推力能够出现明显提高, 但腹鳍的推力状态基本不受尾鳍的影响, 其主要作用机理在于尾鳍所产生的前缘脱落涡强度会在于腹鳍尾迹中脱落涡的融合过程中被强化, 这种鱼鳍间的非定常尾迹相互作用机制在鱼鳍间距较大时会表现得更加稳定。鱼类群体游泳现象与单个鱼的多鳍问题是相似的, 本质上也属于多运动边界问题, Peng等(2018)、Kelly和Menzer (2023)均采用IB方法对鱼类群体游泳现象中的非定常涡相互作用机制及其对推力的影响进行了深入研究, 其中均使用多个运动边界来对鱼群进行模拟, 不同个体的尾迹之间相互干涉所形成的流动特征如图5所示。上述研究也充分体现了IB方法在处理运动边界问题上的优势。由于本文主要关注数值计算方法而非物理机理, 因此仅列举部分IB方法应用的案例, 详细的物理机理读者可参考 Zhang 等 (2022a) 所撰写的文章。

 

● 流固耦合

图5 基于IB方法的叶轮机械转静干涉计算。(a)转静叶片排模型,(b)流场瞬时涡量分布(Duetal,2016a,2016b)

 
 

图6 (a)二维柔性细丝示意图,(b)刚性/柔性细丝组合体

发展现代高阶CFD/CAA方法是处理运动边界复杂流动发声问题的主要途径之一 (Wang etal. 2013, Slotnick et al. 2014), 但基于贴体网格的高阶CFD/CAA方法处理涉及流−固−声多物理场耦合的流致发声问题时一般存在网格重构质量和效率相互矛盾的局限性。为了突破这种局限性, 极有希望的解决思路是利用IB方法构造适合高阶CFD/CAA方法的可压缩体积力模型。

事实上, 高阶CFD/CAA方法结合IB方法的思路早在2011年就已经被Seo和Mittal(2011)提出并尝试进行静止物体绕流的直接噪声模拟 (direct noise computation, DNC), 从而实现流场声场一体解。Mittal 等(2008)直接改进了不可压缩流动下的鬼点IB方法,未引入体积力模型, 通过构造高阶插值方法满足壁面边界条件, 并保持了流固界面尖锐的性质。然而,不同于传统低阶CFD方法, 高阶CFD/CAA方法中的鬼点法必须采用足够多的网格点构造高阶插值,使得靠近壁面处笛卡尔网格上的流动计算必须采用高度非对称或者偏侧形式的空间离散格式以及滤波或者人工黏性模板, 这使得靠近壁面处的网格性质同样高度各向异性。当处理不规则几何边界时, 这种各向异性极易激发起数值伪波, 不仅会使流场及声场数值解降阶, 也极有可能产生数值不稳定性现象。这种由于偏侧格式的使用造成声场解降阶的现象也在Kurbatskii和Tam(1997)关于笛卡尔网格下涉及曲线边界时使用偏侧耗散关系保持(dispersion relationpreserving, DRP)格式时所发现, 并且根据Trefethen(1982)的理论探讨, 发现这种现象在频散性有限差分方法中是广泛存在的。Trefethen通过分析二维有限差分方法的群速度揭示了对于波的传播问题, 任何各项异性如非对称模板数值格式、不均匀介质或者计算网格等均可能导致杂波或者也称为寄生波(parasite waves)的激发。而无论是对于锐利界面IB方法(Seo & Mittal 2011)还是Kurbatskii和Tam(1997)针对CAA发展的笛卡尔方法, 虽然其采用的笛卡尔网格能够最大程度上保持各项同性, 但由于数值格式模板的非对称性总会导致不同强度的寄生波, 这对于运动边界发声问题模拟的影响难以评估。据此不难发现对于涉及声学问题的模拟, 尽量避免各项异性对数值准确性和稳定性可能具有显著好处。

从Sirovich (1967) 的流场延拓理论以及Peskin (1972, 1977) 所提出的体积力方法不难看出,Gilmanov等(2003)、Gilmanov和Sotiropoulos (2005)、Mittal等(2008)、Seo和Mittal(2011) 等研究所采用的锐利界面IB方法实际上已经放弃了 Peskin 所提出的流体域延拓以及原始体积力思想, 这直接导致了壁面附近网格上流体计算必须采用非对称数值格式的结果, 这也是CAA计算中所希望避免的。对比之下, Liu和Vasilyev(2007)、Sun等(2012)、Cheng等(2018, 2021b)的工作大致遵循 Peskin 所提出 IB 方法中连续力源的基本思想, 将固体域视为流体延拓后的伪流体域统一考虑, 并分别针对高阶CAA方法发展了Brinkman罚函数IB方法以及基于Su等(2007) 的半隐式IB方法发展了影响矩阵法(influence matrix method, IMM). 由于体积力模型的采用, 在流固界面附近可以采用统一的高阶中心差分格式求解。之后, Komatsu 等 (2016) 分析发现Liu & Vasilyev(2007)的罚函数IB方法不满足伽利略不变性, 从而不适用于运动边界问题模拟, 为此他们修正了能量方程中的源项模型, 使之能够考虑运动边界发声问题, 如振荡圆柱直接噪声模拟。

相比之下, Sun等(2012)、Cheng等(2018, 2021b)发展的IMM方法具有一定优势, 其体积力模型和原始IB方法一样通过一个近似Dirac函数分布到流固界面周围, 因此不需要区分流固界面内外流场, 并且壁面边界条件可以精确施加在流固界面处, 不需要修改边界形状, 也可以处理无黏流动问题。Sun等(2012)利用 IMM 方法对复杂边界的声散射问题进行了数值模拟,声场结果如图8所示, 通过将数值结果与解析解进行对比来对方法的性能进行评估, 结果表明IMM方法能够对复杂边界的声散射问题进行准确模拟并获得与解析解一致的结果。然而, Sun等(2012)、Cheng等 (2018) 的工作主要针对均匀背景流动下的线性声散射问题, 求解线化Euler方程, 其局限性在于忽略了背景流动、非线性以及流体与声学相互作用, 且其采用的奇异值分解方法 (singular value decomposition, SVD) 求解力源方程组耗时巨大, 难以处理复杂运动边界发声问题。除此之外, 其局限性还体现在方法的收敛性上, 和一般贴体网格方法不同, IMM由于涉及边界网格和笛卡尔背景网格, 其求解的体积力收敛性以及稳定性该如何定义仍然模糊不清。这些问题直接影响其是否能够应用到实际复杂运动边界发声的多物理场耦合问题。

图8 多个圆柱散射下的声场压力分布(Sunetal.2012)

采用IB方法进行CAA计算的一个关键问题仍然在于边界处的处理精度,上文中提到,在连续力源法中采用连续分布的函数来代替Dirac函数导致边界处通常只能达到一阶精度, 这与CAA计算中追求的高阶精度是相矛盾的,但是从Sun等(2012)、Cheng等(2018, 2021a)的研究结果来看, 边界处的一阶精度似乎并不会对声学模拟的总体准确性产生特别大的影响, 但是进一步提高边界处的处理精度始终是研究IB方法所追求的目标之一。Liang等(2008)针对间断问题的数值振荡现象, 提出了一种对间断函数的全局描述思想, 将间断函数表示为光滑函数和修正项之和, 其中修正项由间断处的跃度条件确定, 通过这样的转换, 在求解微分方程时, 一个存在间断的问题可以转换为光滑问题, 进一步可将谱方法应用于间断问题而不会带来Gibbs振荡, 其数值结果表明, 采用此种求解方式可有效提高边界处的求解精度, 精度阶数与所构造的阶跃条件密切相关。这本质上属于一种浸入式界面方法 (immersed interface method, IIM) (LeVeque & Li1994, Xu & Wang 2006a, 2006b, Zhong 2007, Gabbard等 2022), 虽然能够有效地提高间断处理精度, 但是阶跃条件的复杂程度以及相应的计算量都会随精度的提升而大幅提高, 如何进一步提高其数值稳定性并拓展其在CAA中的应用仍有待进一步研究。

 

 

图12  基于SA显式壁面模型的湍流边界层速度分布与壁面律对比(Chenetal.2023)

从上述研究结果来看,如何进一步提高壁面模型的准确性是目前IB方法所面临的核心问题之一。显式壁面律的推导过程中忽略了压力梯度的影响,这一假设对于平板流动是适用的,但是对于某些存在明显压力梯度的复杂流动,显式壁面律可能存在较大的误差。为了提升对近壁面流动的刻画精度,Capizzano(2016)采用了一组包含动量、内能以及湍流输运量的边界层流动控制方程,动量方程中包含了压力梯度的影响,在每个时间步对其进行同步离散求解从而给定边界处网格单元数值,测试结果表明在某些压力变化剧烈的区域,包含压力梯度的处理方式能够使结果更接近贴体网格结果。基于边界层方程的壁面函数还被Shi等(2019)、Ma等(2019,2021)用于发展LES算法以实现IB方法在高雷诺数流动问题中的应用,上述学者的LES模拟结果表明,实现准确预测需在边界层附近对LES亚格子应力模型做适当的修正以平衡壁面模型所提供的摩擦应力。Berger和Aftosmis(2018)采用了一组包含压力梯度的常微分方程对边界层流动进行描述,结果相较于显式的壁面律也表现出一定程度的改善。当压力梯度的影响被包含于边界层方程中时,无法推导得到解析的速度分布,因此需要对边界层控制方程进行同步的离散求解,这一求解过程需要耗费一定的计算时间,虽然模型的准确性及适用性有所提升,但是整体计算效率也会出现下降。Constant等(2021)、Xu和Liu(2021)则通过增加边界单元数量的方式在边界单元的最外层实现几乎一致的分布,对比结果表明采用此种处理方式,即使对于方程(6)所给出的显式壁面律,依然能够获得准确且光滑的压力及摩擦系数分布,并且压力系数的收敛性同样出现明显的改善。如何进一步提高IB方法与壁面模型结合时的准确性以及稳定性,同时兼顾计算效率,目前仍有待深入研究。

图13 (a)内流问题示意图,(b)采用笛卡尔网格结合自适应加密所生成的网格系统,(c)流道贴体网格结合自适应加密生成的网格系统(Wangetal。2023)

Wang等(2023)的分析表明,当流道内的运动边界靠近流道边界时,由于识别流道边界所采用的贴体网格仅沿壁面法向具有较高的分辨率,因此网格的展弦比较大,此时由于运动边界与网格单元的相对朝向未知。传统的锐利界面IB方法是面向均匀正交的笛卡尔网格发展而来,通常采用统一的外伸插值距离,而此时由于网格展弦比较大,统一的外伸插值距离会出现失效的情形。为解决这一问题,Wang等(2023)针对锐利界面IB方法提出了一种自适应外伸插值距离的方法,成功将IB方法拓展至一般的曲线网格,使其能够处理边界与不同展弦比切割的情形,结合上文提到的自适应网格加密以及壁面函数,成功对三维亚声速平面叶栅以及NASA转子37内部的跨音流动进行了模拟。结果表明,该方法能够获得与贴体网格一致的结果,并且能够准确预测叶栅叶片表面的压力分布以及跨音转子的特性曲线,其中叶栅叶片表面的压力分布数值模拟结果及其与实验的对比如图14所示。在此基础上,Chen等(2023)采用上述方法对一低压涡轮叶栅以及涵道风扇流动进行了数值模拟研究,结果表明该方法同样能够准确预测不同工况下叶片表面压力及下游尾迹的分布状况,而计算效率相较于采用传统的正交笛卡尔网格而言大幅提高。上述研究表明,采用自适应外伸插值距离能够提高IB方法与不同类型网格的结合能力从而拓宽方法的适用范围,但值得注意的是,Tamaki等(2017)、Berger和Aftosmis(2018)、Chen等(2023)、Wang等(2023)的研究均表明,对于实际工程领域的高雷诺数流动,流场中可能存在流动分离现象,采用湍流壁面模型在这些分离区域可能出现明显的预测偏差,但是目前IB方法在高雷诺数流动的应用中仍然非常依赖于壁面模型来提高计算效率。因此,如何进一步提高该方法对不同类型复杂流动的适应能力是将其推向工程应用所需要解决的关键问题之一,一种可能地解决手段是采用脱体涡方法或大涡模拟等高保真模拟方法来提高复杂流动的预测精度,但其对计算资源的需求使其目前还无法实现规模化工程应用。

图14 三维亚声速平面叶栅叶片表面压力分布数值模拟结果及其与实验的对比(Wangetal。2023)

文章转载自《力学进展》

作者：王 卓，杜 林，成 龙，孙晓峰