黄晓婷, 孙鹏楠, 吕鸿冠, 殷晓瑞, 董嘉徐

(1.中山大学 海洋工程与技术学院, 广东 珠海 519082;2.南方海洋科学与工程广东省实验室(珠海), 广东 珠海 519082)

近年来,微型飞行器或潜水器的发展受到广泛关注,飞行器或潜水器性能的改进离不开对翼型绕流问题的精确模拟。翼型表面流动如分离、附着等现象对翼型的动力性能具有重要影响。为保证在不同工况下飞行器或潜水器运行的可靠性,需要了解翼型绕流的流动特性和机理。这就要求对翼型绕流现象进行深入研究和分析。

翼型的几何形状对改善力学性能具有重要影响[1]。可以发现,目前对于非对称翼型的绕流现象的研究较少。其中,Anyoji等[2]经过研究发现,非对称Ishii弧形翼型由于后缘附近的弧度,其压力阻力要远远小于其他对称翼型,在同等工况下,可获得比对称翼型更高的升力。因此,展开对Ishii非对称翼型更深入研究,准确获取流场信息及其运动特性,可为未来设计出良好动力性能的新型翼型提供新思路。

现有翼型绕流研究中,主流研究方法中主要有试验法及数值模拟法。针对水翼绕流,季斌等[3]从试验和数值模拟角度综述水翼水动力特性。在飞行气翼绕流探究中,Anyoji等[2]采用风洞试验探究及验证Ishii翼型气动性能。然而在试验法中,试验条件复现存在一定挑战,且边界层湍流细节捕捉对仪器的精密性要求苛刻。此外,试验法的成本高昂,且对多工况无法同时进行研究。因此现有研究大多采取数值模拟法展开研究。数值模拟法中，发展较为成熟的欧拉网格法较为常见。但对于网格法而言,翼型后缘狭窄尖角边界层较薄,网格划分边界层存在一定难度[4],在网格离散中需要较高的网格构建技术来满足复杂的边界条件。值得注意的是,近年来发展的无网格法对边界条件、自由液面及大变形等问题[5]处理具有固有优势。

SPH方法作为无网格法中应用广泛的方法,在流体或固体力学问题处理中逐渐发挥重要作用[6]且对于自由液面变形动边界问题的处理上SPH方法具有显著优势[7]。目前,在水动力学领域上得到较多应用[8-9]。Sun等[10]采用δ+-SPH模型对NACA0012翼型固定及自由运动状态进行模拟,分析其水动力特性。Huang等[11]采用SPH方法对翼型绕流的运动性能进行探究;在空气动力学领域,Shadloo等[12]验证了WCSPH可以实现低雷诺数下不同攻角翼型绕流计算;Huang等[13]采用KGF-SPH模型,在低雷诺数下对气翼的升阻力系数进行监测,验证该模型的准确性。可得出,SPH方法在翼型绕流问题研究上具有一定的优势及特色。但由于其本身的精度及稳定性问题,SPH方法在翼型气动性上的研究受到一定限制,翼型绕流中边界层内部流场特征获取仍需进一步的方法改进和技术讨论。

因此,本文克服以往SPH方法难以模拟黏性边界层流动的不足,采用δ+-SPH模型,实现黏性流动下翼型绕流的模拟,展开翼型特性分析。在SPH计算中,高雷诺数下翼型绕流模拟存在三大难点,主要应用以下方法进行解决:①针对流体绕过翼型时,由于负压产生的张力不稳定现象,采用张力不稳定性控制(tensile instability control,TIC)技术,实现Ishii翼型黏性绕流的压力场和速度场等稳定数值模拟。②为了避免粒子的拉格朗日结构引起的计算精度下降现象,在每个时间步初始时刻,施加粒子修正技术(particle shifting technique,PST),确保粒子的均匀分布。③为实现边界层黏性力的精确计算,采用自适应粒子细化(adaptive particle refinement,APR)技术,减少粒子数量,提高计算效率,实现粒子高分辨率同时确保计算精度。

1 δ+-SPH模型

在δ+-SPH模型中,考虑流体的真实物理黏性,在动量方程中施加黏性项,以精确模拟边界层的流动。再有,在连续方程中添加密度耗散项[14],以避免高频的压力振荡。此外,为了避免传统SPH模型中由于粒子的拉格朗日结构引起的计算精度降低现象,δ+-SPH模型针对粒子进行位移修正,引入粒子位移修正技术(particle shifting technique,PST),确保在计算过程粒子分布的均匀性。因此,基于流体弱可压假设下,δ+-SPH模型[15]最终写为

(1)

式中：i和j表示支持域内目标粒子及其对应的邻近粒子；fi是质量力,本文计算均取fi=0；ρi,mi,Vi,ui,ri分别代表粒子i的密度、质量、体积、速度和位移；光滑核函数Wij采用Wendland的C2核函数[16],光滑长度h=2Δx,其中Δx是初始粒子间距；耗散项Di和πi的具体参数可参考文献[15]。Fij为压力项,在本文中采用张力不稳定控制技术进行处理。

为确保粒子的均匀分布,在δ+-SPH模型中,动量方程计算得出的粒子位移ri需要施加人工位移修正量δri,δri具体定义为

(2)

2 张力不稳定性控制技术

在传统SPH模型中,由于负压会导致张力不稳定现象[18],流场信息将无法很好呈现。基于Sun等[15]的研究结果,为更好实现漩涡附着、分离、脱落细节的捕捉,展开对Ishii翼型流场规律分析,本文采用TIC技术改进压力梯度项。对压力项Fij进行改进[15]，如(3)式所示

(3)

DF表示液面区域自由液面粒子及相邻粒子的集合。

3 拉格朗日拟序结构可视化

在欧拉算法的模拟结果中,流场特征量(如涡量)通常基于欧拉框架下的瞬态速度场计算得出,难以给出流体介质在某个时间段内的输运特征,而本文采用拉格朗日拟序结构对流场中漩涡结构进行可视化,可更加清晰地获取流场细节。拉格朗日拟序结构通过计算有限时间李亚普诺夫指数(finite-time Lyapunov exponent,FTLE)实现。在FTLE云图中呈现的脊线为拉格朗日拟序结构[19]。

在SPH计算程序中采用逆向FTLE展开对FTLE计算,以便能清晰反映漩涡结构边界特征[20]。逆向FTLE计算方程为

(4)

4 多级粒子分辨率技术

翼型的首尾部往往较为狭窄,翼型绕流中表面边界层内湍流流动的捕捉对于网格法和无网格法均具有一定挑战性。因此,本文将利用多级粒子分辨率技术[21],将翼型表面离散成分辨率高的固定虚粒子,实现精确捕捉湍流流动细节,提高边界黏性力的计算精度。多级粒子分辨率技术原理如图1所示。

图1 多级粒子分辨率技术原理示意图[21]

主要原理阐述如下:不同分辨率的粒子定义为不同的粒子集合。首先定义翼型附近区域为细化区,在细化区边缘定义过渡区,如图1b)所示,充当其边界条件,获取进入细化区域母粒子的信息。过渡区域的宽度设置要大于核函数的半径,防止其边缘处核函数截断而造成计算误差。在计算中,参与流场SPH计算的粒子称为SPH粒子;活性关闭,不参与SPH计算的粒子为非活性粒子。在细化区,母粒子分裂,如图1a)所示,在正方形4个顶点处形成更小粒子即子粒子,子粒子活性激活,参与流场的SPH计算,即为SPH粒子。而母粒子在撕裂之后,作为非活性粒子不再参与SPH计算,但仍然随着流场运动,通过从周围细化的子粒子插值获得流场信息如速度、压力等。相反的,在过渡区域,母粒子为SPH粒子,子粒子属性从周围SPH粒子即母粒子插值得出。采用Shepard插值方法,非活性粒子的物理量φ从周围SPH粒子完成插值,插值公式为

(5)

经过此过程,2层分辨率粒子完成了相互之间的信息交换。多级粒子分辨率技术对局部流场粒子加密,在不同的分辨率区域之间,过渡区域及细化区域并不随着粒子的移动而移动。在计算中,当母粒子撕裂而来的子粒子流出细化区及过渡区域时,不发生粒子合并而将直接被删除,而母粒子流出细化区域时,将正常参与SPH流场计算。基于以上粒子细化思想,在多层分辨率粒子间,保持粒子尺度一样,每层粒子的过渡区域流场信息从下一层插值取得,从而实现不同层次间粒子信息交换,实现多级分辨率的SPH数值计算,提高流场计算精度。在下文将首先通过圆柱绕流基准算例验证在δ+-SPH中添加多级粒子分辨率技术的有效性,从而应用于翼型绕流问题的研究中。

5 基准算例验证

圆柱绕流问题是流体力学中的经典问题,是检验数值算法的经典算例。选取雷诺数Re=3 000,人工声速c0=15 m/s时的圆柱绕流模拟算例,与试验结果对比以检验δ+-SPH模型的数值精度。

在该算例中,施加多级粒子分辨率技术,采用5层分辨率不同的粒子。圆柱表面的粒子间距为D/Δx=400,最远处粒子间距为D/Δx=25。在计算中,粒子总数为1.45×106。假设在计算域中粒子分辨率均取为D/Δx=400,粒子总数将达到7.68×107。可见,多级粒子分辨率技术对重点区域流场进行加密,在保证计算精度条件下,有效减少了粒子总数,从而减少计算时间,提高计算效率,对实现SPH方法在工程上应用具有重要意义。

图2 雷诺数Re=9 500圆柱绕流流场特征

图2中选取圆柱上表面的δ+-SPH计算图展示圆柱绕流的漩涡发展过程。通过漩涡流线图与试验结果[22]对比,可以发现试验结果与δ+-SPH结果吻合良好。δ+-SPH模型与多级粒子分辨率技术相结合,漩涡发展细节能较好被模拟,δ+-SPH结合多级粒子分辨率技术可以较为准确模拟绕流现象中分离、附着现象。SPH中拉格朗日拟序结构能清晰展示漩涡演化过程,呈现清晰的内部细节,与试验结果基本保持一致,说明拉格朗日拟序结构获取流场特征的高精度。基于以上讨论,可见SPH方法在漩涡结构捕捉等方面具有独特优势。在下节中,本文将采用多级粒子分辨率技术对计算域粒子进行处理,结合δ+-SPH计算结果中的拉格朗日拟序结构展开对翼型绕流现象研究。

6 Ishii翼型绕流

6.1 Ishii翼型绕流的初始工况和边界条件设置

Ishii翼型具有翼厚比低、前缘半径小等特征,翼型形状[23]如图3所示。本文模拟将Ishii翼型固定在矩形流场中,流场左侧为流入边界,来流速度设为U0=1 m/s,右侧为流出边界,上下边界为自由滑移壁面边界。翼型的上表面施加无滑移边界条件。

图3 攻角α=6°时的Ishii翼型形状

对Ishii翼型表面的无量纲量升阻力系数测量计算,获取翼型的气动力特性。对阻力系数CD,升力系数CL定义为：

(6)

式中，fD和fL分别定义为绕流Ishii翼型的阻力和升力,在SPH程序计算中,这两者的计算基于流体粒子与结构粒子之间相互作用力的平衡算得。在本文翼型模拟中,弦长L=1 m。

6.2 多级粒子分辨率技术对计算效率的提升

为提高计算效率,采用6层不同的粒子分辨率,在翼型表面的最高粒子分辨率为L/Δx=800,最低粒子分辨率即距翼型最远处取L/Δx=25。在本计算中粒子总数为6.40×105,如果采用L/Δx=800对全流场进行计算,粒子总数将达到3.07×107,可见采用多级粒子分辨率技术,粒子总数减少将近43倍,除去不同尺度粒子之间的插值耗时,多级粒子分辨率技术可将计算效率提高约20倍。

6.3 雷诺数Re=3 000下绕流模拟结果和验证

首先,在雷诺数Re=3 000、攻角α=0°条件下,模拟了Ishii翼型绕流特性。阻力系数的计算结果如图4所示。FVM结果与SPH结果基本吻合,验证了SPH方法的计算精度。

图4 雷诺数Re=3 000, 攻角α=0°Ishii翼型阻力曲线

6.4 雷诺数Re=30 000下绕流模拟结果和分析

6.4.1 升阻力系数

随后,在雷诺数Re=30 000、攻角α=0°,3°,6°条件下,模拟了Ishii翼型绕流的特性。α=6°时,翼型阻力及升力系数的SPH计算结果与FVM结果对比如图5所示。在粒子初始阶段到历经剧烈抖动之后,最终CL在0.60～0.85之间动态波动。同样地,曲线可发现阻力CD系数在0.044上下波动。

图5 雷诺数Re=30 000下的升阻力系数时历曲线

表1给出了本文计算结果及文献中试验结果[23]与其他数值结果[24](2-D Lam、2-D RANS)的对比。可以看出,多种算法的升阻力系数值基本吻合,SPH结果与2-D Lam计算结果吻合最佳。通过比较图6中攻角α=0°,3°,6°的升力系数可以看出,攻角变大、升力系数增加。此外,攻角增大引起黏性漩涡的非周期性脱落,使得升力系数随时间产生振荡。

表1 升阻力结果验证

图6 不同攻角α下升力系数时历曲线

6.4.2 SPH计算黏性绕流问题的速度场与压力场

图7展示了雷诺数Re=30 000、攻角α=6°条件下,Ishii翼型绕流的速度场和粒子分布情况。

图7 Ishii翼型绕流的速度场和粒子分布

对比t=0.7 s时刻,施加与不施加PST时的粒子分布,可见：不施加PST时,翼型周围粒子分布十分不规则,且速度场呈现出不稳定;而施加PST,翼型表面边界层内流体粒子分布十分均匀,保证了速度场计算的精度。

在SPH模拟中,实现稳定的压力场计算至关重要。然而在传统SPH算法中,常常存在压力振荡及张力不稳定性现象。在δ+-SPH模型中,通过在连续方程中添加密度耗散项[14],可得到光滑的压力场。针对张力不稳定现象,本文采用张力不稳定性控制技术进行解决。为验证其有效性,图8中给出Ishii翼型绕流模拟中施加TIC与不施加TIC的压力云图。可见,不施加TIC时,翼型周围流场中出现数值空洞,导致计算结果失真;施加TIC时,Ishii翼型周围压力分布均匀,且没有数值噪声。

图8 压力场云图(t=4.076 s,Re=30 000)

以上对比说明,δ+-SPH模型结合TIC技术可以有效消除压力噪声,进而较好实现翼型的压力数值变化预报。此外,观察压力云图可以发现,翼型首部附近区域存在巨大负压区,这是引起张力不稳定现象的主要原因。因此,采取张力不稳定性控制技术是准确获取流场特征的重要前提。

6.4.3α=6°时Ishii翼型绕流中漩涡结构分析

图9 翼型周围流场中的漩涡结构(t=12.326 s,Re=30 000)

图9给出了翼型周围流场中的拉格朗日拟序结构和涡量场。可见,翼型表面附近高分辨率粒子的布置使得边界层内湍流细节能能够清晰地捕捉,反映出多级粒子分辨率技术在高雷诺数下黏性绕流问题中的重要作用。相比于涡量场(网格法较常使用该方法捕捉漩涡结构),在SPH框架中捕捉的拉格朗日拟序结构更能清晰重现涡结构的轮廓和内部细节。此外,相较于涡量场云图,FTLE云图受阈值的影响不敏感,更加便于进行流场细节分析。这也是SPH粒子法处理翼型绕流问题的优势之一。

从图9可见,在翼型形状及压力梯度影响下,Ishii翼型漩涡的生成、附着、分离等主要发生在翼型上缘。基于多级粒子分辨率技术,可以发现,尺度很小的漩涡也能清晰捕捉,再次说明了多级粒子分辨率技术的有效性。

7 拓展应用:拍动翼型绕流问题

为进一步展示SPH在处理运动翼型绕流问题中的优势,本节采用δ+-SPH模型对拍动水翼问题进行模拟和验证。在网格法中,对运动翼型绕流问题的模拟容易引起网格变形而降低计算精度。即使采用重叠网格法,也常常因涉及不同网格间的插值而使算法变得复杂。相比而言,得益于SPH的拉格朗日粒子特性,对于大幅运动边界问题的模拟较为方便。

本节选择头部半圆尾部尖角的拍动翼型(具体细节参考文献[25])进行SPH计算与验证;采用无量纲参数包括雷诺数Re、斯特劳哈尔数Sr、幅值无量纲系数AD表征流体黏性、翼型拍动频率、翼型拍动振幅对翼型绕流模拟的影响。拍动幅值无量纲系数AD=2A/D,其中,A为翼型拍动最大幅值,D取为翼型头部半圆直径。

图10给出了雷诺数Re=440,斯特劳哈尔数Sr=0.035,无量纲幅值系数AD=1.47条件下,SPH模拟的漩涡结构与试验结果[25]的对比。在t1时刻,SPH结果与试验结果[25]吻合较好,体现了δ+-SPH模型在运动翼型绕流模拟上的精度和独特优势。观察SPH不同时刻的结果,可见拉格朗日拟序结构能呈现拍动翼周期拍动时漩涡演化过程,这为未来翼型推进性能与漩涡关系的研究提供了新思路。

图10 拍动水翼绕流尾迹的试验结果[25]与SPH结果

8 结 论

本文采用δ+-SPH方法对不同雷诺数和攻角下的翼型绕流问题进行了数值模拟研究,分析了升阻力特性和流场漩涡特征,并将SPH方法从固定翼绕流拓展应用到了拍动翼绕流的模拟中。研究结果表明:

1) SPH方法在模拟真实黏性边界层和漩涡运动的流体动力学问题中具有竞争力。不同于网格法,SPH方法易于实现拉格朗日拟序结构的捕捉,动态呈现漩涡结构演化过程,因此在未来探究漩涡与流体动力特性关系上具有重要优势。

2) 采用TIC技术及PST技术能够有效提高翼型绕流的速度场和压力场计算精度,并实现对升阻力系数较高精度的预报;采用多级粒子分辨率技术能够有效减少粒子总数,缩短计算时间,从而便于在工程中应用。

3) SPH方法不需要生成网格,能自动追踪运动和变形的流固界面,而且该方法的拉格朗日和无网格粒子特性使其在模拟运动翼型绕流问题上具有先天优势和巨大潜力。在未来研究中,利用该优势,SPH方法在翼型绕流的涡振或颤振问题研究中有望发挥重要作用。