机翼极限环颤振的CFD/CSD全隐式紧耦合方法*

2020-08-06姚卫星姜志平周丹发

振动、测试与诊断 2020年3期

黄杰，姚卫星,2，姜志平，周丹发

(1. 南京航空航天大学机械结构力学及控制国家重点实验室南京，210016)(2. 南京航空航天大学飞行器先进设计技术国防重点学科实验室南京，210016)(3. 中国航空工业集团公司沈阳飞机设计研究所沈阳，110035)

引言

气动弹性是飞行器设计过程中必须考虑的问题，相比静气动弹性，颤振对飞行器带来的危害更严重，若设计不当，飞行动压达到颤振临界动压时，飞行器将发生机毁人亡的事故。早期的飞行器颤振分析基于经典颤振理论，典型的分析方法有v-g法和p-k法[1-2]。经典颤振理论为频域分析方法，非定常气动力的分析通常采用偶极子网格法[3-4]，该气动力模型为线性模型，只适用于亚声速和超声速问题，无法分析高度非线性的跨声速问题。当攻角较大时，流动分离也会导致气动力非线性，故偶极子网格法也无法分析大攻角问题。经典颤振理论的结构动力学通常采用模态叠加法分析，同样也为线性模型。因此，经典颤振理论仅能分析亚声速或超声速、小攻角、线性结构等条件下的颤振问题。

对于大型客机和高速无人机，巡航速度大多处于跨声速区域，其气动力是高度非线性的。而大展弦比无人机机翼柔度较大，在气动力作用下会产生较大的几何变形，机翼必然会产生明显的几何非线性效应，而几何非线性会造成翼面结构在非定常气动力作用下发生极限环振荡(limit cycle oscillation,简称LCO)[5-6]，同时发生颤振问题和LCO现象时可称为极限环颤振[7-9]。LCO颤振分析的关键是准确计算LCO幅值和频率。因此，经典颤振理论无法分析这类高度非线性的LCO颤振问题。随着计算流体力学(computation fluid dynamics,简称CFD)和计算结构动力学(computational structural dynamics,简称CSD)的发展，一些学者开始采用基于CFD/CSD的耦合方法进行LCO颤振的研究，非定常气动力通过有限体积法求解Euler方程或Navier-Stokes方程获得，而结构瞬态响应分析通常采用非线性有限元法，因此这是一种时域方法。目前LCO颤振的研究对象主要集中在壁板，而机翼的LCO颤振研究相对较少。Schairer等[10]通过试验方法研究了切尖三角翼的LCO颤振问题，获得了三角翼的模态及不同来流动压下的翼尖LCO幅值和LCO频率。文献[11-13]通过CFD/CSD耦合方法研究了切尖三角翼的LCO颤振，并将计算结果与试验结果进行了对比。该耦合方法具有松耦合的特点，耦合时间精度较低，具有时间滞后效应，已有学者证明了松耦合方法仅具有一阶耦合时间精度[14]。

针对松耦合方法的时间滞后效应及耦合时间精度较低的缺点，笔者发展了一种跨声速翼面极限环颤振的CFD/CSD紧耦合分析方法，其特点是在传统松耦合方法的基础上增加了内部伪迭代分析过程。采用Schairer等的试验模型进行了紧耦合和松耦合分析，并将计算结果和试验结果进行了对比，以验证紧耦合方法的耦合时间精度。

1 气动力和结构动力学控制方程

不考虑体积力和内热源情况下，直角坐标系下的流体动力学Navier-Stokes控制方程的积分形式为

(1)

其中：W为守恒向量；Fc为对流通量；Fv为黏性通量；∂V为控制体V的边界面；n为∂V的外法线单位向量。

(2)

(3)

(4)

(5)

其中：ρ为密度；p为压强；u，v和w分别为直角坐标系下的速度分量；qx，qy和qz分别为直角坐标系下的热流分量；e为单位质量气体的总能量，其表达式为p/((γ-1)ρ)+(u2+v2+w2)/2；γ为气体比热比；τij为黏性应力分量。

将式(1)按有限体积法进行空间离散可得

(6)

其中：Wi和Vi分别为控制体i的守恒向量和体积；NF为控制体边界面的数目；ΔSN为第N个边界面的面积。

由于对流通量Fc具有高度非线性特点，并且集中体现了流场的对流特征，笔者采用Roe格式[15]对其进行空间离散，为获得单调解，采用完全迎风的二阶MUSCL格式[16]离散分裂后的无黏通量,并采用minmod限制器使空间离散格式达到空间二阶精度，且非定常问题的求解采用双时间步长法[17]。

以上是气动力分析的控制方程和数值计算方法，而通过有限元法离散后结构动力学运动方程可表示为

(7)

本研究的机翼极限环非线性颤振仅考虑结构的几何非线性，不考虑其材料非线性，且采用Newmark法进行结构非线性瞬态动力学的求解，其属于隐式求解方法。将式(7)写为结构非线性动力学的平衡方程为

(8)

Fi体现了结构的几何非线性，需要将其线性化以方便求解

Fi,n+1=Fi,n+KT,nΔu

(9)

其中：Δu=un+1-un。

将式(9)带入式(8)，即可得到线性化的结构动力学的平衡方程

(10)

由于式(9)在一定程度上属于一阶线性化，故一阶线性化的式(10)求解精度较差，可引入子迭代Newmark方法来求解式(10)，其平衡方程为

(11)

基于子迭代Newmark方法的结构瞬态非线性动力学求解的详细过程可参考文献[14]。

2 全隐式紧耦合分析策略

本研究非定常气动力采用CFD方法求解，而结构非线性动力学采用基于有限元的CSD方法求解，且CFD和CSD均采用隐式时间推进，相比显示格式，隐式格式能采用更大的时间步长且具有更好的稳定性。

传统的机翼非线性颤振分析采用基于CFD/CSD的串行松耦合方法，其在相同的时间步进行气动力和节点位移的数据交换，并进行迭代推进求解。松耦合方法简单明了，实现起来较容易，但在任意n到n+1时间步内CFD(或CSD)求解过程中节点位移(或气动力)不变，即冻结边界条件，这样会造成时间滞后的效应，且随着分析的进行，累积误差会越来越大。文献[14]已经证明了即使CFD和CSD分析均达到了二阶时间精度，由于耦合边界上气动力和节点位移在交换时间上的滞后，传统松耦合方法只能达到一阶时间精度。这明显不满足机翼极限环颤振对耦合方法时间精度的要求。

基于传统松耦合方法时间精度较低的缺点，笔者采用全隐式紧耦合方法进行机翼极限环颤振的研究。其在任意一个时间步内CFD与CSD进行反复的内迭代，通常将其称之为伪迭代。当流体与结构满足精度要求，跳出伪迭代进行下一个时间步的求解，即进入物理迭代的分析。相比松耦合方法，紧耦合方法增加了伪迭代，并要求伪迭代收敛后才进行物理迭代的计算，这在一定程度上能消除松耦合方法在时间推进过程中累积的误差，即降低松耦合方法存在的时间滞后效应，具有更高的耦合时间分析精度。由于伪迭代完全收敛后再进入物理迭代，造成了计算量大大增加。因此在实际的紧耦合分析中通常需要设置一个最大伪迭代次数，即当伪迭代达到最大次数时即使其还未完全收敛，也跳出伪迭代进入物理迭代。

如图1所示，紧耦合方法的主要分析步骤如下：

图1 紧耦合方法分析流程Fig.1 Analysis process of tightly coupled method

1) 进行CFD和CSD建模，并进行CFD定常流场的计算，计算得到的定常气动力结果作为耦合分析的初始条件；

2) 进行ti时刻的非定常气动力的计算，并将壁面压力热Fi传递给CSD模型，进行CFD与CSD的伪迭代分析，直到伪迭代达到最大迭代次数后结束该时间步内的伪迭代分析；

3) 进行下一个时间步的分析，即进入物理迭代的计算；

4) 若物理迭代分析时间ti到达分析总时间ttotal，结束分析，否者返回步骤2，直至ti=ttotal为止。

以上耦合流程中涉及到气动力和节点位移在耦合面上的数据交换，而CSD的网格尺寸通常远大于CFD网格尺寸，故在耦合面上不能直接通过节点一一对应的方式进行以上耦合变量的传递，需要采用插值算法完成数据的传递。笔者采用三维薄板样条插值方法进行气动力和节点位移的插值计算。此外，结构变形会导致流场网格的重新调整，笔者采用基于无限插值方法的动网格技术实现流体网格的变形，其详细计算方法可参考文献[18]。

3 跨声速切尖三角翼极限环颤振分析

3.1 数值模型

笔者采用跨声速切尖三角翼试验模型[11]进行研究，通过试验可知，该模型的非线性因素来自于气动载荷作用下的结构大变形和跨声速非线性气动力。此剪切三角翼为等厚度薄板，其厚度为0.889mm，材料为钢，弹性模量为200GPa，泊松比为0.25，密度为7 850kg/m3，其平面形状和几何尺寸如图2(a)所示，且翼根为固支边界条件。来流马赫数在0.86和0.879之间，攻角为0°。为了与试验来流参数相匹配，来流动压的单位采用psi，其与Pa的变换关系为1psi=6 895Pa。图2为CFD和CSD数值模型，由于采用松耦合方法分析该模型极限环颤振的文献中未考虑流场的黏性效应，为了对比分析，本研究CFD模型仅求解Euler方程。此外CSD模型采用壳单元模拟，且考虑结构大变形带来的几何非线性效应。

图2 流体和结构模型网格Fig.2 Fluid and structural meshes

笔者进行了结构模态分析，结构振型的试验和分析结果如图3所示，机翼前3阶振型分别为一阶弯曲、一阶扭转和二阶弯曲，且分析结果与试验结果吻合得很好。表1列出了机翼前3阶固有频率的试验和分析结果[11]，固有频率的计算值与试验值最大相对误差为1.82%，这说明本研究的CSD模型能够反映切尖三角翼的基本动力学性能。

图3 机翼振型对比Fig.3 Comparison of modal shape of wing

表1 机翼固有频率的对比

3.2 机翼几何非线性的影响

为了研究结构几何非线性对颤振的影响，笔者采用紧耦合方法同时分析了来流动压q=2.58psi下考虑和不考虑几何非线性时的翼尖瞬态位移响应，其中耦合时间步长Δt取0.000 4s，每个时间步内的最大为迭代为10。图4为线性结构(不考虑结构几何非线性)和非线性结构下的结果对比图，可知考虑结构几何非线性时机翼在初始扰动下振动幅度逐渐增加，并且在最终呈等幅振荡状态。这是因为翼面为刚度渐硬的系统，即随着翼面变形的增加，翼面几何刚度也逐渐增加，机翼抵抗变形的能力得到了加强。线性刚度的机翼在分析过程中刚度保持不变，当受初始扰动时，机翼在初期呈现较小幅度的振动，随后进入剧烈的动态发散阶段，此阶段机翼的变形远大于考虑几何非线性情况下的结果。据此可知机翼几何非线性对时域分析结果产生了重要影响，必须考虑其影响才能获得正确的颤振分析结果。

图4 翼尖线性和非线性位移响应(动压q=2.58psi)Fig.4 Linear and nonlinear displacement responses of wing tip(q=2.58psi)

3.3 非线性临界颤振动压

笔者采用紧耦合方法进行了切尖三角翼的临界颤振动压分析，同时考虑了非线性跨声速气动力和结构几何非线性。图5为不同动压下机翼翼尖的法向位移时域分析结果，根据结果可知：动压q为2.42psi时，在外界扰动下翼尖的振动呈收敛的趋势；动压q为2.44psi时，翼尖的振动逐渐发散；动压q为2.43psi时，翼尖呈近似等幅振动的状态，故本研究的切尖三角翼模型非线性临界颤振动压为2.43psi。

图5 不同动压下翼尖位移响应Fig.5 Displacement responses of wing tip under different dynamic pressure

表2为p-k法[13]与本研究紧耦合方法计算得到的临界颤振动压与试验结果的对比情况，可知紧耦合方法获得的临界颤振动压值2.43psi与试验值2.40psi吻合良好，而p-k法计算得到的临界颤振动压为2.75psi，结果较差。这是因为p-k法为工程算法，其气动力和结构分析均为线性，无法考虑切尖三角翼的非线性跨声速气动力和结构几何非线性。表2中的数据同样也验证了紧耦合方法在求解非线性临界颤振动压中的高精度。

表2 非线性颤振动压比较

3.4 LCO分析结果和讨论

在切尖三角翼的颤振试验中观察到了LCO现象，其中的非线性来自于非线性跨声速气动力和结构几何非线性，且结构几何非线性是主要影响因素。笔者分析了不同来流动压下的翼尖法向位移的时域结果和相应的相位图，如图6所示。从图中结构可知，在初始扰动下翼尖振动在初始阶段呈现轻微的发散趋势，随后翼尖振动进入剧烈发散阶段，翼尖位

图6 不同动压下翼尖位移响应及相位图Fig.6 Displacement responses and phase diagrams of wing tip under different dynamic pressure

移迅速增加，并且最终作等幅振荡的状态，即表现出LCO现象。此外随着来流动压的增加，翼尖会更剧烈且更快地进入LCO状态。

分析获得了不同来流动压下的翼尖LCO幅值和频率，如图7所示，图7(a)中的纵坐标为LCD与机翼半展长b(其值为203.2mm)之比。为了对比分析，图中还给出了一些采用松耦合方法获得的结果。为了保证对比分析的合理性，本研究紧耦合方法和文献中松耦合方法的气动力分析控制方程均为Euler方程，且耦合时间步长均为0.000 4s。由图可知，紧耦合和松耦合方法获得的LCO幅值在较低来流动压下吻合较好，随着来流动压的增加，计算值与试验值误差越来越大，尤其在来流动压q=3.45psi时LCO幅值的试验值远大于所有计算值，但可以明显观察到采用紧耦合方法获得的LCO幅值更靠近试验值。随着来流动压的增加，LCO频率的试验值逐渐增加，紧耦合方法计算得到的LCO频率能清晰地反映这一规律，而松耦合方法却不能。此外，紧耦合方法获得的LCO频率也与试验值吻合得更好。

图7 不同动压下翼尖LCO幅值和频率Fig.7 LCO amplitude and frequency of wing tip under different dynamic pressure

以上分析说明了紧耦合方法比传统松耦合方法能更好地分析LCO幅值和频率，具有更高的分析精度。但是在较高来流动压下，紧耦合方法获得的LCO幅值及频率与试验值还存在较大的误差。主要原因有：①三角翼的扭转变形会造成局部大攻角，三角翼在大攻角下会形成前缘分离涡和前缘低压区，如图8所示，这难以通过Euler方程进行数值模拟，即存在气动力不精确的问题；②试验过程中机翼的振动容易引起翼根的松动[11]，降低约束刚度，造成振动幅值增加；③翼根部为高应力区域，有部分结构已进入塑性阶段，而本研究的分析不考虑材料非线性，这也将引起振幅增加。