俯仰振荡翼型流场的离散涡数值模拟

傅翔， 李高华， 王福新

(上海交通大学 航空航天学院， 上海 200240)

0 引言

拍动翼运动是自然界生物最常见的机动形式之一，鱼类，鸟类以及昆虫都是利用拍动翼来获得升力或者推力。一个完整的拍动翼过程可以分解成沉浮振荡和俯仰振荡两种最基本的运动。科学研究表明拍动翼的气动效率远远高于喷气飞机或者螺旋桨的气动效率[1]，所以吸引了越来越多空气动力学家的关注。

国内外已经获得了一些针对拍动翼的研究结果。1973年，Weis-Fogh[2]进行了大量的昆虫实验，提出了一种理论，认为升力的产生主要在于翅膀合拢后的快速打开过程。该发现在研究昆虫飞行的非定常流动中起到了先导性的作用，并且被认为是历史上发现的第一个非定常效应。Weis-Fogh的实验是通过将昆虫拴在风洞中来进行观察和测量的，这与自然界中的实际情况显然是不符合的，不能够真实地反映昆虫的翅膀在拍动过程中的真实受力情况，但是，正是在Weis-Fogh等人的理论提出之后，人们开始从非定常流动角度去研究拍动翼的复杂流动机理。Triantafyllou[3,4]等人对NACA0012翼型进行了水动力试验，翼型做沉浮俯仰组合运动，他们选取了一定范围的参数，对推力系数和效率进行了测量。Anderson[5]利用PIV对二维拍动翼进行了研究，总结出了尾涡形态与振荡频率之间的关系。

而随着数值模拟技术的兴起和发展，越来越多的研究者选择使用数值模拟的方法来研究振荡翼型的运动机理。Guglielmini[6]等人利用数值模拟的方法研究了俯仰沉浮组合运动的振荡翼型周围的流场，认为推力随时间的变化主要受制于翼型的运动细节。Wang[7]等采用基于同位网格的SIMPLE算法对进行沉浮俯仰组合运动的二维振荡翼型周围的粘性流场进行了数值模拟，分析了翼型所受的力和尾涡结构之间的关系，并提出一种用来确定涡的脱落频率的方法。然而传统的数值计算方法在研究拍动翼问题时还存在问题，一方面是因为振荡翼型流场过于复杂，涡的尺度大小不一，导致计算精度难以保证；二是传统数值计算方法在计算非定常问题时所需的计算量太大，常常需要几天甚至几个月的时间才能算完一个拍动周期。

与传统的数值计算方法不同，离散涡方法作为一种拉格朗日求解方法，不需要设置网格，计算量小，在计算拍动翼问题时有着得天独厚的优势。事实上，已经有学者利用离散涡方法对拍动翼问题进行了研究，但由于原有的离散涡方法无法准确模拟前缘涡的脱落，导致计算结果精度较差[8]。本文将利用边界层涡量和分离点处脱落涡量的关系，推导出了高攻角情况下前缘涡脱落位置以及环量的计算公式，再对俯仰振荡翼型流场进行离散涡数值模拟，并结合流场变化对翼型所受升力的演化过程进行分析。

1 离散涡方法的基本理论

二维不可压流动由连续性方程和Navier-Stokes方程控制[9]，即式(1)、(2)。

(1)

(2)

对式(2)取旋度，得式(3)。

(3)

(4)

如果略去粘性扩散效应，则可得式(5)。

(5)

式(5)即等价于无粘欧拉方程。

离散涡方法，首先是离散，把平面上的连续涡量离散化为一系列的点涡。点涡的数量分布和位置坐标事先确定，一种分布方法是根据需要的离散涡数将翼型按照弦长分段，在每一段的四分之一长处设置点涡。

离散后，利用点涡的速度环量值算出离散涡的诱导速度分布。每个点涡的速度由两部分组成,一部分是流场本身的速度，另一部分是由其它离散涡共同诱导出的速度，该速度可以通过毕奥萨伐尔定律算出式(6)。

(6)

最后根据点涡前一时刻的速度和位置，可以计算下一时刻的新位置，采用欧拉近似，有式(7)。

(7)

由此可以模拟出点涡每一时刻的位置，进而模拟出整个涡量场，从而获得流场每一时刻的特征量。

假设离散的涡元总数为N，则在时刻t场点r处涡量的近似值可以表示为式(8)。

(8)

其中δ(r)是狄拉克函数，χj(t)表示t时刻第j个涡元的位置，于是第i个涡元的运动可以描述如式(9)、(10)。

(9)

(10)

式(9)和式(10)即为经典的点涡运动微分方程。

2 离散涡方法的数值化处理

2.1 非定常流场计算

假设一个点涡被放置在(x0,y0)处，那么根据毕奥萨伐尔定律，在xy平面上任意(x,y)处的速度势为式(11)。

(11)

将其分解到x方向和y方向上有式(12)、(13)。

(12)

(13)

用矩阵的形式来表达的话，点涡Γj(xj,yj)在控制点i(xi,yi)处的诱导速度为式(14)。

(14)

式(14)中的r表示这两个点之间的距离，为式(15)。

(15)

将Γj提取出来(亦即代入Γj=1)，可以得到单位点涡的诱导速度公式，将其与控制点i处的法向向量相点乘，可以得到位于j处的单位点涡在控制点i处所产生的诱导速度在当地法向的分量，称其为影响系数aij，有式(16)。

aij=(ui,wi)·ni

(16)

注意到不可穿透边界条件的实现需要满足在翼型表面法向速度为零，对于定常情况，这里的法向速度由两部分构成，一部分是自由来流在当地法向分量，另一部分是所有点涡的诱导速度的法向分量，可以用下面的公式来表达为式(17)、(18)。

(17)

令Vi=-U∞·ni,有

(18)

式(18)即为点涡的边界条件公式，对于控制点1，会受到所有的n个点涡的影响，则可以得到关于控制点1的边界条件公式为式(19)。

a11Γ1+a12Γ2+a13Γ3+…+a1nΓn=V1

(19)

对所有的n个控制点都有这样的方程，将其列在一起，可以得到一个矩阵形式的方程组为式(20)。

(20)

注意到式子中a矩阵由点涡和控制点的位置确定，是已知的，右边的V矩阵由自由来流和各控制点法向量决定，也是已知的，未知的是布置在翼型表面的n个点涡的强度，通过求解这个矩阵方程可以得到初始的定常情况下的翼型表面的点涡强度。

在非定常情况下，会有新生涡的脱落，这时因为要考虑到新生涡与翼型表面涡之间的相互影响，以及新生涡对控制点处产生的诱导速度，所以需要对矩阵方程a进行修正。

当翼型的攻角很小的时候，通常认为没有前缘涡的脱落，仅有尾缘涡脱落，尾缘涡在脱落以后强度不会再发生改变，所以对于每个脱涡的瞬时来说，有N+1个未知数，那么就需要N+1个方程才能求出结果，这里第N+1个方程使用环量守恒来构建，即式(21)。

(21)

同时对矩阵a做出相应的修改，将尾缘涡的诱导速度考虑进去，比较方便的处理方法是在方程右边的V矩阵中添加尾缘涡项为式(22)。

V(t)i=-[U(t)∞+U(t-Δt)w]·ni

(22)

那么对于每一个控制点处的方程也相应的进行修改为式(23)。

ai1Γ1+ai2Γ2+…+ainΓn+aiwΓW=V(t)i

(23)

注意这里的ΓW是新生的尾缘涡的强度，而在环量守恒中的尾缘涡强度Γ尾缘涡则是所有尾缘涡的强度之和。可以得到新的矩阵方程为式(24)。

(24)

其中Γ(t-Δt)是上一时刻流场中的总涡量。

当翼型的攻角较大或者攻角变化较大时，还要考虑前缘涡的脱落，关于前缘涡的脱落还没有统一的处理方式，一种方案是由实验或者经验确定分离点的位置，然后由其他的条件决定新生涡的强度，另一种方案是先通过边界层理论决定初生涡的强度，然后由其他的条件来决定其位置。本文所采用的方法是，将前缘涡布置在翼型前缘附近的弦线上，这样就确定了前缘涡脱落的位置，接下来就是确定前缘涡的强度。如图1所示。

图1 前缘涡脱落示意图

可以将虚线方框内包含的部分当作是当前时刻从分离点脱落的前缘涡，方框内的涡量的增量可以通过公式得出式(25)。

(25)

在该分离点处，有ql≈0，而qu则可以从上一时间点求出，此时有式(26)。

(26)

另外由于前缘上的单位长度的涡量持续地脱落，可以得到式(27)。

(27)

其中ΓLE为前缘上的表面涡的强度，Δl为该涡所在微元段的长度，由此求得前缘涡的强度为式(28)、(29)。

ΓLEW=ΨΓ1

(28)

(29)

将式(28)并入矩阵a，可以得到新的矩阵方程为式(30)。

(30)

至此所有的求解点涡强度的方程已经全部给出。

2.2 气动力计算

以往离散涡方法中都是用非定常伯努利方程来求解气动力，求解过程复杂且难以保证计算精度。这里采用吴镇远提出的涡量矩定理[10]来计算物体所受的气动力，即为式(30)。

∬RsvdR

(31)

其中F为流体施加在浸没在其中并相对于流体运动的固体上的作用力。Rs为固体所占据的区域，R∞为流体和固体所共同占据的无限大的空间，v是速度场，ω是涡量场，ρ是流体的密度，r是位置矢量。

将F正交分解到x,y方向上即可得出这两个方向上的分力。

2.3 计算流程

在本离散涡算法的计算中，均采用了突然启动的模式，即流动从静止状态突然启动而后保持设定好的速度匀速前进。

在求解的过程中，首先设定好并初始化所有参数，然后根据需求在平板或者NACA0012翼型上布置好表面点涡(平板100个，NACA0012翼型300个)。再依据前文给出的诱导速度计算公式，计算出定常情况下的a矩阵，然后进入时间推进，每一个时间点根据涡脱落的设置确定每一个新生的前缘涡和尾缘涡的初始位置，接下来将定常情况下的n阶a矩阵转化为包含前缘涡尾缘涡脱落的非定常情况的n+2阶矩阵。下一步就是根据翼型的运动情况和流体速度计算rhs矩阵，然后解矩阵方程即可获得每一时刻的离散点涡的涡量大小。得到了涡量大小后就可以根据式(12)和式(13)计算出每一个点涡的运动速度，从而确定其下一时刻的位置，这里需要注意的是，翼型表面点涡满足边界条件所以始终保持在翼型上相同的位置，每一步位置会发生变化的是前缘涡和尾缘涡，但是表面点涡的涡量大小会收到脱落的前缘涡和尾缘涡的影响每一时刻都在变化。在每个瞬时获得了每个点涡的涡量大小，涡量变化率和位置，速度这些参数之后，就可以利用涡量矩定理进行力的计算了。重复时间推进后的步骤直到时间达到设定值，程序运行结束。整个离散涡程序计算流程，如图2所示。

图2 离散涡程序计算流程

3 计算结果分析

3.1 算例验证

首先利用离散涡方法计算出了平板翼型从0°攻角匀速抬起至45°攻角的流场。平板翼型弦长为1 m，来流速度为1 m/s，时间步长取为0.05 s。

平板翼型附近的点涡分布及演化过程，图3所示。

可以看到在攻角较小的时候，流场类似定常算例中的流场，有明显的启动涡的卷起。当攻角达到设定值时，前缘涡开始脱落，可以看到前缘涡的脱落对尾缘涡造成了影响，尾缘涡的流线出现了凸起。前缘涡在前缘不断聚集，呈现顺时针旋转，而受到前缘涡影响的尾缘涡也在尾缘聚集，呈逆时针旋转。当前缘涡发展到一定大小时，逐渐被尾缘涡夹断脱落，尾缘涡亦如是，如此往复，形成了如图所示的卡门涡街。

同时程序还用涡量矩定理计算了这一过程中产生的升力，与Darrel K. Robertson[11]的实验结果对比，误差在可以接受的范围内，如图4所示。

图3 翼型附近点涡分布及演化过程

图4 平板翼型抬起过程中的升力

3.2 俯仰振荡翼型流场分析

在验证了程序的可靠性后，对俯仰振荡翼型的非定常流场进行模拟。这里所采用的翼型为弦长1 m的NACA0012翼型，来流速度均取为1 m/s，翼型上分布了300个离散点涡。

翼型采用的振荡运动方式为正弦振荡，这是自然界及工程领域最普遍存在的一种周期性振荡方式，其运动方程可以用下面的式子来表示为式(31)。

αt=α0+αmaxsin(ωt)

(31)

式子中αt表示t时刻翼型的攻角，α0为翼型初始攻角，αmax为翼型攻角的振幅。ω为俯仰频率，这里给出缩减频率k，有式(32)。

(32)

本文主要关心的是俯仰运动的振幅对于升力的影响，因此对所有算例设定相同的参数α0=0，k=0.1，dt=0.02 s。翼型攻角的振幅则取了两个不同的值：αmax=5°和αmax=30°。具体的数值模拟结果和分析如下。

(1) αmax=5°

首先从较小的攻角振幅开始模拟，从程序输出的结果，如图5所示。

图5 αmax=5°时的点涡分布及演化

从图5可以看到攻角振幅为5°时，流场与翼型静止不动时的流场很类似。流动刚开始时，尾缘出现了很明显的启动涡，但是前缘涡没有出现卷起，而是顺着翼型表面向下游平顺的运动。当前缘涡运动到尾缘附近时，可以看到尾缘涡受到了前缘涡强烈的干扰，尾缘涡的脱落不再平顺。当前缘涡进入尾缘涡的范围内时，与尾缘涡交错形成一个小的涡对，由于时间步长的设置，这里的图像中每个小涡对里面的离散涡点数不是很多，看起来不够明显。同时可以看到在初始的时候，翼型正攻角，前缘涡从翼型上表面脱落，而当翼型摆动到负攻角时，前缘涡从翼型下表面脱落。另外可以注意到，即使是从下表面脱落的前缘涡，依然保持着沿着翼型表面向下游流动的趋势，没有出现聚集和卷起。

(2) αmax=30°

攻角振荡幅度增大到30°的流场与αmax=5°的流场有了明显的不同。一个振荡周期内的点涡分布及演化，如图6所示。

图6 αmax=30°时的点涡分布及演化

从模拟出来的结果中可以看到，在攻角较小的时候，流场与αmax=5°的流场比较类似，尾缘涡逆时针卷起，而前缘涡则顺着翼型上表面向下游运动。随着攻角逐渐增大，前缘涡先于尾缘涡开始聚集，顺时针卷起，在翼型表面上方形成了一个巨大的涡团，而后受到前缘涡的影响，尾缘涡不再直接向下游运动，而是在尾缘处逆时针卷起。尾缘涡聚集而成的涡团逐渐发展壮大，然后切断了前缘涡的涡团。而后被切断的前缘涡团又切断了尾缘涡团，形成了一个苹果型的涡对，如此交叉往复，就在流动中生成了一个个较大的涡对形成的涡街结构。可以看到，在翼型攻角逐渐变小的过程中，涡对的大小明显变小了，而且脱落的频率更高。而进入后半个周期后，攻角变为负值，前缘涡开始从翼型下表面脱落，并随着负攻角的变大而出现了逆时针旋转的涡团，相应的，尾缘涡形成了顺时针旋转的涡团，与前半个周期的运动类似，前缘涡和尾缘涡的涡团互相夹断，成对脱落。

3.3 俯仰振荡翼型的升力分析

本文主要关心的是翼型的升力，利用涡量矩定理，对两种不同攻角振幅下的算例进行了力的计算，首先是翼型的总升力系数。在较长周期的模拟中，总升力系数都表现出了周期性，这与翼型周期性振荡的运动本质是相符合的。翼型振荡一个周期的升力系数时程曲线，如图7所示。

图7 翼型振荡一个周期的升力系数时程曲线

从图7中可以看到5°的攻角振动幅度下，升力系数很小，而30°的攻角振动幅度下，升力系数出现了明显的起伏，基本呈现在前半个周期大于零，后半个周期小于零的形式，而且曲线类似于正弦函数的形状，升力系数的极值出现在大约四分之一周期和四分之三周期处，此时翼型的攻角处于最大值，但是角速度为最小值零。30°的攻角振动幅度下，升力系数先达到一个最大值，然后突然变小。这个最大值点出现在t=4s左右，回顾图6，发现此时前缘涡在翼型表面聚集形成了巨大的前缘附着涡团。随后，前缘涡团与翼型表面发生了分离，导致升力突然急剧下降。下降的最低点对应的时间点，流动逐渐从翼型表面分离，流场中尾缘涡卷起，并有切断前缘涡的趋势，翼型上表面出现向前运动的尾缘点涡，与自由流动方向相反，此时翼型上表面存在逆流，之后流动再次附着到翼型表面，升力开始增长。故认为是翼型表面的前缘涡脱落和尾缘涡的逆向流动导致了升力下降，而前缘涡的再次附着又导致了升力重新增加。

4 总结

利用边界层涡量和分离点处涡量的对应关系，推导出了前缘涡脱落位置及环量的计算公式，对离散涡方法进行了修正。用修正后的离散涡方法对俯仰振荡翼型的非定常流场进行了数值模拟，分析了高低两种攻角振幅对应的前缘涡和尾缘涡生长及脱落过程，结合流场对翼型所受升力的演化过程进行了定量化研究，主要结论有：

(1) 攻角振幅越大，升力系数的极值越大，通常出现在攻角最大值附近。

(2) 在一个振荡周期内，对于前后缘涡来说，前半个周期前缘涡的影响较大，后半个周期尾缘涡的影响较大。

(3) 大攻角下流动分离更加剧烈，导致升力系数会出现骤降。

[1] Shyy W, Aono H, Chimakurthi S K, et al. Recent progress in flapping wing aerodynamics and aeroelasticity[J]. Progress in Aerospace Sciences, 2010, 46(7): 284-327.

[2] Weis-Fogh T. Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production[J]. Journal of Experimental Biology, 1973, 59(1): 169-230.

[3] Triantafyllou M S, Triantafyllou G S, Gopalkrishnan R. Wake mechanics for thrust generation in oscillating foils[J]. Physics of Fluids A: Fluid Dynamics, 1991, 3(12): 2835-2837.

[4] Triantafyllou G S, Triantafyllou M S, Grosenbaugh M A. Optimal thrust development in oscillating foils with application to fish propulsion[J]. Journal of Fluids and Structures, 1993, 7(2): 205-224.

[5] Anderson J M. Vortex control for efficient propulsion[D]. Massachusetts Institute of Technology, Cambridge, MA, 1996.

[6] Guglielmini L, Blondeaux P. Numerical experiments on the transient motions of a flapping foil[J]. European Journal of Mechanics-B/Fluids, 2009, 28(1): 136-145.

[7] Zhi-Dong W, Xiao-Qing Z, Yu-Min S, et al. Analysis of the Caudal Vortices Evolvement around Flapping Foil[J]. Journal of Bionics Engineering, 2(04): 195-201.

[8] Lin H, Vezza M, Mcd. Galbraith R A. Discrete Vortex Method for Simulating Unsteady Flow[J]. AIAA Journal, 1997, 35(3): 494-499.

[9] Ginevsky A S, Zhelannikov A I. Discrete Vortex Method[J]. Vortex Wakes of Aircrafts, 2009: 13-32.

[10] Wu J C. Theory for aerodynamic force and moment in viscous flows[J]. AIAA Journal, 1981, 19(4): 432-441.

[11] Robertson D, Reich G, Joo J. Vortex particle aerodynamic modelling of perching manoeuvres with micro air vehicles[C]//51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference 18th AIAA/ASME/AHS Adaptive Structures Conference, 2010: 2825.