用于跨/超声速壁板颤振精确分析的流-固耦合有限元算法

梅冠华,杨树华,张家忠,孙旭,陈嘉辉

(1.西安交通大学能源与动力工程学院, 710049, 西安; 2.沈阳鼓风集团股份有限公司, 110869, 沈阳)

用于跨/超声速壁板颤振精确分析的流-固耦合有限元算法

梅冠华1,杨树华2,张家忠1,孙旭1,陈嘉辉1

(1.西安交通大学能源与动力工程学院, 710049, 西安; 2.沈阳鼓风集团股份有限公司, 110869, 沈阳)

为了精确和定量分析超声速与跨声速壁板的颤振特性,提出了一种基于有限元方法的流-固耦合算法,并用其研究了二维壁板颤振问题。首先,给出了壁板的von Kármán几何大变形运动方程,以及高速气流的欧拉控制方程。然后,采用标准有限元方法对壁板方程进行空间离散,而对流动控制方程的离散则运用双时间步长推进的特征线分裂有限元方法,从而有效地消除了流场数值解的振荡问题。随后,采取松耦合算法实现了流体与固体间的数据传递。最后,运用所提出的算法对超声速和跨声速气流作用下壁板的气动弹性特性进行了分析,考察了归一化动压、预紧力和厚度比对系统特性的影响,并将该算法的分析结果与采用线性/非线性活塞理论和线性化势流理论的经典壁板颤振结果进行了对比,证明该算法可以在较宽广的马赫数范围内给出气动力的精确描述,尤其适合于分析跨声速气流下的壁板气动弹性响应。

壁板颤振;流-固耦合;特征线分裂算法;有限元方法;气动弹性

在特定参数下,飞行器表面暴露于高速气流中的壁板结构将产生自激振动现象,常称作壁板颤振。这类振动是典型的气动弹性问题,往往是由惯性力、弹性力和气动力的共同作用而激发的,通常具有很强的非线性动力学特性,并将严重影响飞行器的疲劳寿命、飞行性能、飞行安全和乘坐品质。因此,深入研究壁板颤振问题对于高速飞行器的壁板设计、颤振抑制、疲劳寿命估计和颤振的合理利用都具有十分重要的意义[1-3]。

作为经典问题,壁板颤振已由众多学者进行了深入细致的研究[4-10]。在大多数壁板颤振的研究中,壁板一般采用von Kármán几何大变形理论描述,以便准确捕捉系统的各类复杂响应,气动载荷大多采用简化气动力模型近似表达,如线性/非线性活塞理论[4,11]、线性化势流理论[5-6]等,由此可推导出系统的偏微分控制方程。对于此类控制方程,常采用Galerkin方法、Rayleigh-Ritz方法、有限元方法(FEM)等将其离散为常微分控制方程,并运用各种频域或时域分析工具进行研究。虽然这些简化气动力理论使用起来较为方便,然而其适用范围和精度都较为有限,一般而言,线性活塞理论适用于1.45时,由于非线性效应变得十分显著,气动力常采用非线性活塞理论来描述。目前,尚没有一个普遍适用于亚声速、跨声速、超声速和高超声速流动的简化气动力理论。此外,简化气动力理论给出的气动载荷是与固体位移、速度及其偏导数相关的函数,而并非是通过计算流体动力学(CFD)所得出的,因此无法给出流场的细节信息,如黏性边界层、旋涡、流动分离、激波等,进而无法准确分析流体与固体间的耦合机理。事实上,这些正是壁板颤振定量分析中的开放性问题。

近年来,伴随着计算机软、硬件水平的飞速发展,流-固耦合算法被逐渐用于分析壁板颤振问题。其中,气动力由基于流体欧拉或Navier-Stokes控制方程的CFD方法获得,同时将CFD与计算结构动力学(CSD)耦合起来研究壁板的气动弹性响应。流-固耦合算法不仅广泛适用于亚声速、跨声速、超声速和高超声速流动,可以在几乎所有马赫数范围内给出气动力的精确表达,而且还可以给出流场的细节描述,这些都是简化气动力理论所不能实现的。使用流-固耦合算法,Davis求解了二维壁板颤振问题,发现了许多不同于简化气动力理论的新现象[12-13];随后,Gordiner研究了二维和三维壁板颤振问题,发现边界层的存在推迟了颤振的发生[14];Hashimoto分析了跨声速三维壁板颤振问题,在与他人实验结果取得一致的同时,还发现湍流边界层不仅对颤振可以起到稳定作用,在特定条件下还能引起不稳定效应[15]。

在目前的壁板颤振的流-固耦合研究中,对流场的求解大多采用有限差分法(FDM)和有限体积法(FVM),而采用FEM的较少。虽然FEM的推导过程相对烦琐,计算时间也稍长,然而FEM有其独特的优势:①由于无须构造交错网格,FEM的网格生成和数据结构较为简单;②在已经生成的粗网格上,FEM可以通过选用高阶单元来提高精度,而FDM和FVM则往往需要通过加密网格来提高精度;③由于FEM的普适性,可被广泛用来处理多场耦合问题,尤其是流-固耦合问题。若能将流体和固体采用统一的FEM求解,就有可能实现流-固同步耦合求解,并在后续研究中引入系统自由度缩减方法[16],不仅能加快求解速度,更能提取出符合物理意义的气动弹性模态信息。

由于流体控制方程中的非线性对流项往往会引起数值解的振荡,因此标准Galerkin FEM不能直接用于求解流体控制方程。基于此,众多改进的FEM被不断提出,如流线迎风Petrov-Galerkin(SUPG)方法、Taylor-Galerkin(TG)方法、Galerkin最小二乘(GLS)方法等,其中特征线分裂(CBS)方法由Zienkiewicz首先提出,并用于求解对流扩散方程[17]。经过众多学者的不断发展和完善,迄今为止,特征线分裂有限元方法(CBS-FEM)已被广泛应用于求解各类流动问题,例如不可压缩流动、可压缩流动、含激波的高速流动、湍流流动等。近年来,将CBS-FEM发展应用于求解动边界问题,如流-固耦合问题和自由液面流动问题等,已成为新的研究热点。

为了精确和定量分析超声速与跨声速壁板颤振特性,本文提出一种基于FEM的流-固耦合算法,并用来在时域内对二维壁板颤振问题进行分析。首先,分别采用von Kármán几何大变形方程和欧拉方程来描述壁板变形和高速气流。然后,对壁板控制方程采用标准FEM进行空间离散,对流体控制方程则运用双时间步推进的CBS-FEM进行离散。随后,采取松耦合算法实现流体与固体间的双向耦合。最后,运用所提出的算法对超声速和跨声速气流下壁板的气动弹性响应进行分析,考察归一化动压、预紧力和厚度比对系统特性的影响,并对壁板的稳定性边界进行计算。

1 壁板颤振模型及其数值解法

若三维壁板的展向尺寸远大于其弦长,则可忽略展向效应,将其简化为二维壁板,这给问题的分析带来了很大便利。作为初步研究,先对二维壁板颤振问题进行分析,其模型如图1所示:二维平板置于刚性平面上,其上表面处于沿x方向的高速来流中,下表面对应着空腔。流体的来流速度、马赫数和密度分别为U∞、Ma∞和ρ∞,平板的长度、厚度、单位长度质量、弹性模量和泊松比分别为a、h、ρm、E和μ,平板两端初始拉伸面内力为N0,上表面由流场所施加的气动载荷为p,下表面空腔压力与来流压力相同,皆为p∞。

图1 二维平面壁板颤振示意图

对于厚度远小于长度的壁板,即常见的薄板,可以采用Kirchhoff-Love假设,仅考虑其横向运动w(x,t)。应变-位移关系采用von Kármán几何大变形理论表示,则平板的控制方程为[4]

(1)

其中平板的弯曲刚度D=Eh3/[12(1-μ2)],大变形产生的中面拉伸载荷

(2)

气动载荷

Δp=p-p∞

(3)

壁板的边界条件可以为简支或固支。对于简支端,其数学表达式为w=0,∂2w/∂x2=0;对于固支端,其数学表达式为w=0,θ=∂w/∂x=0。

N1=1-3ξ2+2ξ3;N2=l(ξ-2ξ2+ξ3)

N3=3ξ2-2ξ3;N4=l(ξ3-ξ2)

采用标准Galerkin FEM对控制方程(1)进行变分以及空间离散,将其转化为关于时间二阶导数的常微分方程组

(4)

各单元矩阵的具体表达式为

q={w1,θ1,w2,θ2,…,wN+1,θN+1}T为壁板的整体位移。

对于该动力系统,只要给定初始条件,采用自适应Runge-Kutta积分方法数值求解即可。

2 流动控制方程及其数值解法

在高速流动问题中,黏性的影响相对较小,而可压缩性的影响十分显著,故为方便研究,可以采用可压缩无黏Navier-Stokes方程,即欧拉方程来描述流场,其守恒形式为

连续方程

(5)

动量方程

(6)

能量方程

(7)

补充方程

(8)

式中:ρ为流体密度;p为流体压强;e为单位质量流体的内能;Q为单位质量流体的总能量;γ=1.4为空气的比热容比;xi、ui和Ui=ρui分别为坐标、速度和守恒速度。在二维流动中,下角标i=1,2分别对应着x和y方向。

流动方程的求解采用基于特征线分裂算法的有限元方法,将流场变量沿着特征线运动方向进行变换处理,从而可消除对流项,并有效避免数值解的振荡。当前,在可压缩流动问题求解中广泛采用的是单时间步长推进的CBS-FEM,由于其大多采用完全显式格式,使得时间步长的选取受限于收敛条件。此外,为方便计算,在求解过程中常进行质量集中处理,这种简化对稳态问题没有任何影响,因为随着稳态解的收敛,时间变化项也会随之消失。然而,对于瞬态问题而言,该近似所带来的误差相当严重,此时为获得准确解,可引入附加的迭代过程。基于此,本文采用双时间步长CBS-FEM来求解欧拉方程。事实上这是隐式格式算法,在每个物理时间步的推进过程中,通过引入一系列伪时间步的迭代过程,来获得局部的伪定常收敛解[17-18]。这样,真实时间步长的选取就不再受限于收敛条件,而且质量集中带来的误差也可有效消除。双时间步长CBS-FEM算法的具体流程如下。

第1步

(9)

第2步

(10)

第3步

(11)

第4步

(12)

修正速度、密度和能量

(ρQ)m+1=(ρQ)m+Δ(ρQ)

计算压力

(13)

式(9)～式(13)中:上角标n和n+1分别代表当前时刻和下一时刻,m和m+1分别对应前、后2个伪时刻;Δt为物理时间步长;Δτ为伪时间步长;θ1和θ2为松弛因子,一般选取θ1=1,θ2=0。

采用线性三角形单元剖分计算区域,则流体变量可表达为形函数和节点值相乘的形式

运用标准Galerkin有限元方法,将式(9)～式(12)进行变分和空间离散,并忽略高阶项,即可得到双时间步长推进CBS-FEM的离散格式。

一般而言,对于存在激波的流动问题,如果直接进行数值模拟,则解在激波附近会引起强烈的振荡,影响收敛性。因此,必须添加合理的人工黏性项,选用合适的激波捕捉方法来抑制振荡并捕捉激波。本文采用基于压力二阶导数的激波捕捉方法[17,19]。

壁板上方流场计算区域及网格剖分如图2和图3所示,远场边界位于25倍壁板弦长处。为了生成高质量网格并获得精确的气动载荷,将计算区域划分为2个子区域:子域1为壁板上方0.5a高度的矩形区域,其余为子域2。对于子域1,流-固耦合界面和高度方向分别均布48和24个网格;对子域2,布置环绕子域1的48层C型网格。流场网格共计5 881个节点和11 520个三角形单元。

图2 壁板上方流场计算区域(非等比例绘制)

(a)整体网格

(b)局部网格

远场边界条件由Riemann无反射边界条件给定[12]。壁板上游和下游皆为刚性壁面,在无黏流动中将其法向速度设为0即可。流-固交界面虽然也为壁面边界条件,然而由于壁板运动的影响,将其给定为(u-us)·n=0,其中u为流体速度,us为壁板运动速度,n为壁板节点的法向量。

3 动网格及流-固耦合方式

随着壁板的上下运动,流-固交界面也在不断运动。在该类动边界流动问题的计算过程中,如果规定动边界上的节点随边界同步运动,将使得动边界附近的网格单元发生扭曲变形甚至重叠,导致无法计算流场;如果规定流场网格固定不变,不随边界运动而移动,就会给边界位置的捕捉带来很大不便。常用的处理方法是采用动网格技术,让动边界上的网格节点随着边界同步运动,求解区域内部的节点也随之相应调整,从而不仅可以准确地描述动边界位置,还可以保证流体网格的质量。

由于该流-固耦合问题的特殊性——壁板仅在y方向产生位移,而且在壁板上方的子域1内均布了24层结构化流体网格,故仅需将壁板节点的位移在子域1内沿着y方向均分即可,这样既能保证网格质量,又不会在动网格上耗费大量的计算时间。作为验证,令壁板以w=0.05asin(2πx/a)形式弯曲,移动后的网格如图4所示。整体和子域1内网格的最大单元偏斜度分别为0.568 2和0.485 7,远小于0.9,说明虽然壁板发生了较大的变形,然而网格的质量仍得以较好地保持,这证实了动网格方法的可靠性。

图4 动网格方法验证

(14)

式中:Δxi为旧网格节点至新网格节点的位移。

流-固耦合的方式为松耦合,在每个时间步流体与固体相互交换一次数据,也称为单步耦合,具体实施方法如下:

(1)将气动载荷Δp传递给壁板;

(2)将壁板响应推进至下一时刻;

(4)移动流体网格;

(5)更新流场初值;

(6)将流场变量推进到下一时刻;

(7)继续下一个时间步直到计算结束。

与紧耦合和同步耦合算法相比,松耦合使用起来较为方便,仅需对CFD和CSD的程序进行适当修改即可,同时它的计算速度较快,是一种应用较广的双向耦合算法。然而,由于在一个时间步长内流体和固体仅相互交换一次数据,且该数据并非是同一时刻的,因此两者间存在着时间滞后,如果时间步长选取不当,很容易造成误差累积并导致解的发散,故需要慎重选择时间步长。经过不断尝试,最终将时间步长选取为ΔtU∞/a=0.01,这样既能保证解的准确性,又不至于消耗过长的计算时间。对于该时间步长的独立性验证,将在下一节进行阐述。

4 数值模拟和结果分析

4.1 网格无关性分析和时间步长独立性验证

先分析流场网格的无关性,计算Ma∞=2时壁板在固定变形位置y=0.1asin(πx/a)上的稳态扰流问题。计算3种不同疏密的网格,即第2节给出的标准网格、疏网格、密网格,在其壁板边界分别均布48、24和72个单元,得到壁面上的相对压力分布情况,如图5所示,可见在壁板前缘和后缘存在着2个明显的激波。标准网格与密网格所得出的结果并无明显差异,说明采用标准网格可以得到足够精度的流场。

图5 流场网格无关性分析

图6 固体网格无关性分析

图7 流-固耦合算法中时间步长的独立性分析

4.2 超声速气流下的壁板响应

当归一化预紧力R=0时,在较小的归一化动压λ下,初始扰动系统经过瞬态振动将恢复到初始位置,表现为平面稳定状态。如图8所示λ=300时的壁板响应,经过大约30至40个周期的瞬态振动,初始扰动被完全衰减,系统恢复到了平面稳定状态。当λ逐渐增大并越过临界值λcr,系统将发生Hopf分岔,由平面稳定状态转变为极限环颤振状态。λ=500时的系统响应如图9a所示,稳态极限环颤振的最大值和最小值分别为0.660和-0.679,后者的绝对值略大于前者,表明该振动关于中性面是非对称的,这意味着壁板的正向运动和负向运动对于流场的影响是不同的。

作为该振动的详细分析,图9b和图9c分别给出了在一个周期内不同时刻的壁板变形和气动载荷。从图9b中可以看出,虽然不同空间质点的振幅不同,然而各质点在时间上的运动规律是相同的,即它们同时达到最大或最小位移,这是驻波的典型特征,因此可认为壁板振动呈现出了驻波的形式,其中最大振幅位于x=0.75a附近,最小振幅位于x=0.3a附近。图9c所示的压力波动也呈现出驻波的运动形式,压力最大值约位于壁板末端x=a附近,且随着壁板的上下运动在壁板尾缘附近表现为激波与膨胀波的交替变换。由于激波是流场中典型的非线性现象,在激波附近的流体被强烈压缩,因此激波的存在对壁板的颤振具有重大影响。

(a)位移时间历程

(b)一个周期内的瞬时壁板位移

(c)一个周期内的瞬时壁板表面压力分布

为考察系统对归一化动压的敏感性,研究了λ对极限环颤振的影响,结果如图10所示。由图中可以看出,在λ>λcr时壁板发生颤振,振幅与壁板厚度量级相同,且随λ的增大而增大。图10中还展示了线性和非线性活塞理论所得的结果[21]:线性活塞理论所得出的振动正负对称,而非线性活塞理论所得振动的负向峰值绝对值大于正向峰值,这与本文算法的结果是一致的。定量来说,本文算法所得出的临界动压稍大,而振幅稍小。由于活塞理论没有考察壁板不同时刻和不同空间质点间的相互影响,是一种基于当地流动和当前时刻的简化气动力理论,而基于欧拉方程的流-固耦合算法则直接从流体控制方程着手,综合考虑了上述因素,利用CFD方法可以精确计算气动载荷,因此这是两者间存在差异的主要原因。

(a)位移时间历程

(b)屈曲变形和压力分布

在施加了预紧力的情况下,壁板更多形式的响应将被激发出来。在较小的λ下,当归一化预紧力R逐渐增大并跨过临界值时,系统将发生静态分岔,从平面稳定状态转变为屈曲状态。例如:当λ=50和R=2时,壁板的位移时间历程如图11a所示,在动压和预紧力的联合作用下,系统最终到达稳定屈曲变形位置;壁板的最终变形和表面压力分布如图11b所示,最大变形位置大约位于x=0.6a附近,壁板的前缘和尾缘存在2个明显的激波,其中尾缘处的激波比前缘处的激波更为剧烈,这可能是屈曲最大变形位置略向后偏移的结果。

(a)位移时间历程

(b)相图

(c)频谱分析

在特定的λ和R下,系统可表现出非线性复杂振动。例如:在λ=140和R=5时,壁板的稳态响应如图12a所示,虽然振动为周期形式,然而振动方式却表现出了极强的非线性特性;在图12b所示的相图中,不难发现该振动具有倍周期运动的特性;在图12c所示的该振动的频谱中包含了6阶频率分量,分别约为15、30、45、60、75和90 Hz,其中前3阶频率所占比重较大,后3阶频率的比重较小,由于频率间呈现出1至6的整数倍关系,因此系统表现出了倍周期运动的特征。

在R和λ所构成的参数平面上,对系统的平面稳定、屈曲、振动响应区域进行了划分,如图13所示:随着R从0开始逐渐增大,壁板的颤振临界动压λcr逐渐降低,即预紧力的增加使得壁板更易发生颤振;随着λ从0开始逐渐增大,壁板的屈曲临界预紧力R逐渐增大,即气动载荷的增加使得壁板不易发生屈曲。此外,本文算法所得结果与线性活塞理论所得结果是一致的。

4.3 跨声速气流下的壁板响应

归一化动压λ*对于颤振幅值的影响如图15所示:当λ*越过临界值后,最大振幅随着λ*的增加而增大,本文算法所得结果与文献[5]采用线性化势流理论所得结果吻合较好,然而随着λ*的增大,两者间的差距逐渐增大,这可能是由于线性化势流理论是从扰动速度和势流理论的角度来推导气动力的,并不能准确反映整体流场的状况,与欧拉方程算得的气动力仍存在一定差距的缘故。

对于Ma∞<1(如Ma∞=0.8和Ma∞=0.9)的跨声速气流,系统失稳后壁板将偏离初始位置,根据初始条件的不同到达正向或是负向屈曲变形位置。以位于海平面的铝板为例,当Ma∞=0.9和h/a=0.004时,壁板的屈曲响应如图16所示,可以看到其正向屈曲变形略大于负向屈曲变形。

关于系统对厚度比的敏感性,图17给出了在Ma∞=0.8和Ma∞=0.9时h/a对壁板正向和负向屈曲变形的影响:最大屈曲变形随着h/a的增大而减小,这是因为较小的h/a意味着壁板相对刚度较小,更容易发生屈曲;Ma∞=0.8时正、负屈曲变形近乎对称,而Ma∞=0.9时正向变形要大于负向变形,这说明Ma越接近1,系统变得越复杂。此外,本文算法所得结果与文献[13]的结果相当一致。

图18显示了Ma∞=0.8～1.4时跨声速气流作用下壁板的稳定性边界。在Ma∞=1.0的声速气流附近,临界动压明显降低,稳定性曲线在此处形成了一个“凹坑”,通常称为“跨声速凹坑”,因此对壁板颤振的研究设计者来说,声速附近的临界动压是需要仔细考察的。本文算法所获得的稳定性边界与文献[13]的结果亦相当一致。

(a)正向屈曲

(b)负向屈曲Ma∞=0.9, h/a=0.004, 海平面铝板

(a)Ma∞=0.8

(b)Ma∞=0.9

图18 壁板的稳定性边界

5 结论与展望

为了精确定量分析超声速和跨声速壁板的颤振特性,本文提出了一种基于有限元方法的流-固耦合算法,并用其研究了二维壁板颤振问题。该算法成功消除了流体方程中对流项所引发的数值振荡,并实现了流体和固体间的双向数据交换。通过对数值结果进行分析,所得主要结论如下:

(1)在马赫数为2的超声速气流下,当归一化动压越过临界值后,系统由平面稳定状态转变为极限环颤振状态,振幅随动压增大而增大,且负向峰值的绝对值略大于正向峰值。对典型极限环颤振的分析发现:一个周期内壁板变形和气动载荷皆呈现驻波的运动形式;在归一化预紧力与归一化动压的共同作用下,系统可表现为屈曲或非简谐复杂振动;激波对壁板颤振特性的影响较为显著。

(2)在马赫数为1.2的跨声速气流下,壁板响应与马赫数为2的系统类似,只是此时系统对外界扰动反应较快,达到稳态响应所需的瞬态过程较短。

(3)在马赫数为0.8和0.9的亚声速气流下,壁板失稳后表现为正向或负向的屈曲变形,正、负位置由初始条件所决定,且屈曲变形随着壁板厚度比的减小而增大。

(4)对马赫数在0.8至1.4之间的系统的稳定性边界进行了计算,精确捕捉到了“跨声速凹坑”现象。

(5)通过与采用线性/非线性活塞理论和线性化势流理论的经典壁板颤振结果进行对比,表明本文算法可在更为宽广的马赫数范围内给出气动载荷的精确描述,尤其适合于分析跨声速气流下的壁板气动弹性特性。

作为初步探索,本文对流动采用了欧拉方程描述,下一步将采用完全Navier-Stokes方程,并考察湍流的影响,以更加准确地描述流场和气动载荷。此外,流-固耦合的计算过程耗费了大量的计算时间,未来将引入有效的系统自由度缩减方法,在保证计算精度的同时实现节省计算时间的目的。本文的工作仅对固体的气动弹性特性进行了分析,今后将对流动的稳定性与失稳机理展开深入研究。

[1] 梅冠华, 张家忠. 时滞惯性流形在三维壁板颤振数值分析中的应用 [J]. 西安交通大学学报, 2011, 45(9): 40-46. MEI Guanhua, ZHANG Jiazhong. Numerical analysis of 3-D panel flutter by inertial manifolds with delay [J]. Journal of Xi’an Jiaotong University, 2011, 45(9): 40-46.

[2] 杨智春, 夏巍, 孙浩. 高速飞行器壁板颤振的分析模型和分析方法 [J]. 应用力学学报, 2006, 23(4): 537-542. YANG Zhichun, XIA Wei, SUN Hao. Analysis of panel flutter in high speed flight vehicles [J]. Chinese Journal of Applied Mechanics, 2006, 23(4): 537-542.

[3] MEI Guanhua, ZHANG Jiazhong, WANG Zhuopu. Numerical analysis of panel flutter on inertial manifolds with delay [J]. Journal of Computational and Nonlinear Dynamics, 2013, 8(2): 021009.1-021009.11.

[4] DOWELL E H. Nonlinear oscillations of a fluttering plate [J]. AIAA Journal, 1966, 4(7): 1267-1275.

[5] DOWELL E H. Nonlinear oscillations of a fluttering plate: II [J]. AIAA Journal, 1967, 5(10): 1856-1862.

[6] DOWELL E H. A review of the aeroelastic stability of plates and shells [J]. AIAA Journal, 1970, 8(1): 385-399.

[7] OLSON M D. Finite element approach to panel flutter [J]. AIAA Journal, 1967, 5(12): 226-227.

[8] OLSON M D. Some flutter solutions using finite element [J]. AIAA Journal, 1970, 8(4): 747-752.

[9] MEI Chuh, ABDEL-MOTAGLY K, CHEN R. Review of nonlinear panel flutter at supersonic and hypersonic speeds [J]. ASME Applied Mechanics Reviews, 1999, 52(10): 321-332.

[10]CHENG Guangfeng, MEI Chuh. Finite element modal formulation for hypersonic panel flutter analysis with thermal effects [J]. AIAA Journal, 2004, 42(4): 687-695.

[11]ASHLEY H, ZARTARIAN G. Piston theory: a new aerodynamic tool for the aeroelastician [J]. Journal of the Aeronautical Science, 1956, 23(12): 1109-1118.

[12]DAVIS G A, BENDIKSEN O O. Unsteady transonic two-dimensional Euler solutions using finite elements [J]. AIAA Journal, 1993, 31(6): 1051-1059.

[13]DAVIS G A. Transonic aeroelasticity solutions using finite elements in an arbitrary Lagrangian-Eulerian formulation [D]. Los Angeles, USA: University of California, 1994.

[14]GORDINER R E, FITHEN R. Coupling of a nonlinear finite element structural method with a Navier-Stokes solver [J]. Computers and Structures, 2003, 81(2): 75-89.

[15]HASHIMOTO A, AOYAMA T. Effects of turbulent boundary layer on panel flutter [J]. AIAA Journal, 2009, 47(12): 2785-2791.

[16]ZHANG Jiazhong, REN Sheng, MEI Guanhua. Model reduction on inertial manifolds for N-S equations approached by multilevel finite element method [J]. Communications in Nonlinear Science and Numerical Simulation, 2011, 16(1): 195-205.

[17]ZIENKIEWICZ O C, TAYLOR R L, NITHIARASU P. The finite element method for fluid dynamics [M]. 6th ed. Singapore: Elsevier Pte Ltd., 2009: 195-211.

[18]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.

[19]NITHIARASU P, ZIENKIEWICZ O C, SATYASAI B V K, et al. Shock capturing viscosities for the general fluid mechanics algorithm [J]. International Journal for Numerical Methods in Fluids, 1998, 28(9): 1325-1353.

[20]孙旭, 张家忠. 具有运动边界的不可压缩黏性流动的CBS有限元解法 [J]. 西安交通大学学报, 2011, 45(1): 99-104. SUN Xu, ZHANG Jiazhong. A characteristic-based split-FEM scheme for incompressible viscous flow with moving boundaries [J]. Journal of Xi’an Jiaotong University, 2011, 45(1): 99-104.

[21]EASTEP F E, MCINTOSH S C. Analysis of nonlinear panel flutter and response under random excitation or nonlinear aerodynamic loading [J]. AIAA Journal, 1971, 9(3): 411-418.

(编辑 葛赵青)

AFluid-StructureCouplingAlgorithmBasedonFiniteElementMethodforPreciseAnalysisofTransonicandSupersonicPanelFlutter

MEI Guanhua1,YANG Shuhua2,ZHANG Jiazhong1,SUN Xu1,CHEN Jiahui1

(1. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2. Shenyang Blower Works Group Corporation, Shenyang 110869, China)

To analyze the supersonic and transonic panel flutter behavior quantitatively and accurately, a fluid-structure coupling algorithm based on the finite element method (FEM) is proposed for the two-dimensional panel flutter problem. First, the von Kármán’s large deformation theory is adopted to model the panel, and the high speed air flow is approached by the Euler equations. Then, the equation of panel is discretized spatially by the standard FEM, and the equations of fluid are discretized by the characteristic-based split finite element method (CBS-FEM) with dual time stepping, thus the numerical oscillation often encountered in numerical simulation of fluid flow can be eliminated. Furthermore, a loose coupling algorithm is applied to the data exchange between the fluid and the structure. Finally, the proposed algorithm is used to investigate the aeroelastic behavior of the panel in supersonic and transonic air flows and the influences of the non-dimensional dynamic pressure, pre-tightening force and thickness ratio on the system. The results are compared with those of the classical panel flutter analyses using linear/nonlinear piston theory and linearized potential flow theory. It shows that the proposed algorithm enables to obtain accurate aerodynamic pressure in a wide range of Mach numbers, especially for the analysis of panel aeroelasticity in transonic air flows.

panel flutter; fluid-structure interaction; characteristic-based split method; finite element method; aeroelasticity

10.7652/xjtuxb201401013

2013-03-29。 作者简介: 梅冠华(1984—),男,博士生;张家忠(通信作者),男,教授,博士生导师。 基金项目: 国家“973计划”资助项目(2012CB026002);国家“863计划”资助项目(2012AA052303)。

时间: 2013-10-22 网络出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131022.1113.003.html

O323

:A

:0253-987X(2014)01-0073-11