NNW-FSI软件静气动弹性耦合加速策略设计与实现

2021-10-20孙岩王昊江盟岳皓孟德虹

航空学报 2021年9期

孙岩，王昊，江盟，岳皓，孟德虹

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

航空工业飞行器设计一般需要经历概念、初始和详细3个关键阶段[1]。在概念设计阶段，设计师利用半经验方法或低可信度的气动分析工具在众多设计方案中做出快速选择[2]。在初始设计阶段，气动/结构布局方案将冻结，并确定飞行器目标设计载荷。目标设计载荷是飞行器全机或部件设计必须满足的极限载荷，所有细节设计必须在此载荷框架下进行[3]。在详细设计阶段，将根据目标设计载荷，确定飞行器的全部结构细节。

飞行器设计有朝着复合材料等轻质结构发展的趋势[4]，以减轻全机结构重量，提高燃油效率和航程，如波音787和空客A350XWB飞机的全机复合材料比例已经超过50%。飞行器变得更加柔性，在全部设计阶段考虑机翼气动/结构之间的耦合作用能够更好地改善飞行器的性能、降低设计制造成本[5]。

近十年，随着计算机能力和并行计算技术的快速发展，高可信度的计算流体力学(CFD)/计算结构力学(CSD)耦合方法已经逐步应用于飞行器的初始和详细设计阶段，而在概念设计中，因为存在大量的计算状态，仍然采用快速的半经验或低可信度分析工具。但Rizzi认为飞行器80%的全生命周期成本来源于概念设计阶段的决定，且该阶段由于设计工具自身误差导致的错误在后续阶段难以修复，引起的代价和风险极高[6]。因此，发展静气动弹性耦合加速收敛技术对于促进高保真度CFD方法在航空飞行器设计过程中的应用具有潜在的价值。

当前，基于CFD的静气动弹性耦合模拟方法主要采用气动与结构直接耦合的方式[7-11]，CFD计算气动载荷，通过耦合界面数据传递将载荷传递给结构，结构计算变形，再通过耦合界面数据传递将变形传递给流场，如此反复，直至气动和结构均达到收敛。当前的耦合策略在大多数常规飞行器变形问题中可以取得良好的收敛效果，但在一些大柔性机翼结构的静气动弹性耦合模拟中，作者所在课题组发现当前耦合策略可能引起机翼反复振荡，降低收敛效率，甚至极端情况下会引起结构变形发散。

针对不同类型算例表现出的收敛特性差异，在国家数值风洞(National Numerical Windtunnel, NNW)工程的资助下，依托NNW-FSI软件平台，研究静气动弹性耦合加速策略的设计与实现，探索耦合加速的理论机制，为工业应用中参数选择提供理论依据。

本文第1节对资助整个研究的国家数值风洞工程做简要的介绍，第2节对NNW-FSI软件平台的概况、组成、功能、性能、开发进展和规划进行介绍，第3节介绍静气动弹性耦合加速策略的设计与实现，第4节开展两种外形算例的静气动弹性耦合模拟，测试耦合参数对收敛效率的影响，第5节根据静气动弹性模拟结果分析耦合加速的理论机制，第6节为结论。

1 国家数值风洞(NNW)工程

国家数值风洞工程是中国空气动力研究与发展中心(CARDC)牵头，国内外高校、企业和科研院所协同参与建设，致力于发展以CFD为核心的数值模拟软件体系，构建CFD产业应用生态，推动国产制造业转型升级的大型自主CFD软件研发项目[12-13]。

NNW工程于2018年10月批复启动，将通过两期进行建设。一期包括基础科学问题研究、CFD软件及集成框架、多学科耦合应用软件、验证与确认、高性能计算机5大系统，将在2022年底基本建成完全自主知识产权的CFD软件群；二期将持续优化软件体系结构，建成满足更高品质工业设计需求的新型国家数值风洞。

2 NNW-FSI软件

2.1 软件概况

NNW-FSI软件是NNW工程多学科耦合应用软件系统(简称多学科系统)下的流固耦合模拟软件，主要针对航空航天飞行器快速精细发展的迫切需求，兼顾船舶、桥梁、高层建筑以及风力机等产业领域的设计需要，依托NNW工程CFD软件及集成框架的建设成果，通过流动/结构多学科耦合模拟关键技术的突破以及耦合数据传递、网格变形、耦合策略等核心功能模块的开发，形成面向不同工程问题提供快速高保真度流固耦合数值解决方案的能力。

2.2 软件组成

图1给出了NNW-FSI软件的模块组成，主要包括流场解算、结构解算、耦合策略、网格变形和数据插值5个核心模块。其中流场解算将依托NNW的CFD软件开发成果，结构解算采用模态、柔度和有限元求解3种方式，耦合策略有面向静态问题的松耦合和面向动态问题的紧耦合[14]，网格变形有RBF[15-16]、RBF_TFI[17]、RBF_Delaunay[18]和RBF_HBGM[19]等多种策略，数据插值采用平板样条(IPS)、薄板样条(TPS)和径向基函数(RBF)插值等[20]。

图1 NNW-FSI软件模块组成Fig.1 Module components of NNW-FSI software

2.3 软件功能及性能

NNW-FSI是多学科系统的核心功能软件之一，开发完成后将能够开展静态、动态流固耦合问题的数值模拟，具备动态流固耦合问题的快速分析能力。

针对航空航天飞行器研制需求，静态流固耦合功能主要解决飞行器变形后的气动载荷重新分配、机翼弹性扭转发散边界预测、舵面效率与舵面反效等问题。动态流固耦合功能主要解决颤振边界预测、非线性极限环响应、突风载荷下结构响应、非线性抖振等问题。动态流固耦合快速分析主要基于非定常流动建模，实现颤振边界和突风载荷响应的快速预测。

性能方面，NNW-FSI软件预期可以实现亿量级流场网格、千万量级结构自由度流固耦合模拟，并行计算规模达万核量级。效率方面，单次静态流固耦合模拟时间与定常流场计算时间量级相当，能够在数小时内完成，基于高性能计算机集群系统的容量计算能力，大量状态的常规静气动弹性计算分析可以在数天内完成。时域动态流固耦合模拟的时间消耗在工程型号问题可接受范围，多个状态的颤振分析可以在一周内完成。精度方面，典型构型典型状态下，NNW-FSI软件预期气动和结构载荷预测偏差不高于5%，安全边界预测偏差不高于10%。

2.4 软件开发进展及预期规划

2018年NNW工程启动后，NNW-FSI软件在早期in-house代码TRIP-StAE开发工作的基础上[10-11,21-23]，持续进行代码重构优化和功能扩展完善，目前已经完成了全部核心功能模块初始版本的开发，部分功能模块已经进行了多轮迭代完善，效率精度均优于同等水平的商业软件，基本具备了开展超大规模网格静/动态气动弹性数值模拟的能力，但程序设置仍然较为繁琐，学习成本较高，软件易用性水平与商业软件还有非常大的差距。

NNW-FSI软件后续将进一步完善功能、改进精度、提升效率，开展软件预测精度验证与确认研究，提高软件核心质量，同时在输入参数配置、计算文件准备、输出数据定制化方面开展流程优化，提升用户的软件使用体验，打造提供高品质流固耦合数值解决方案的仿真软件平台。

NNW-FSI软件在2020年底形成了面向用户使用的串并行执行版本，能够开展高精度的静态流固耦合模拟和部分动态流固耦合模拟，并在单位内部试用和推广；2021年底将形成带完整界面的串并行版本，能够开展全部的静/动态流固耦合数值模拟，具备设计和评估能力，向主机设计单位和科研院所进行推广应用；2022年6月，将形成最终完整版本，向全社会推广使用。

3 静气动弹性耦合加速策略

3.1 设计思想

静气动弹性耦合数值模拟中，CFD首先根据流动参数计算出给定气动外形的流场，并积分得到作用在飞行器表面的气动载荷，通过耦合界面数据传递转换为施加在结构加载点的输入载荷，计算结构力学(CSM)根据输入载荷和边界条件计算出飞行器结构的静变形，再通过耦合界面数据传递转换为流场物面网格点的位移，利用网格变形技术生成新的外形下的流场计算网格，CFD根据新的流场网格计算新外形下的流场，并重复上述流程，直至流动和变形均趋向收敛。

从上述静气动弹性耦合流程可以知道CFD与CSM之间需要经过多次反复迭代才能收敛至结构静平衡位置。而从静气动弹性耦合的物理机制可以得出，在气动外形和结构形式给定的前提下，结构的静平衡位置由来流和初始条件唯一确定，与CFD和CSM之间迭代的中间过程没有关系。

因此，可以根据静变形过程中飞行器外形变化规律，采用不同的耦合迭代策略，加速耦合迭代收敛的过程。根据静气动弹性耦合的物理机制，建议耦合加速迭代策略设计的两个原则为：① 初次结构变形尽可能地接近结构静平衡位置；② 耦合迭代过程尽可能避免结构变形的反复振荡。其中，原则 ① 是为了获得接近收敛解的初值，降低初始误差，原则 ② 是为了加快初始解向平衡位置的收敛。

3.2 加速策略实现

根据两个设计原则，耦合加速策略通过采用增量叠加方式实现:

un+1=un+λΔun+1

(1)

(2)

将式(2)代入式(1)，可以得到

(3)

需要说明的是，采用的耦合迭代加速策略没有将初次结构变形u1和结构变形增量Δun+1的调整通过不同参数来实现，而是统一采用了参数λ来调整。主要基于以下两个原因;一是从物理机制的角度考虑，由于最终结构静平衡位置是未知的，因此缺少调整初次结构变形u1值的参考，工程上通常根据不同飞行器结构的静变形规律和使用经验来确定一个取值范围;二是从软件开发角度考虑，采用单个参数可以统一程序结构设计，降低软件开发的复杂度，并减少用户自由参数的选择，降低软件使用的难度。

从后续的算例测试也可以进一步看出，初次结构变形u1的大小对收敛效率的影响要远低于结构变形增量Δun+1的影响，但对整个静气动弹性耦合模拟的稳定性有较大影响。

4 静气弹耦合算例模拟

采用超大展弦比无人机和CHN-T1模型两类外形对NNW-FSI软件的静气动弹性耦合加速策略的参数影响和加速效果进行测试和评估，其中超大展弦比无人机模型使用超长平直机翼加细长机身的翼/身/平尾/立尾组合体构型，CHN-T1模型使用大展弦比后掠机翼加圆形机身的翼身组合体构型。

4.1 超大展弦比无人机耦合模拟

图2给出了超大展弦比无人机模型外形，该模型平面布局参考国外高空长航时太阳能无人机形状进行设计，采用翼/身/平尾/立尾组合体构型，主要用于NNW-FSI软件的算法测试和性能演示。

图2 超大展弦比无人机外形Fig.2 Ultra-high respect-ratio unmanned aerial vehicle configuration

该无人机模型采用超长平直机翼，翼剖面使用NACA0012翼型，翼根位置弦长为1 m，平直段长度为12 m，为了降低气动载荷、减小翼根弯矩，在靠近翼梢附近做了后掠20°的修形处理，机翼展长为27 m。机身采用流线型细长水滴形状，长度为6.312 m。平尾和立尾翼剖面也均采用NACA0012翼型，平尾展长为1.5 m，立尾高度为0.8 m。

无人机模型的计算条件设定为来流速度V=30 m/s，流动攻角α=3.0°，基于平均气动弦长的雷诺数Re为3.0×106，来流动压q=1 000 Pa。流场离散采用多块结构网格，共有84个网格块，5 653 272 个网格点，5 280 896个六面体网格单元。机翼结构有限元模型采用四边形壳单元进行离散，共有495个网格点，432个壳单元，壳单元厚度δ统一为0.05 m，结构材料采用金属铝，弹性模量E为70 GPa，泊松比为0.3。结构有限元模型通过翼根固支的方式进行约束。

静气动弹性耦合模拟中，数据传递采用平板样条插值IPS，网格更新采用RBF_TFI，流场每计算100步进行一次耦合数据传递和位移更新。

图3给出了不同松弛因子λ下超大展弦比无人机模型静气动弹性耦合过程中直机翼最大弯曲位移随耦合迭代步的收敛曲线，图中横轴表示耦合迭代步数，纵轴dymax表示机翼的最大弯曲位移。可以看出，不同松弛因子λ的耦合迭代过程均收敛于相同的值，通过实例验证了静气动弹性耦合的最终收敛解与初次结构变形大小和迭代过程无关。

对于图3超大展弦比无人机算例，从收敛过程来看，λ≤0.30时，位移的收敛过程是单调的，随着λ的增加，收敛显著加快；λ≥0.50时，位移的收敛过程开始出现振荡，且λ越大，振荡越明显，例如λ=1.50时存在多次反复振荡，但不同λ的收敛效率相差不大，均在15个耦合迭代步，即CFD流场计算约1 500步时，位移达到收敛。相比单调收敛过程(λ≤0.30)，带振荡的收敛过程(λ≥0.50)反而要更加快速，说明较小的λ值时，位移收敛缓慢，需要多次迭代才能逼近收敛解，而较大的λ值时，残差被快速衰减，计算结果迅速逼近收敛。

从无人机初次结构变形dymax(迭代步数为1)的大小来看，λ=1.0时，初次结构变形最接近收敛解，如图3所示，但逼近最终变形的初次结构变形并没有明显加速耦合迭代的收敛过程，对整体收敛效率的影响很小。

图3 最大弯曲位移收敛曲线(直机翼)Fig.3 Convergence curves of maximum bending displacement (straight wing)

图4给出了变形前后超大展弦比无人机机翼的升力载荷分布。可以看出，变形前后机翼升力载荷分布变化很小，说明无人机机翼变形以弯曲为主，不存在明显的扭转变形，机翼变形对气动载荷的影响很小。

图4 超大展弦比无人机机翼载荷分布Fig.4 Wing load distribution of ultra-high respect-ratio unmanned aerial vehicle

4.2 CHN-T1翼身组合体耦合模拟

CHN-T1模型是中国空气动力研究与发展中心面向国产新一代窄体干线客机而设计研制的大展弦比飞机标模[24]，采用翼/身/挂架/吊舱/垂尾/平尾组合体布局，旨在通过不同类别的高质量风洞试验数据，探索窄体客机气动布局特性和评估国产CFD软件的可信度水平[25]。

为了进一步研究大展弦比柔性机翼气动弹性特性和评估国产气动弹性计算软件可信度，基于CHN-T1模型翼身组合体构型设计制造了柔性风洞试验缩比颤振模型，模型缩放比例为1∶10，缩放后的机翼展长b为2 986.9 mm。图5给出了CHN-T1模型翼身组合体构型，有关CHN-T1模型的几何参数和详细信息可以参考文献[24]。

图5 CHN-T1模型翼身组合体构型Fig.5 Wing-body configuration of CHN-T1 model

颤振风洞试验在CARDC的2.4 m暂冲式跨声速风洞开展，为了匹配风洞的运行速压，颤振模型设计偏软，在气动载荷作用下存在显著的变形。需要通过静气动弹性耦合获取模型变形后的气动载荷，评估模型的结构强度安全。

CHN-T1模型的计算条件设定为来流Ma=0.70，流动攻角α=1.5°，基于平均气动弦长的雷诺数Re为8.565×106，来流动压q=30 kPa。流场离散采用多块结构网格，共有69个网格块，8 669 469 个网格点，8 166 208个六面体网格单元。机翼结构有限元模型采用混合单元进行离散，包括三角形和四边形的壳单元、四面体和六面体的体单元，共有18 056个节点，13 077个单元，梁、肋板和蒙皮结构材料均为复合材料。

静气动弹性耦合模拟中，数据传递采用薄板样条插值TPS，网格更新采用RBF_TFI，流场每计算200步进行一次耦合数据传递和位移更新。

图6给出了不同松弛因子λ下CHN-T1模型静气动弹性耦合过程中机翼最大弯曲位移随耦合迭代步的收敛曲线。不同λ值下CHN-T1模型的收敛过程再次证明收敛解与初次结构变形大小和迭代过程无关。

对于图6的CHN-T1模型算例，从收敛过程来看，仅λ=0.50时收敛过程是单调的，且收敛速度最快，在8个耦合迭代步，即CFD流场计算约1 600步时，位移达到收敛；λ≤0.30时，收敛过程缓慢，随着λ的增加，收敛速度加快；λ=0.75时，位移出现反复的振荡，且收敛十分缓慢，50个耦合迭代步后，位移仍存在小幅振荡，没有达到完全收敛；λ=1.0时，结构位移大幅振荡，在第3步位移即过大引起网格破裂，导致耦合计算失败。

从CHN-T1模型初次结构变形dymax(迭代步数为1)的大小来看，λ=0.2时，初次结构变形最接近收敛解，如图6所示，但逼近最终变形的初次结构变形也没有明显加速耦合迭代的收敛过程，对整体收敛效率的影响很小。

图6 最大弯曲位移收敛曲线(后掠机翼)Fig.6 Convergence curves of maximum bending displacement (swept wing)

图7给出了变形前后CHN-T1模型的升力载荷分布，其中图7(a)为未变形前的升力载荷分布，图7(b)为耦合迭代1步的升力载荷分布，

图7 CHN-T1模型机翼载荷分布Fig.7 Wing load distribution of CHN-T1 model

图7(c)为耦合迭代2步的升力载荷分布，图7中耦合迭代松弛因子λ的值取为0.75。

可以看出，变形前后机翼升力载荷分布发生显著变化，说明CHN-T1模型机翼不仅发生了弯曲，还产生了明显的扭转。这是由于CHN-T1模型机翼采用了后掠布局，在正弯曲变形的情况下，流向翼剖面后缘位移大于前缘位移，诱导翼剖面产生负的扭转角，导致机翼升力减小，如图7(b) 所示，在翼梢附近产生了正变形，却诱导产生了负的升力。随着机翼载荷的减小，机翼结构变形反转，产生了负的变形，又诱导产生正的翼剖面扭转角，机翼升力增加，如图7(c)所示。

5 耦合加速机制分析

第4节对两种类型大展弦比机翼外形开展了静气动弹性耦合模拟，两种机翼外形表现出了不同的收敛特性。本节将基于上述模拟结果，通过进一步的理论分析，弄清不同收敛过程的产生机制，以及松弛因子λ对收敛过程的加速效果，为工程应用中λ的选取提供依据。

5.1 结构变形收敛性

定义结构变形的最终收敛解为uconv，则第n步结构变形计算的误差εn可以定义为

εn=un-uconv

(4)

定义变形误差εn的增长率γn为

(5)

式中：m为结构变形点的数量。将式(3)代入到式(5)，可以得到

(6)

根据迭代收敛理论，当误差增长率γ的绝对值恒小于1.0时，变形误差会随着迭代步数逐渐减小，结构变形最终趋向收敛，当γ值恒等于1.0时，变形误差将保持不变，当γ值恒等于-1.0时，变形误差将持续等幅振荡，当γ的绝对值恒大于1.0时，变形误差将逐渐增大，结构变形将趋向发散。

(7)

式中：f(εn)表示结构变形相比平衡位置变形的偏离诱导的气动力变化引起的结构变形计算结果的变化量，且f(0)=0。将式(7)代入式(6)，可以得到

(8)

式中：fi表示第i个结构点变形的变化量。

5.2 收敛性讨论

式(8)中，γi的值决定了结构变形的收敛特性，当耦合迭代松弛因子λ给定时，γi值仅通过fi(ε)/εi的值决定。fi(ε)/εi的值受到来流参数、气动特性、结构特性等多种因素的影响，可以大致分为3类情况。

1) |fi(ε)/εi|≤δ，δ为某个小的正数。说明结构比较刚硬或结构变形引起的气动力变化很小，例如超大展弦比无人机算例中，机翼弯曲变形后，气动载荷变化很小。这种情况下可以忽略fi(ε)/εi的影响，耦合迭代松弛因子λ取值在0.5～1.5之间，均可以获得较好的加速收敛效果，是静气动弹性耦合计算中经常遇到的情况。

2)fi(ε)/εi<-δ。这种情况下结构通常较软，结构变形偏离平衡位置，会诱导产生反向偏离平衡位置的计算结构变形。例如CHN-T1后掠机翼算例中，机翼弯曲变形诱导翼剖面负扭转，机翼载荷降低，结构变形计算产生负的弯曲变形，见图7(b)、图7(c)；或者机翼扭转刚度较低，力矩引起的扭转变形为主导形式，扭转角增加，力矩减小，计算的扭转角减小。对于此类情况，应当适当降低耦合迭代因子λ的取值(<1.0)，避免结构变形出现反复振荡或发散，如图6中λ≥0.75的情况。

3)fi(ε)/εi>δ。这种情况下结构通常较软，结构变形偏离平衡位置，会诱导产生正向偏离平衡位置的计算结构变形。例如前掠机翼，机翼弯曲变形诱导翼剖面正扭转角，机翼载荷增加，结构变形计算产生更大的弯曲变形；或者机翼扭转刚度较低，力矩引起的扭转变形为主导形式，扭转角增加，力矩增加，计算的扭转角继续增加。对于此类情况，可以适当增加耦合迭代因子λ的取值(>1.0)，提高耦合迭代的收敛速度。需要说明的是，增加耦合迭代因子提高收敛效率的前提是结构在气动载荷作用下能够稳定到静平衡位置，即静气动弹性耦合具有真实的物理解。

对于上述3类情况，λ的取值均不宜过小，当λ低于0.2时，λfi(ε)/εi项的影响变小，变形误差的增长率主要由1-λ决定，耦合收敛缓慢。另外，最优松弛因子λ的值与机翼形状、流动参数和结构刚度分布均有密切的关系。目前只能根据耦合类型给出松弛因子合理选择的范围建议，达到一定的加速收敛效果，避免数值迭代策略不当引起的耦合发散。定量给出最优松弛因子的取值，对于提升软件的使用体验，降低用户参数选择的难度，具有重要的工程应用价值，但仍然需要开展深入的研究和系统性的测试。

5.3 收敛影响因素

图8给出了不同初始值和不同增长率下的误差衰减曲线。可以看出，当误差衰减缓慢时(γ=0.95)，初始误差对收敛效率有一定影响，初始误差越大，需要更多的迭代步数达到收敛，而当误差衰减迅速时(γ=0.50)，初始误差对收敛效率的影响很小，不同初始误差达到收敛需要的迭代步数相差很小。

图8 误差衰减曲线Fig.8 Error damp curves

因此，针对静气动弹性耦合模拟，当耦合迭代松弛因子λ选取合理时，结构会迅速收敛至平衡位置，初次结构变形的大小对收敛效率的影响将很小。但初次结构变形受到解算器流场网格变形能力的限制，如果过大，容易引起网格破裂出现负体积，导致耦合模拟失败。所以对于特别柔性的机翼初次结构变形过大的情形，需要通过技术方式降低初次结构变形值，保证结构变形平稳收敛至平衡位置。