先进风洞是支撑飞行器自主研发，引领空气动力学及其相关学科创新发展的战略性基础设施^[1-2]。柔壁喷管是确保超声速风洞获得设计马赫数下均匀气流的重要部段，喷管型面质量的优劣很大程度上影响着风洞流场品质^[3]。在柔壁喷管静态调试阶段，可以通过激光跟踪仪等仪器对型面进行测量^[4]，从而实现型面调节。但是在风洞流场校测及运行阶段，喷管内流道无法放置测量设备，难以实时获取喷管型面信息，进而影响喷管流场调试效率与运行安全。开展柔壁喷管变形重构研究，实时精确感知型面信息，不仅有助于提高风洞流场品质，而且对喷管健康监测和安全运行具有重要意义。

为确保柔壁喷管运行安全，国内外早期超声速风洞多采用曲率保护装置限制柔性壁板弯曲最大曲率，如中国空气动力研究与发展中心(CARDC)的1.2 m超声速风洞^[5]和加拿大国家航空机构(NAE)五英尺三声速风洞等^[6]。随着新型传感器技术的发展，电阻式应变传感器、光纤传感器等被应用到风洞健康监测领域，聂旭涛等^[7]通过动力学仿真给出柔性壁板驱动位移对弯曲应力的影响关系。高川等^[8]在柔性壁板应力较大区域设置应变片，作为应力监测点，防止机构运动不同步导致柔性壁板应力过大。马列波等^[9]通过驱动机构同步监测、应变超阈值监测等手段，利用控制系统实现柔壁喷管安全监测。但是以上研究都是将传感器作为安全连锁控制手段，缺乏柔壁喷管型面变形重构相关研究。

国内外学者在柔性梁变形重构方面研究较多，NASA的William等^[10]提出了用于机翼变形重构的Ko理论，并发展出针对多种不同截面梁的变形重构方法。王寅等^[11]采用光纤光栅传感器在翼梁上测得应变数据，利用Ko理论研究了适用于柔性无人机机翼形变的测试方法。Ding等^[12]基于Ko理论研究了碳纤维层合板的变形重构问题。该方法基于梁的小挠度变形理论，认为梁沿轴线方向不发生变化，存在一定缺陷，不适用于大挠度变形重构问题。美国航空航天局(NASA)兰利研究中心基于最小二乘变分方程提出了逆向有限元法，张科等^[13]开展了基于逆向有限元法的机翼鱼骨变形重构研究，吴懋琦等^[14]研究了基于绝对节点坐标法的平面梁有限变形下变形重构问题。但是采用逆向有限元法需要构造出合适的单元，单元类型和节点数量影响变形重构的精度，因而存在单元构造复杂、计算效率不高等缺点。谭跃刚等^[15]针对薄板变形重构问题，结合应变与曲率、曲线与曲率的关系，提出了一种带误差补偿的曲率积分递推算法，但是该方法的重构精度依赖于分段数量和补偿因子的选取。

综上，目前已发展出多种基于应变监测的变形重构方法，为梁或板结构的健康监测提供了有利工具，但是也存在模型构造复杂、重构精度依赖于应变测点数量等问题。本文从柔性壁板应变-弯矩关系出发，针对超声速风洞多支点半柔壁喷管大挠度变形问题，建立型面组件力学模型，提出基于应变-弯矩的变形重构方法，推导柔性壁板型面曲线的微分方程，给出柔性壁板变形的椭圆积分表达式及解析解，并基于薄板变形理论对重构算法进行修正，在此基础上开展多支点半柔壁喷管的实验验证。

1 多支点半柔壁喷管力学模型 1.1 工作原理

如图 1所示，超声速风洞喷管内部气体流动分3个区域：亚声速段、喉部、超声速段。在亚声速段气流马赫数小于1.0，为使气流加速，截面积需要呈收缩状；气流加速至喉部，马赫数达到1.0；超声速段则需要截面积扩张，才能保证气流继续加速到给定马赫数。

图 2为典型的超声速半柔壁喷管三维模型，包括喉块、柔性壁板、驱动电缸、喷管框架、侧壁等部件。驱动电缸输出位移带动喉块和柔性壁板运动，从而成型出所需马赫数型面。

喉道与出口截面积比决定了喷管出口的马赫数，所以喷管喉块的位置对于出口的气流速度至关重要。柔性壁板弯曲成型出喷管的消波区，精确的型面可以极大降低喷管出口的气流脉动，提高流场品质。增加柔性壁板支撑点数量可以提高型面成型精度，因此多支点半柔壁喷管被广泛应用于国内外超声速风洞。

1.2 力学模型

图 3为单侧柔性壁板与喉块组成的型面组件机构简图，柔性壁板在喷管出口端与框架固定连接，自由端与喉块连接，柔性壁板与喉块背气流面设置驱动装置，受驱动装置强制位移约束。型面组件成型包括悬臂矩形薄板大变形和喉块刚体运动2类问题，应首先求解柔性壁板成型曲线，然后根据柔性壁板自由端坐标和转角，利用刚体平面运动理论得到喉块成型后的位置。

对于超声速风洞, 喉道下游随着截面扩张, 气流速度迅速增大, 动压上升, 型面组件所受气体静压降低, 本文忽略该区域内气体压力对柔性壁板变形的影响, 主要考虑支撑点的集中力和力矩作用。以柔性壁板固定点建立坐标系, 图 3中S₁~S₆为应变传感器, F₁~F₃为柔性壁板在各支撑点受力, M₁~M₃为力矩。

尽管柔性壁板厚度远小于长度与宽度, 属于薄板变形问题, 受横向剪切作用, 大型柔壁喷管存在沿宽度方向的翘曲问题。但是在柔性壁板靠近中轴线的核心区域, 相对两侧边缘处翘曲较小, 本文主要针对气流核心区域的柔性壁板变形开展分析, 因此假设沿柔性壁板宽度方向各处弯曲角度一致。在计算时可以先假设柔性壁板符合欧拉梁模型, 简化为多点支撑的悬臂梁模型并给出解析解, 然后基于薄板弯曲理论对抗弯刚度进行修正, 从而获得半柔壁喷管型面变形曲线。因此多支点半柔壁喷管柔性壁板可简化为图 4中的悬臂梁结构。

2 柔性壁板变形重构算法 2.1 应变-弯矩关系

根据欧拉梁弯曲的几何方程, 结构应变与曲率之间存在以下关系

(1)

依据平面假设与胡克定律, 在梁发生平面弯曲时, 曲线曲率与弯矩之间关系为

(2)

由(1)和(2)式可以得到弯矩与应变关系

(3)

式中：M(x)为梁结构微元所受弯矩；I为截面惯性矩；h为柔性壁板厚度。当获取柔性壁板上布置的应变传感信息后, 通过线性插值计算出柔性壁板弯矩分布。

2.2 大挠度变形推导

对于第i段柔性壁板, 根据悬臂梁模型假设, 各支点间梁的应变沿轴线呈线性分布, 其斜率为

(4)

式中：ε_i为应变值；x_i为第i个应变传感器对应的x坐标。

由(4)式可知第i段柔性壁板任意x位置的应变为

(5)

联合(3)和(5)式可得到第i段柔性壁板任意x位置的弯矩为

(6)

对公式(2)进行微分, 得到

(7)

(8)

两端同时乘以dθ得到

(9)

积分得到

(10)

对于第一段柔性壁板, S=l₁时边界条件为

可以求解出

(11)

对(10)式两端开平方处理, 并且根据梁变形的形状选取正负号, 梁变形形状下凹为负, 上凸为正, 然后进行积分, θ为0~θ₁; s为0~l₁, 得到

(12)

令, , (12)式的椭圆积分形式为

(13)

式中，, f=F(γ₁, t)-F(γ₂, t)，F为第一类不完全椭圆积分

求解公式(12)可以得到第一段柔性壁板末端转角θ₁, 对于柔性壁板任意一点, 转角为θ∈[0, θ₁], 对应坐标x和y为

(14)

(15)

式中，f′=F(γ′₁, t)-F(γ₂, t), e′=E(γ′₁, t)-E(γ₂, t), ；E(·)为第二类不完全椭圆积分。

对于之后的第i段, 可以考虑成旋转一定角度的悬臂梁来计算, 旋转角度为前一段自由端的转角, 为已知参数。模型如图 5所示, x_ny_no_n为第i段梁自身坐标系, 与全局坐标系xyo之间夹角为θ_i-1。

在自身坐标系下计算出梁弯曲曲线坐标x′_i, y′_i, 然后对该曲线进行旋转和平移, 组合后得到完整的变形曲线。

(16)

式中，x_i-1和y_i-1分别为前一段梁自由端转角和坐标。

2.3 基于薄板理论的抗弯刚度修正

根据薄板弯曲理论, 忽略板内部沿x和z方向的拉压力以及平错力的情况下, 薄板受外载荷弯曲变形后的微分方程为

(17)

(18)

式中：q为薄板单位面积所受载荷；D为抗弯刚度；μ为柔性壁板材料泊松比。

由(18)式可知, 与梁模型抗弯刚度EI不同, 薄板模型需要考虑泊松比影响。本文在不同柔壁喷管基础上，考虑了柔性壁板厚度与宽度比例因素, 构造修正系数β, 如(19)式所示。该系数体现了柔性壁板宽度、厚度和泊松比3个方面信息, 其物理含义是描述柔性壁板介于板壳和梁2种模型之间的程度。厚度和宽度比越小, 越接近板壳模型, 当厚度宽度比无限小时, 抗弯刚度与板壳模型一致。而厚度和宽度比越大, 泊松比对抗弯刚度的影响越小, 越接近梁模型, 考虑到柔壁喷管的柔性壁板厚度和宽度比一般不会大于1/5, 构造的修正系数仅对厚度和宽度比小于1/5的情况适用。将重构算法中抗弯刚度EI修正为βEI, 其中当厚度与宽度比小于1/5时，；当厚度与宽度比大于1/5时, 无需对抗弯刚度进行修正, 即β=1。

3 变形重构算法的仿真分析 3.1 与Ko理论对比

对图 4中结构，以n=2为例, M₁=0。利用ABAQUS有限元仿真软件开展两段悬臂梁结构变形数值分析, 求解器开启几何非线性选项。对比Ko理论与本文算法的变形重构效果。其中悬臂梁的结构参数如表 1所示, 各工况见表 2。

图 6a)为3种方法计算的变形曲线, 图 6b)为2种变形重构算法与有限元分析结果的偏差, 其中FEM代表有限元分析结果, Ko代表Ko理论计算结果, ELL代表本文算法结果。从图中可以看出, 在3种工况下, Ko理论和本文算法均与有限元仿真结果具有较好的重合性。其中在工况3下, 本文算法结果与有限元的最大偏差约为-0.014 mm, 为Ko理论偏差的21.13%。结果表明对于悬臂梁结构, 本文算法相比于Ko理论具有更加精确的变形重构能力。

3.2 修正效果分析

利用图 7中薄板模型对本文重构算法的修正效果进行验证, 长度尺寸l₁=10 mm, l₂=10 mm; 宽度尺寸b=10 mm; 选择2, 1, 0.5, 0.25 mm 4种厚度规格, 对应宽厚比分别为5, 10, 20, 40。支撑点1强制位移为D₁=1 mm, 支撑点2强制位移为D₂=3 mm。

选取板壳单元S4R划分网格, 网格尺寸为0.2 mm, 求解器设置中打开几何非线性选项。图 8为该薄板的有限元仿真应力分布结果和y向变形结果, 实线为固定约束位置, 虚线为对称约束区域, 也是要提取应变和变形的曲线。

图 9给出了4种规格薄板修正前后中心线的重构结果。从图中可以看出, 在不同宽厚比条件下, 与修正之前相比，修正后型面偏差均有明显降低, 重构型面与有限元仿真型面的偏差平均降低为修正前的2.97%, 修正后最大偏差约为修正前的10.78%, 重构误差平均从0.225 mm降低为0.006 mm, 型面重构精度明显提高。

4 实验与讨论

利用CARDC 0.3 m量级超声速风洞半柔壁喷管实验平台开展实验验证, 实验装置如图 10a)所示。该实验平台可实现声速至马赫数3.0范围的型面成型。柔性壁板设置4个驱动点, 喉块设置2个驱动点, 喷管出口为柔性壁板固定点, 柔性壁板与喉块之间采用螺钉连接。图 10b)为柔性壁板及背气流面传感器布置情况。

柔性壁板与支撑点连接处厚度为15 mm, 其余位置厚度均为10 mm, 因此本文假设柔板为10 mm均匀厚度。柔性壁板材料为07Cr17Ni7Al, 材料弹性模量为203.4 GPa, 屈服强度为982 MPa, 抗拉强度为1 274 MPa。

共计开展6种马赫数型面测试, 分别为1.0(声速型面), 1.3, 1.5, 2.0, 2.5, 3.0。采用TST5933便携式动态信号测试系统进行应变测试, 图 10为全行程工作时应变测量结果, S₁~S₁₀对应10个应变传感器。型面成型顺序为: 零位→马赫数1.3→马赫数1.5→马赫数2.0→马赫数2.5→马赫数3.0→零位→声速型面→零位, 在每个型面停留约20 s时间。

图 11为全行程应变测试结果, 可以看出, 马赫数1.0工况下, S₁~S₁₀均为负值, 即压应力状态。应变传感器S₃在马赫数1.5和2.0工况下出现轻微的压应力现象, 其余均为拉应力。柔性壁板在马赫数2.5工况下应变水平最高, 最大应变位置对应传感器为S₈, 最大值约为1.27×10^-3, 对应应力值约为258.32 MPa, 为材料屈服强度的26.3%。

以图 10b)中柔性壁板出口固定点为坐标原点, 建立空间坐标系。沿风洞逆气流方向为x轴正方向, 垂直于柔性壁板表面为y轴正方向, 垂直于x和y轴指向测量仪器侧为z轴正方向。

在风洞无气流情况下, 利用Leica AT960绝对位置激光跟踪仪测量喷管型面, 测量时由操作人员手持激光跟踪仪靶球, 紧贴柔性壁板及喉块内气流面移动, 获取靶球中心移动轨迹后, 再根据靶球半径对测量轨迹进行偏置, 得到被测型面。

图 12a)为实测型面与重构型面曲线对比, 图 12b)为二者偏差, 可以看出, 6种马赫数工况下, 本文算法与实测结果均可以良好吻合。其中柔性壁板范围内偏差在-1~0.5 mm。柔性壁板与喉块连接处最大变形达到80 mm, 为柔板厚度的8倍, 属于典型的大挠度变形问题。该位置变形重构结果与实测结果偏差约为-1.15 mm, 仅为变形量的1.44%。

表 3给出了6种马赫数型面下喉道位置偏差详细结果, 可以看出最大绝对误差发生在马赫数1.0工况下, 绝对值为3.56 mm; 相对误差最大发生在马赫数1.5工况, 约为喉道高度变形量的6.12%。

本文变形重构算法对超声速风洞半柔壁喷管柔性壁板力学模型进行了一定的简化, 包括应变分布的线性近似估计、柔性壁板出口刚性固定约束、柔性壁板等厚度假定等, 这些都会带来一定的变形重构误差。除此之外, 应变与型面测量误差也是影响变形重构偏差的一个重要因素。尽管存在一定的偏差, 但是本文算法变形重构结果与实验结果最大偏差仅为6.12%, 表现了良好的变形重构能力。

5 结论

1) 针对超声速风洞多支点半柔壁喷管，发展了基于应变-弯矩的大挠度变形重构方法。建立了半柔壁喷管力学模型，基于欧拉梁理论推导了柔性壁板大挠度变形微分方程和解析解，并利用薄板理论对模型进行了修正。

2) 采用有限元仿真手段初步分析了本文变形重构算法的重构效果，针对梁模型本文算法求解误差为Ko理论的21.13%；针对矩形薄板模型，修正后的变形重构误差降低为修正前的10.78%。

3) 利用0.3 m超声速半柔壁喷管实验平台开展实验研究，实验结果表明，在喷管喉道位置本文变形重构结果与实测型面最大偏差仅为结构变形的6.12%。与以往超声速风洞应力安全监测方式相比，本文方法极大提高了半柔壁喷管的健康监测能力。