2. 空军工程大学 航空航天工程学院, 陕西 西安 710043
流动的稳定性是翼型绕流的壁面流动分离,近尾迹的漩涡脱落,涡系的演化与相互作用等复杂物理现象的发生根源。由于不同的流动稳定性决定了流场的主要流动特性和翼型的受力情况,因此对翼型绕流的流动稳定性深入理解,对飞行器设计,旋转机械的稳定运行等工程应用具有重要的指导意义。
对翼型绕流流动稳定性的理论研究主要采用线性稳定性理论进行分析。其思路是忽略扰动引起的非线性项,得到关于扰动的线性控制方程,进而研究扰动在线性方程中的发展规律。这种方法忽略了非线性的耦合作用,模态和模态之间相互独立。对于复杂的非线性流动现象(如非线性饱和、转捩等)无法给出合理解释。针对这一问题,Stuart[1]在Landau研究成果的基础上从线性稳定性理论提出了著名的弱非线性方程Stuart-Landau方程。在线性振幅方程中引入非线性项,考虑模态之间的耦合作用,这一模型被广泛用于流动稳定性分析中[2, 3]。然而,该模型仅考虑了模态自身的耦合作用,并未计及不同模态之间的联系。
随着对外流流动稳定性问题研究的深入,大量研究结果表明外流流动中存在着不同形式的模态作用。1996年William[3]对圆柱尾迹流动稳定性的理论和实验研究进行了综述,并指出尾迹的演化过程中存在2种不同的模态(Mode A与Mode B),分别对应流动发生二维涡脱落、流动出现三维不稳定性效应的流动稳定性现象。同年Barkley与Henderson[4]采用Floquet乘子分析了圆柱尾迹的稳定性。他们同样发现了对应长波长的模态A与涡脱落相关联,对应短波长的模态B与涡的二次失稳相关联。1999年Blackburn和Henderson[5]对振动圆柱绕流的尾迹结构进行了研究。他们发现在不同振动频率下流场具有特有的流场结构,其中包括周期模态,反对称周期模态等。进一步地,Govardhan与William[6]于2000年对自激振动圆柱绕流的尾迹斑图进行了细致研究。结果表明自激频率与流体涡脱落频率相关时流场出现独有的流动结构。在不同的自激频率下尾迹呈现出特有的涡对(P)与单涡(S)之间的耦合作用。近年来,proper orthogonal decomposition(POD)模态被引入到流动稳定性和控制的研究中。Noack等人[7, 8]在圆柱绕流的低维建模中提出了区别于常规模态的偏移模态,但这一模态是为了逼近原动力系统而人为引入的对平均流场的偏移修正。2010年Sengupta等人[2]对圆柱绕流的流体动力系统进行了模态分析。他们从流场拓扑结构角度发现了在圆柱尾迹中存在的3种模态:成对出现,具有周期性的常规模态;单一出现,也具有周期性,但非对称非常规模态;成对出现,呈现波特性,非周期性非常规模态。在尾迹中这3种模态的相互作用与尾迹的流动稳定性相联系。2011年Sengupta等人[9]指出在内流和外流中都存在具有相同性质的模态特征。这些模态的耦合作用与流动稳定性相关,并给出了多模态耦合的Landau-Stuart equation(LSE)模型。但是文中仅比较了利用模态获得的降维模型与direct numerical simulation(DNS)结果的接近程度,未对该模型进行深入研究。目前对模态间相互作用的分析主要集中在圆柱绕流研究上。对翼型绕流流场中特征模态的相互作用机理以及流动稳定性和流场典型模态间的联系并未得到广泛的关注,有待深入研究。
本文采用非线性动力学理论对翼型绕流的多模态耦合机制进行研究。首先,通过特征线有限元方法对翼型绕流问题进行数值计算,建立非定常流场数据库。利用POD理论提取大迎角下低雷诺数翼型绕流流场的特征拓扑结构,即POD模态。分析所得模态的特点以及模态间的耦合作用,并探索流动稳定性与POD模态之间的联系。
1 控制方程与数值算法 1.1 控制方程由于在本文中考虑流动马赫数小于0.3,流体可以认为是不可压缩流体。则无量纲黏性不可压流体的控制方程可以表示为
式中,u∞为自由来流速度,c为翼型的弦长,Re= ,v为运动黏性系数。 1.2 数值方法采用常规的有限元方法求解N-S方程时,由于在动量方程中的非线性对流项的作用,数值解会发生振荡。为解决这一问题,本文采用特征线有限元方法求解N-S方程。其思路是根据特征线的物理意义,建立沿特征线方向的新坐标系,来消去对流项。然后利用Taylor级数在空间坐标系上展开来得到抑制数值振荡的稳定项,最后采用有限元方法进行求解。特征线有限元characteristic based split(CBS)方法是由Zienkiewicz等人[10, 11]于1995年提出的用于求解可压缩和不可压流动问题的通用算法。但是该算法针对不可压缩流动时需要对压力进行隐式求解,这样就导致在求解格式上与可压流体的显式求解算法并不统一。2003年,Nithiarasu[12]通过引入人工压缩项改进该算法,应用于不可压流动的求解中,提高了算法的计算效率同时使得在求解格式上与可压流体算法一致。之后,康伟等人[13, 14, 15]通过引入arbitrary lagrangian eulerian(ALE)坐标系发展了CBS算法,使之应用于流-固耦合求解。该算法的具体实施步骤可参见文献[10, 11, 13, 14, 15]。
2 翼型绕流的低维建模 2.1 POD理论POD方法的目标是构造一组最优正交基{φi|i=1,…,∞},使得函数空间{qk∈L2(Ω)|k=1,…,M}投影到正交基上产生的误差最小[16],即
利用变分法,该最大值问题可以转化为一个特征值问题(3),求解该特征值问题,可以得到相应的特征向量,即为POD模态。
式中,C(x,y)=<u(x,t)⊗u(y,t)>为u的自相关函数,也称为POD的核,(φi,φj)=δij。由于具体应用中,流场瞬照都是在离散空间点上的实验数据或数值模拟结果。这时的自相关函数可以表示为对于计算流体力学问题,流场瞬照包含的空间点数Np远大于所选取的流场瞬照数M。如果直接求解特征值问题(3),待求解的矩阵阶数为Np×Np,这使得求解过程变得困难,且造成较大的误差。为了解决这一问题,本文采用了Sirovich提出的瞬照法[17],使得在流体动力系统中能够更高效快捷地提取POD模态集。
2.2 翼型绕流的低维模型在流体空间域Ω上,速度场所在的空间为Hilbert空间 的子空间,其内积可表示为
对于黏性不可压流体的流动控制方程(1),速度场可以在POD模态张成的空间中表示为
式中,ai(t)表示第i阶模态的时变系数,u0为速度时均场,a0≡1,ui,i≥1为正交的POD模态。需要指出的是传统POD模态提取方法是对线性化的Navier-Stokes方程进行求解得到瞬时流场,其模型忽略了流场内在的非线性因素。在本文中,由于对Navier-Stokes方程求解得到的瞬时流场作为POD数据库,提取得到POD模态,这样模态就包含了系统的非线性信息,并进行Galerkin投影得到系统的非线性降维模型,得到的模型保留了系统内在的非线性因素。将(6)式代入方程(1)可以得到流体动力系统在POD模态张成的空间中的降维动力系统模型(7)。对流体系统的流动稳定性研究转化为在模态张成的空间中,对应各模态的稳定性问题[19]。
其中线性项和二次项分别对应Navier-Stokes方程中的黏性项和对流项,且 3 多模态耦合与流动稳定性分析本文选取了Re=800下NACA0012翼型在迎角为20°时绕流问题进行研究。图1给出了Re=800的时均流场的涡量与流线。在Re=800时翼型近尾迹区具有清晰的涡结构。漩涡之间顺时针旋转的涡对应负涡量集中在上半部,逆时针旋转的涡对应正涡量集中分布在翼型的后缘附近,且占主导作用。
为了分析流场的特征拓扑结构,图2和图3给出了所得模态所占的能量分布、以及前4阶POD模态的流场拓扑结构。从模态所占的能量分布中可以看出POD经验在不同数量级的能量比例中成对出现。前2阶模态占流场能量比例最大,其中第1阶较第2阶POD模态所占流场总能量比例稍大,但仍分布在一个数量级上,第3阶POD模态与第4阶POD模态所占能量比例基本接近,其他的高阶模态也都如此。
从图3中可以看出成对出现的模态之间的流场拓扑关系。第1阶模态与第2阶模态在近翼面附近涡量分布相似。在尾迹中单涡与单涡正负交替排列出现(记为2S),这与卡门涡街结构相同。两模态在相同位置处涡量正负相反,这说明第1阶模态与第2阶模态间存在相位上的滞后。而第3阶模态与第4阶模态的涡量分布在近翼面附近涡量分布基本一致。这时涡结构较第1、2阶模态明显不同,在尾迹中涡的尺度有所减小,且出现涡对结构(记为P)。此外,第3、4模态之间同样存在相位差。
由POD模态得到的单涡与单涡(2S)、涡对(P)的流场拓扑结构是流场中特有的流动结构[6]。为了研究这些特征结构的相互作用如何影响翼型绕流的流动稳定性,本文分别选取前2、4、8、12阶模态求解不同模态数对应的动力系统(7),得到流场中各阶模态对应的系数ai变化情况。第i阶模态在流场的脉动能 为第i阶模态的系数。因此,模态系数ai的变化,反映了该阶模态在流场中脉动能的变化情况。
由于POD模态成对出现,图4成对给出了不同模态数时各阶模态系数幅值的时间历程。前2阶模态是流场中占主导地位的模态,因为它们占流场总能量的50%以上(见图2)。当仅选取前2阶模态表示流场时,从图4a)可以看出系数的幅值一直线性增大,当时间增加至143 121时,幅值发生突然下降,之后迅速增大,直至计算溢出。利用分岔理论分析由前2阶模态构成的动力系统(7),该系统的Jacobian矩阵的特征值的实部为0.267 562×10-0.3>0,这说明该模态对应的系统是不稳定的。当选取前4阶模态表示流场时,前2阶模态的幅值最终稳定至439.227,这时前2阶模态的系数随时间周期性变化,且存在约π/2的相位差(图5实线),所对应能量周期脉动。这说明引入第3、4阶模态后,前2阶模态的稳定性发生了改变。由于模态之间的耦合作用,流场的能量从前2阶模态中部分地转移到第3、4阶模态,使得前2阶模态的脉动能变得稳定。当考虑前8阶和前12阶模态时,前2阶模态的幅值分别收敛到49.336,33.12。随着模态数的增加,振幅逐渐收敛到一个稳定值。这说明更多模态的相互作用,使得系统的低阶模态的能量向高阶模态转移,同时由于非线性耦合作用,能量的转移逐渐趋于饱和,最终达到稳定。另外,前2阶模态系数稳定性的改变也表明流场中2S涡结构对扰动极其敏感。2S涡结构在流场中是不稳定的流动结构,在高阶模态作用下该流场结构的稳定性会发生改变。
对于第3、4阶模态,当取模态数为4时,第3、4阶模态的振幅随时间发生周期性变化,如图4b)所示。对应的模态系数出现多个频率叠加的复杂运动(图5虚线)。当模态数进一步增加时,对比图4b)的粗线对应的振幅,可以看出第3、4阶模态的振幅(见图4b)的细线与粗线)的大小有所变化,但振动形式并未发生实质性改变。第3、4阶模态对应了流场中的P涡结构。模态系数的振动形式对高阶模态作用未发生实质性的改变,说明流场中的P涡结构在流场中不易受到扰动的影响,是流场中稳定的流动结构。
图6给出了在Re=800附近以雷诺数为分岔参数的翼型绕流动力系统分岔图。通过计算可以得到当雷诺数小于624.5时,系统的Jacobian矩阵在平衡位置处的特征根的实部都是负的,这时系统的平衡位置渐进稳定,对应的模态系数最终收敛到一点上。当雷诺数等于624.4时,系统的Jacobian矩阵在平衡位置处存在1对纯虚特征根,且其他特征根的实部不等于零。这时系统发生Hopf分岔。系统的平衡位置失稳(对应图6的标记点H),系统变为周期运动。当雷诺数大于624.4时,系统的Jacobian矩阵在平衡位置处的特征根实部均大于零,这时平衡位置是不稳定性的(对应图6的虚线部分)。从相图7中可以看出在雷诺数为Re=800时系统为周期运动。
4 结 论本文通过对翼型绕流的POD模态耦合作用的分析,得到如下结论:
1)在翼型绕流流场中存在具有典型特征的涡结构,如2S涡结构,P涡结构。
2)在流场中模态之间存在非线性的耦合关系。流场的能量能够从低阶模态向高阶模态转移,起到稳定流场作用,但最终能量转移达到饱和。
3)2S涡结构在流场中是不稳定的流动结构,对扰动极其敏感。该涡结构在高阶模态作用下其稳定性会发生改变;流场中的P涡结构在流场中不易受到扰动的影响,是流场中稳定的流动结构。
4)翼型绕流动力系统在雷诺数为624.4时,系统的Jacobian矩阵在平衡位置处存在1对纯虚特征根,且其他特征根的实部不等于零。这时系统发生Hopf分岔。系统的平衡位置失稳,变为周期运动。
[1] | Stuart J T. On the Non-Linear Mechanics of Wave Disturbances in Stable and Unstable Parallel Flows Part 1: The Basic Behaviour in Plane Poiseuille Flow[J]. Journal of Fluid Mechanics, 1960, 9(3): 353-370 |
Click to display the text | |
[2] | Sengupta T K, Singh N, Suman V. Dynamical System Approach to Instability of Flow Past a Circular Cylinder[J]. Journal of Fluid Mechanics, 2010, 656(82): 115 |
Click to display the text | |
[3] | Williamson C. Vortex Dynamics in the Cylinder Wake[J]. Annual Review of Fluid Mechanics, 1996, 28(1): 477-539 |
Click to display the text | |
[4] | Barkley D, Henderson R D. Three-Dimensional Floquet Stability Analysis of the Wake of a Circular Cylinder[J]. Journal of Fluid Mechanics, 1996, 322: 215-242 |
Click to display the text | |
[5] | Blackburn H, Henderson R. A Study of Two-Dimensional Flow Past an Oscillating Cylinder[J]. Journal of Fluid Mechanics, 1999, 385(1): 255-286 |
Click to display the text | |
[6] | Govardhan R, Williamson C. Modes of Vortex Formation and Frequency Response of a Freely Vibrating Cylinder[J]. Journal of Fluid Mechanics, 2000, 420(85): 130 |
Click to display the text | |
[7] | Noack B R, Afanasiev K, Morzynski M, et al. A Hierarchy of Low-Dimensional Models for the Transient and Post-Transient Cylinder Wake[J]. Journal of Fluid Mechanics, 2003, 497(1): 335-363 |
Click to display the text | |
[8] | Noack B R, Tadmor G, Morzynski M. Low-Dimensional Models for Feedback flow Control. Part I: Empirical Galerkin Models[R]. AIAA-2004-2408 |
Click to display the text | |
[9] | Sengupta T K, Vijay V, Singh N. Universal Instability Modes in Internal and External Flows[J]. Computers & Fluids, 2011, 40(1): 221-235 |
Click to display the text | |
[10] | Zienkiewicz O, Codina R. A General Algorithm for Compressible and Incompressible Flow——PartⅠ the Split, Characteristic-Based Scheme[J]. International Journal for Numerical Methods in Fluids, 1995, 20(8/9): 869-885 |
Click to display the text | |
[11] | Zienkiewicz O, Morgan K, Sai B, et al. A General Algorithm for Compressible and Incompressible Flow——PartⅡ: Tests on the Explicit Form[J]. International Journal for Numerical Methods in Fluids, 1995, 20(8/9): 887-913 |
Click to display the text | |
[12] | Nithiarasu P. An Efficient Artificial Compressibility (AC) Scheme Based on the Characteristic Based Split (CBS) Method for Incompressible Flows[J]. International Journal for Numerical Methods in Engineering, 2003, 56(13): 1815-1845 |
Click to display the text | |
[13] | 康伟, 张家忠. 翼型局部弹性自激振动的增升减阻效应研究[J]. 西安交通大学学报, 2011(05): 94-101 Kang W, Zhang J Z. Numerical Analysis of Ligt Enhancment and Drag Reduction by Self-Induced Vibration of Localized Elastic Airfoil[J]. Journal of Xi'an Jiaotong University, 2011(05):94-101 (in Chinese) |
Cited By in Cnki (9) | |
[14] | Kang W, Zhang J Z, Feng P H. Aerodynamic Analysis of a Localized Flexible Airfoil at Low Reynolds Numbers[J]. Communications in Computational Physics, 2012, 11(4): 1300-1310 |
[15] | Kang W, Zhang J Z, Lei P F, et al. Computation of Unsteady Viscous Flow around a Locally Flexible Airfoil at Low Reynolds Number[J]. Journal of Fluids and Structures, 2014, 46: 42-58 |
Click to display the text | |
[16] | Holmes P, Lumley J L, Berkooz G. Turbulence, Coherent Structures, Dynamical Systems and Symmetry[M]: Cambridge Univ Pr, 1998 |
[17] | Sirovich L. Turbulence and the Dynamics of Coherent Structures. Part I: Coherent Structures[J]. Quarterly of Applied Mathematics, 1987, 45 (3): 561-571 |
Click to display the text | |
[18] | Kang W, Zhang J Z, Ren S, et al. Nonlinear Galerkin Method for Low-Dimensional Modeling of Fluid Dynamic System Using POD Modes[J]. Communications in Nonlinear Science and Numerical Simulation, 2015, 22(1/2/3): 943-952 |
Click to display the text |
2. School of Aeronautics and Astronautics Engineering, Air Force Engineering University, Xi'an 710043, China