基于广义Richardson外插方法的颤振模拟耦合时间精度研究

2021-05-04王运涛孟德虹

空气动力学学报 2021年2期

王昊，王运涛，孟德虹，王毅

(中国空气动力研究与发展中心计算空气动力研究所，绵阳 621000)

0 引言

20世纪80年代，美国、欧洲等发达国家和地区已对数值模拟可信度开展了大规模、有组织、有计划的研究工作[1]，在研究过程中提出了验证与确认的概念作为数值模拟可信度研究的重要内容。1998年美国航空航天学会（AIAA）发布了CFD验证与确认研究的指南[2]，提出了该领域的一些基本问题和基本概念。在此基础上，CFD验证与确认工作不断深入，国内外机构通过发布基准模型构建研究和比较平台，系统性地开展了一系列试验与计算工作。

可信度研究可以分为验证与确认两个过程。根据文献[2]的定义，验证过程是评估计算模型的求解精度，确认过程则是评估计算模型反映实际物理问题的精度，即本文对计算精度的研究属于验证过程。

目前CFD验证工作主要集中在定常问题，主要内容是空间离散以及网格尺度带来的影响，其中一项重要手段是网格收敛性研究，通过Richardson外插方法来获得网格无关解、分析计算空间精度和不确定度。但是对于非定常问题，时间离散带来的影响同样需要考虑。时间离散精度目前得到的关注较少，一般使用算法的理论分析精度来替代实际计算精度，没有进行相关的不确定度分析。对于气动弹性问题的时域求解，涉及到流场和结构的耦合计算，相较于单纯CFD问题更为复杂，其精度不仅取决于CFD以及CSD求解精度还和耦合算法有关，实际计算精度与理论分析精度难免有所差异。同时作为一种更为精细的气动弹性问题预测方法，时域模拟需要耗费大量计算资源，对其结果进行可信度分析，不仅能够得到可靠的计算结果，而且能够选取更为经济的时间离散步长，提高计算效率。

由于非定常问题的时间离散在数学上与空间离散并无本质上的差别，因此可以把时间离散看作一维空间离散问题，把空间离散问题的验证方法加以推广进行时间离散验证。本文参照网格收敛性分析方法，对颤振问题的时域模拟进行了时间步长收敛性分析，采用广义Richardson外插分析计算时间精度并获得时间步长无关解，采用网格收敛指数IGC（Grid Convergence index，GCI）估计时间离散不确定度。建立了颤振问题时域分析时间精度的验证过程，并检验了该方法的可行性。

1 Richardson外插方法及验证方法

1.1 Richardson外插方法

Richardson外插方法又称“h2外插法”、“迭代外插法”，于1910年[3]由Richardson首次使用，并于1927年[4]加以完善。该方法使用离散步长h将离散解f与精确解fexact的关系表示为级数的形式：

其中g1、g2等参数与离散过程无关。

对于二阶精度方法g1= 0，此时只需离散步长分别为h1和h2时的离散解f1、f2即可求得离散无关解fexact。

根据离散解f1、f2计算方法的不同，外插fexact为三阶或四阶精度。

1.2 广义Richardson外插方法

采用Roache[5]的方法Richardson外插可推广至p阶方法，可称之为广义 Richardson外插方法。此时离散解与精确解之间的关系可以表示为：

对于非定常问题时间离散，选取时间步长h3＞h2＞h1及其对应计算结果f1、f2、f3，代入公式（3）可得：

由式（4）可得精度p：

其中：

为保证分析结果准确性，时间步长选取时需要保证ε32、ε21不能太过接近于0，经验估计r需要大于1.3[6]。

外插得到“精确解”：

在使用Richardson外插方法进行分析的过程中，需要注意所选结果离散步长不能够太大，此时计算过程没有充分收敛，分析过程中忽略的高阶项会对分析结果产生影响，使分析结果具有很强的不确定性。

1.3 验证方法

采用上述结果可继续进行误差及不确定度估计。

近似相对误差：

外插相对误差：

网格收敛指数[7]：

本文选取安全系数Fs=1.25。

2 气动弹性分析方法

本文研究工作基于TRIP软件气动弹性模块[8-9]开展。该模块主要包括流场计算、结构计算、耦合界面插值和动网格四个主要功能部分，依照耦合计算方法的流程依次调用。

2.1 计算方法

流场求解器采用课题组自主开发的亚跨超声速CFD软件平台TRIP。经过课题组多年的发展，TRIP已经成为一个非常成熟的数值模拟软件平台，大量验证工作[10-11]表明TRIP软件具有较高的计算精度和效率。

结构求解模块采用基于模态坐标的气动弹性运动方程，可使用龙格库塔法、线性多步法、双时间步方法等多种方法求解。为避免CFD/CSD耦合求解中的质量不相似问题，本文采用变刚度方法[12-13]确定颤振边界。

动网格模块集成了目前工程应用中常用的TFI、Delaunay以及RBF插值方法，以及基于这些基础方法开发的一些复合型动态网格变形方法，如RBF_TFI、RBF_Delaunay等。界面插值模块集成了无限平板样条IPS、薄板样条TPS和径向基函数插值RBF三种插值方式。本文二维算例Isogai Wing模型为二元翼型，满足刚体运动假设，动网格方法采用刚性旋转法。为简化计算，该算例中结构求解部分与流场求解采用相同网格，无需使用界面插值技术。

2.2 耦合方法

目前实际应用中的耦合计算方法主要分为松耦合和紧耦合两类。松耦合方法通过交替求解流场和结构方程进行时间推进，不在单场求解过程中进行流场与结构的数据交换。该方法能够充分利用现有的流场和结构求解器，实现过程简单，应用广泛。但是松耦合方法为了计算简便往往造成流场和结构的时间不同步，耦合精度只有一阶。为提高松耦合方法计算精度，部分学者采用预估校正思想对松耦合方法进行改进，提出了具有二阶精度的改进的松耦合方法[14-16]。紧耦合方法则是将流场和结构方程均构造为子迭代形式[17-18]，在每个物理时间步内，对流场和结构进行多次数据交换，来消除流场与结构信息交换不同步的问题，时间精度能够达到二阶。

为充分验证所建立精度分析方法的可靠性，本文分别选取松耦合方法、改进的松耦合方法和紧耦合方法的结果进行分析。松耦合方法采用冻结气动力的龙格库塔方法[14]，即在松耦合流程中结构求解使用简化的龙格库塔方法。作为对照，使用相同耦合流程，使用双时间步方法作为结构求解方法，构造了使用双时间步方法的松耦合方法，两种松耦合方法理论时间精度均为一阶。改进的松耦合方法采用二阶杂交的线性多步法[15]，理论精度为二阶。本文所使用紧耦合方法根据文献[17]构造，理论精度为二阶。同时作为对照，本文将龙格库塔方法作为结构求解嵌入流场求解子迭代过程中，自行构造了一种一阶精度的紧耦合方法。本文选取上述多种精度的共五种耦合方法进行计算，并对计算结果进行精度分析。五种耦合方法的对比如表1所示。

表1 耦合方法对比Table 1 Comparison of coupling methods

3 Isogai Wing模型颤振计算

3.1 计算模型及网格

Isogai Wing[19-21]是由Isogai提出的后掠机翼颤振问题的典型截面，是检验气弹不稳定性预测方法的二维标准算例。翼型截面采用NACA010翼型，属于二元机翼，翼型结构如图1所示。图中来流速度为U∞，机翼半弦长b= 0.5 m，在弹性轴（对应于刚心）处固定一个垂直方向的线弹簧以及一个扭转弹簧，弹簧刚度分别为Kn和Kα，弹簧阻尼分别为Ch和Cα。选取机翼后缘一侧为正方向，弹性轴在距翼弦中点（C.L.）ab处，其中a= -2.0，表明弹性轴位置在翼弦中点前方，重心在距弹性轴bxα处，其中xα= 1.8，重心位置在翼弦中点后方。机翼具有绕弹性轴俯仰运动和上下平移的浮沉运动两个自由度，俯仰运动由转角α表示，前缘向上为正，浮沉运动由弹性轴的上下位移ξ表示，向下为正。机翼无量纲回转半径rα2= 3.48，浮沉运动固有频率ωh= 100 rad/s，俯仰运动固有频率ωα= 100 rad/s，质量比μ= 60。计算网格为ECERTA Project网站[22]提供的适用于Euler方程计算的结构网格。

图1 Isogai Wing结构模型[21]Fig. 1 Structural model of Isogai Wing[21]

3.2 颤振计算结果

选取Ma= 0.6，颤振临界状态下无量纲来流速度vf= 1.710，不同时间步长下五种耦合方法的结构响应对比如图2所示，横轴为时间t，纵轴为结构一阶广义位移x1。

由结果对比可知，时间离散步长对计算结果具有显著影响。由图2（a）可知，时间步长较大时，五种耦合方法的结构响应存在较大差异，表明此时不同耦合方法的颤振边界存在较大差异；随着时间步长的减小（图2（b）），结果趋向一致，颤振边界趋向收敛；时间步长进一步减小（图2（c））五种方法所得结果基本重合，此时五种方法所得颤振边界基本一致，即可推断随着时间步长减小各耦合方法计算颤振边界均趋向收敛，时间步长足够小时五种耦合方法均能够得到正确结果。

图2 不同时间步长下机翼时域结构响应对比Fig. 2 Comparison of time domain structural response with different time step sizes

通过图2可以对五种耦合方法进行初步精度分析：两种二阶精度耦合方法计算所得结构响应曲线基本重合，并且不受时间步长的影响，即两种二阶精度方法在计算中表现出的精度基本一致；三种一阶精度耦合方法结构响应存在一定差异，即实际计算精度存在差异，结构求解采用双时间步方法的松耦合方法精度略小于冻结气动力的龙格库塔方法、一阶精度的紧耦合方法计算精度小于两种二阶精度方法。由于随时间步长减小结构响应曲线趋向收敛的方向不同，这种定性分析方法并不能够直接对比一阶精度的松耦合方法与二阶精度方法的计算精度，这也说明了本文所建立的精度分析方法的必要性。

4 Isogai Wing模型数值精度分析

为系统地研究时间离散对计算结果的影响，在Ma= 0.6条件下选取一系列时间离散步长计算模型的颤振边界。颤振边界的确定采用对数减幅法，选取颤振临界速度和颤振临界频率作为目标变量进行分析。

4.1 时间收敛性

为确定计算结果的时间收敛性，如图3所示，对比五种方法在不同时间步长下的颤振临界速度和颤振频率。五种方法颤振临界速度和颤振频率均随时间步长的减小而单调变化，最终收敛于同一点。即五种方法均具有良好的时间收敛性，可以使用Richardson外插方法进行分析。

图3 不同时间步长下Isogai Wing模型颤振边界计算结果Fig. 3 Flutter boundary simulation results of Isogai Wing with different time step sizes

4.2 广义Richardson外插方法分析

首先选取合适时间步长及其计算结果。根据1.2节分析，所选取时间步长不应过大，同时应使r大于1.3，另考虑对数减幅法判断颤振临界状态带来的误差，同样需要避免选取过小的时间步长，最终选取时间步长0.05 、0.1、0.2 ms，即每周期约800、400、200个时间步。

使用广义Richardson外插方法分析所选计算结果实际精度，并计算不确定度。基于颤振临界速度分析可得表2。结果表明，分析所得精度p与理论精度基本相符。两种松耦合方法和两种二阶精度耦合方法实际精度均略高于理论精度，一阶精度紧耦合方法精度略低于理论精度。冻结气动力的龙格库塔方法精度略高于采用双时间步方法的松耦合方法，改进的松耦合方法和传统紧耦合方法精度均高于一阶精度紧耦合方法，改进的松耦合方法和传统紧耦合方法精度及误差基本相同，与3.2节的定性分析结果一致，说明该方法得到了可信的分析结果。五种方法插值所得时间离散无关解之间最大误差为0.15%，在精度允许范围内，五种方法分析得到了一致的“精确解”。三种误差和的对比则说明相同时间步长下二阶精度方法计算误差更小，具有明显的精度优势。

表2 基于颤振临界速度的时间精度分析结果Table 2 Time accuracy analysis results based on the flutter critical velocity

如表3所示，基于颤振临界频率与颤振临界速度的分析结果基本一致。但是一阶精度紧耦合方法、改进的松耦合方法和传统紧耦合方法精度与颤振临界速度分析结果均有一定差异，表明使用广义Richardson外插方法进行精度分析时，所选取目标量对分析结果有一定影响。

若计算要求颤振临界速度误差小于1%，各耦合方法能够取的最大时间步长如表4所示。二阶精度耦合方法能够大幅提高计算效率，体现了进行精度分析的重要意义。

表3 基于颤振临界频率的时间精度分析结果Table 3 Time accuracy analysis results based on the flutter critical frequency

表4 外插相对误差小于1%时的最大时间步长Table 4 Maximum time step for ＜1%

Coupling method 1st order loosely coupled by RK 2nd order tightly coupled Δtmax /ms 0.16 0.08 0.16 0.5 0.5 Step per period 240 480 240 80 80 1st order loosely coupled by dual 1st order tightly coupled 2nd order loosely coupled

4.3 Richardson外插方法分析结果验证

基于广义Richardson外插方法进行时间精度分析需要选取三个计算时间步长的计算结果，因此需要考虑特定时间步长是否对分析结果产生影响。为避免所选时间步长结果带来的偶然性，对4.1节使用的计算结果进行进一步处理，以确定该方法进行时间精度分析的可行性。

图4 Isogai Wing模型颤振边界计算结果时间精度曲线Fig. 4 Time domain accuracy results for Isogai Wing with flutter boundaries

如1.2节所述，在时间步长较大时，分析结果具有很强的不确定性，不适用于采用广义Richardson外插方法进行精度分析，因此图4中时间步长大于0.6 ms（每周期约60个时间步）时曲线线性度较差。在时间步长较小时，计算误差相对较小，使用对数减幅法判断颤振边界引入的误差对精度分析带来较大干扰，因此图4中时间步长小于0.04 ms （每周期约1 000个时间步）时，曲线线性度相对较差。若去除时间步长过大和过小的部分，曲线则具有良好的线性度。为更直观地说明剩余点拟合曲线良好的线性度，对时间步长在0.04～0.6 ms之间的点进行线性回归分析。线性回归所拟合直线如图5所示，直线斜率k和可决系数R2在表5和表6中给出。

拟合直线斜率k即耦合计算精度，可决系数R2则反映了拟合直线与选取数据点的接近程度。线性回归分析所得计算精度与4.2节分析结果基本一致，R2非常接近于1，数据点与拟合直线相当接近，即所选取数据点具有良好的线性度。线性回归结果表明在恰当的时间步长区间内，使用广义Richardson外插方法分析所得结果与具体时间步长选取无关，也进一步验证了本文所提出时间精度分析方法的可行性。

图5 Isogai Wing模型颤振边界计算结果线性回归拟合直线Fig. 5 Linear regression fitting for Isogai Wing with flutter boundaries

表5 基于颤振临界速度的线性回归分析结果Table 5 Linear regression analysis results based on the flutter critical velocity

表6 基于颤振临界频率的线性回归分析结果Table 6 Linear regression analysis results based on the flutter critical frequency

5 AGARD 445.6 Wing颤振计算及精度分析

5.1 计算模型

AGARD 445.6 Wing[23-24]被广泛应用于颤振计算程序的验证工作。机翼材质为均匀的桃花心木薄板，1/4弦长处后掠角为45°，机翼截面为NACA65A004对称翼型。本文计算网格单元数为321万，结构计算选取前四阶模态，模态频率为试验所得频率，计算马赫数0.499。

5.2 时间精度分析

选取颤振临界速度和颤振临界频率为分析变量，分析三维颤振标模AGARD 445.6的计算精度。由于计算量的限制，相较于二维算例，所选计算时间步长及耦合方法较少。

如图6所示，该算例具有良好的时间收敛性，可以使用Richardson外插方法进行精度分析。根据1.2节及4.2节所述时间步长选取原则，考虑一阶精度与二阶精度耦合方法收敛性的差异，一阶精度松耦合方法选取时间步长为0.05、0.1、0.2 ms （每周期约900、450、225个时间步），改进的松耦合方法选取时间步长为0.2、0.4、0.8 ms （每周期约225、113、57个时间步）。

采用上述时间步长结果进行精度分析，结果如表7、表8所示。改进的松耦合方法基于颤振临界速度的精度与理论精度差异稍大，其余分析精度与理论精度基本一致。由于三维问题更为复杂，相较于二维算例分析结果，两种耦合方法实际计算精度均有所下降。对于两个目标变量，两种耦合方法得到的时间离散无关解误差分别为0.07%、0.18%，均得到了一致的“精确解”。由以上分析结果可知，该精度分析方法在该三维算例中得到了合理的分析结果。

为进一步研究三维算例中时间步长选取对分析结果的影响。采用4.3节所述方法，得到如图7所示精度曲线。由于计算结果已经收敛，计算误差相对较小，改进的松耦合方法精度曲线在小时间步长处线性度较差。由于精度较低，松耦合方法精度曲线在较大时间步长处线性度较差。两种耦合方法结果分别去除较小时间步长或较大时间步长的点进行线性回归分析，分析结果如表9、表10所示，拟合直线如图7虚线所示。

图6 不同时间步长下AGARD 445.6 Wing模型颤振边界计算结果Fig. 6 Flutter boundary simulation results of AGARD 445.6 Wing with different time step sizes

表7 基于颤振临界速度的时间精度分析结果Table 7 Time accuracy analysis results based on the flutter critical velocity

表8 基于颤振临界频率的时间精度分析结果Table 8 Time accuracy analysis results based on the flutter critical frequency

图7 AGARD 445.6 Wing模型颤振边界计算结果时间精度曲线Fig. 7 Time domain accuracy results for AGARD 445.6 Wing with flutter boundaries

表9 基于颤振临界速度的线性回归分析结果Table 9 Linear regression analysis results based on the fluttercritical velocity

表10 基于颤振临界频率的线性回归分析结果Table 10 Linear regression analysis results based on the flutter critical frequency

线性回归分析所得计算精度k与表7、表8分析结果基本一致；R2均大于0.99，即数据点与拟合直线相当接近，具有良好的线性度。线性回归结果表明，在恰当的时间步长区间内，使用广义Richardson外插方法分析所得结果与具体时间步长选取无关，进一步验证了本文所提出时间精度分析方法针对三维颤振问题的可行性。