风洞天平实时载荷监测方法研究

金 锋，徐 越，李 聪，徐传宝

(1.中国航空工业空气动力研究院，哈尔滨 150001；2.低速高雷诺数气动力航空科技重点实验室，哈尔滨 150001)

0 引言

风洞天平是风洞试验的核心测量设备之一，用来测量作用在风洞模型上面的气动载荷，它的安全与稳定是风洞试验顺利进行的前提[1-3]。目前针对风洞天平的研究多是从精准度方面下手，不断提高天平的性能指标，而安全方面则主要依赖于理论计算、有限元分析以及设计经验，预留比较大的安全系数，无法准确评估风洞试验过程中一些非常规因素(例如：振动等)给天平带来的影响。因此，风洞天平实时载荷监测方法的提出将有助于全面监测天平使用过程中的实时受载情况，降低危险情况的发生概率。

国外已经开展了一些风洞天平载荷监测方面的研究，法国宇航研究中心[4]基于应力叠加原理，提出了一种通过应力矩阵计算天平瞬时力和力矩的方法，并建立了专用的PC监控系统实时监测天平的真实负载(动态和静态)，取得了良好的应用效果。波音公司[5]于2004年开始积极开发基于应力的天平监测系统，并于2005年应用到跨音速风洞的使用中，该系统已经为波音风洞和美国宇航局艾姆斯研究中心的数十个试验项目提供了服务。国内方面，中国空气动力研究与发展中心[6]首次提出了需要对风洞试验过程中的一些关键试验参数进行监视测量，例如：风速、温度、压力、天平各元信号输出等。通过对这些关键试验参数的监视和测量可以有效的降低风洞试验风险，当某些参数超过设定的阈值时会自动报警提示。但也只是提出了一种思想却并未提供有效的监测手段。中国核电工程有限公司[7]提出了一种大批量标准部件在不同载荷下的应力计算评定方法，给出了详细的评定方法和可靠的实施路径，实现了对标准部件的真实应力评定，减少了后续评定的工作量。

本文基于应力分量的叠加原理，通过对天平实施单位载荷加载，计算得到应力影响系数矩阵，再通过天平六元信号输出叠加计算载荷监测区的Von-Mises应力值，通过与阈值进行对比判断出是否进行超载报警。为验证方法的准确性，分别开展了有限元法对比验证和静态校准对比验证。

1 载荷监测方法

风洞天平在使用过程中可能会受到较高的动态载荷或冲击载荷，从而使得天平受不可逆损伤，为了避免这种情况，就需要确定一种监测应力计算方法[8-10]来实时地计算天平选定节点的应力值，并与应力极限值作对比，以便风洞试验时实时判断天平是否超载。

本文所提出的风洞天平实时载荷监测方法的中心思想是分别计算6种单位载荷分别单独作用情况下，各选定节点的六个应力分量，每个节点对应36(6×6)个应力分量值。对于不同的加载载荷，36个应力分量值可以计算得出6个新的对应特定载荷的应力分量及MISES应力值，此时设定一个应力判断阈值就可以以此监测天平节点是否超载[11]。

先以二维平面为例，假设此二维平面单独受到某载荷A的作用，那么其某点W应力状态如图1所示；假设此二维平面单独受到某载荷B的作用，那么其某点W应力状态如图2所示；最后，假设此二维平面受到某载荷A与B的共同作用，那么其某点W应力状态如图3所示。

图1 单独受载荷A作用下应力状态

图2 单独受载荷B作用下应力状态

图3 受载荷A与B同时作用下的应力状态

因为始终处于线弹性，所以根据叠加原理可知，有：

(1)

(2)

(3)

假定测出天平中心所受载荷为：

X=Fx×1 N、Y=Fy×1 N、Z=Fz×1 N、L=Myz×1 NM、M=Mxz×1 NM、N=Myx×1 NM。则可推测出A点所受的应力状态(6个应力分量，监测应力值σA及Mises应力)为：

σx=σx1Fx+σx2Fy+σx3Fz+σx4Myz+

σx5Mxz+σx6Mxy

(4)

σy=σy1Fx+σy2Fy+σy3Fz+σy4Myz+

σy5Mxz+σy6Mxy

(5)

σz=σz1Fx+σz2Fy+σz3Fz+σz4Myz+

σz5Mxz+σz6Mxy

(6)

τxy=τxy1Fx+τxy2Fy+τxy3Fz+

τxy4Myz+τxy5Mxz+τxy6Mxy

(7)

τxz=τxz1Fx+τxz2Fy+τxz3Fz+

τxz4Myz+τxz5Mxz+τxz6Mxy

(8)

τyz=τyz1Fx+τyz2Fy+τyz3Fz+

τyz4Myz+τyz5Mxz+τyz6Mxy

(9)

(10)

因此，只要得到六种单位载荷分别单独作用时，对应的6个应力分量(6×6=36，一个点36个值)，就可以准确计算线弹性情况下某个指定节点在任意载荷下的MISES应力及应力分量。这种方法从理论上确保了监测应力值替代MISES应力来评估天平安全性能的可靠性和安全性。

需注意的是，一个点虽有36个系数，但是并不是在计算时都会用到。比如计算X=2 000 N时某点的应力，只需X=1 N时6个应力分量乘以2 000即得X=2 000 N时的六个应力分量。又比如计算Y=200 N，L=10 Nm时某点MISES应力，只需Y=1 N时六个应力分量乘以200，对应加上L=1 Nm时的六个应力分量乘以10的值，就可以得到最终的应力分量，最后再根据叠加原理计算得到该点的MISES应力值。

因此，对于一些特殊情况，基于应力分量的载荷系数计算方法是需要简化的。在6种载荷Fx，Fy，Fz，Myz，Mxz，Mxy作用下，应力分量及MISES应力如式(1)～式(10)所示。假设，6种载荷等比例扩大到K倍，即6种载荷为KFx，KFy，KFz，KMyz，KMxz，KMxy。此时的应力分量为：

σx'=σx1KFx+σx2KFy+σx3KFz+

σx4KMyz+σx5KMxz+σx6KMxy=Kσx

(11)

即：

σx'=Kσx

(12)

同理可得：

(13)

(14)

(15)

(16)

(17)

因此，监测应力值σA'及MISES应力为：

(18)

这就是基于MISES应力的载荷系数计算方法。因此，若知道基元载荷(Fx，Fy，Fz，Myz，Mxz，Mxy作用)下的MISES应力，就可以知道等比例加载到K倍后的MISES应力值。值得注意的是，若基元载荷为单个方向载荷，即只有一项存在时(单个载荷加载)，式(4～10)也是成立的。

同时，式(4)～(9)可写成矩阵形式。

{σ}=[S]{P}

(19)

其中：{P}为六元载荷向量，[S]为6阶方阵，其系数Sij表示第j个载荷单独作用引起结构某一点的第i个应力分量。上式的分量形式为

(20)

其中：系数S11，S21，……，S61为Fx=1，Fy=Fz=Mx=My=Mz=0时的6个应力分量，可通过有限元仿真分析获得。[S]的其他5列系数也是如此。

方便起见，以下称[S]为"应力影响系数矩阵"。该矩阵也可被视为一种特殊的柔度矩阵，反映了6个应力分量与六元单位载荷的传递关系。在实际应用时，可将实时测量到的六元载荷作为输入量{P}代入(20)，将其与[S]相乘后得到天平各处的应力响应，实现对风洞天平的实时载荷监测。

2 天平载荷监测系数提取

想要提取天平载荷监测系数，首先要确定需要重点监测的节点位置。一般有如下几个原则：

1)应变计粘贴位置，一般位于天平梁根部；

2)应力集中位置，一般位于天平梁与模型端或支杆端交界位置；

3)最大应力位置；

4)其他重点关注位置；

以直径40 mm杆式天平作为研究对象，使用ANSYS Workbench仿真平台对其进行网格划分，天平体网格大小为2 mm，四柱梁、支撑片等关键部位网格细化为1 mm，网格类型为Tet10四面体单元，单元数为313 181个，节点数为464 392个。基于上述重点监测节点提取原则，分别选取四柱梁根部、阻力梁根部、支撑片根部以及应变片粘贴位置作为载荷监测区域，如图4～7所示。在所选定的载荷监测区域中挑选若干个典型节点用于载荷监测法与其他方法进行对比分析。

图4 四柱梁根部监测点

图5 阻力梁根部监测点

图6 支撑片根部监测点

图7 应变片粘贴位置监测点

图8 Y元应变片粘贴处载荷监控点应力

图9 应变梁根部载荷监控点应力

图10 应变梁根部载荷监控点应力

利用ANSYS Workbench仿真平台分别计算六分量单元载荷作用下的节点应力信息，组合得到各应力监测点的应力影响系数矩阵。表1列出了其中位于四柱梁根部的第4 043号节点的[S]矩阵。

表1 四柱梁根部(4 043号节点)[S]矩阵

3 有限元法对比验证

有限元法[12]是上世纪中叶发展起来的一种通用性强、功能强大的数值分析方法。随着计算机技术的飞速发展，有限元法也迎来了高速发展期，广泛应用于工程结构、流体力学、电磁学、热力学等各个学科领域。有限元法是求解复杂数学问题的有效解决方法，它将复杂的求解域视为由很多个简单的单元组合而成，而这些单元之间通过若干个节点互相连接，巧妙地将原本求解复杂的连续体问题转变为求解相对简单的单元内部的近似函数问题，实现了复杂问题的简单化。因其在力学、声学、热力学和电磁学等领域的快速发展，我国也逐渐重视并发展有限元法在实际工程中的应用，目前已经全面应用于机械制造、土木工程、航空航天等诸多领域。

ANSYS软件[13-14]是一款功能强大、通用性强、集成度高的有限元计算和分析软件，同时也是世界范围内应用最广泛的一款商用有限元计算分析软件。用户可在ANSYS软件上实现实体建模以及有限元网格划分，同时可以调用其分析计算模块实现对力学、电磁学、声学等问题的求解，最后通过后处理模块实现对计算结果的可视化展示，例如：图表、曲线、梯度显示等。其中，ANSYS Workbench仿真平台能对复杂机械系统的结构静力学、结构动力学、刚体动力学、流体动力学、结构热、电磁场以及耦合场等进行分析模拟。本文所有有限元计算都是在ANSYS Workbench仿真平台上完成。

基于上一节所建立的有限元模型以及所选取的载荷监测点，使用ANSYS Workbench仿真平台进行有限元计算分析[15-17]，得到了所有载荷监测点的矩阵[S]。基于此分别开展单元加载验证、等比例加载验证以及混合加载验证等工作。

1)单元加载对比验证：

以Y元为例，分别使用ANSYS Workbench仿真平台和载荷监控系数法计算500 N升力作用下的天平(升力测量应变片粘贴位置)应力输出，对比其结果如表2所示。

表2 载荷为Y=500 N时,应力值对比(抽取5点)

由表2中结果可知，单元载荷加载情况下，Workbench计算结果与载荷监测系数法计算结果吻合情况良好，最大差值小于0.03%，说明单元加载情况下载荷监测系数法具有较高的计算精度。

2)等比例加载对比验证：

对天平同时施加六元载荷：X=50 N，Y=50 N，Z=50 N，Mx=50 Nm，My=50 Nm，Mz=50 Nm。分别使用ANSYS Workbench仿真平台和载荷监控系数法计算天平(应变梁根部)应力输出，对比其结果如表3所示。

表3 等比例载荷作用时应力值对比(抽取5点)

由表3中结果可知，等比例载荷加载情况下，Workbench计算结果与载荷监测系数法计算结果吻合情况良好，最大差值约为0.03%，说明等比例加载情况下载荷监测系数法具有较高的计算精度。

3)非等比例载荷加载对比验证：

对天平同时施加六元载荷：X=500 N，Y=1 800 N，Z=800 N，Mx=70 N·m，My=70 N·m，Mz=100 N·m。分别使用ANSYS Workbench仿真平台和载荷监控系数法计算天平(应变梁根部)应力输出，对比其结果如表4所示。

表4 非等比例载荷作用时应力值对比(抽取5点)

由表4中结果可知，非等比例载荷加载情况下，Workbench计算结果与载荷监测系数法计算结果吻合情况良好，最大差值约为0.04%，说明非等比例加载情况下载荷监测系数法具有较高的计算精度。

至此可以得到结论，在各种不同的载荷条件下，使用基于应力分量的载荷系数计算方法得出的MISES应力(监测应力值)都与有限元法计算得出的MISES应力值近似相等，这种方法可以很稳定的衡量天平是否超载。

4 静态校准对比验证

天平静态校准验证与一般的天平静态校准过程一致，只不过将校准过程中的实际加载值与电压信号输出数据提取出来用来与基于应力分量的载荷监测系数法计算得出的MISES应力(监测应力值)进行对比验证。具体实施过程如下。

天平静态校准[18-19]使用BACS500自补偿静校台，如图11所示。可实现0～5 000 N载荷量级的杆式应变天平的高准确度静态校准。使用空气静压轴承滑轮和激光全息透镜分光技术结合激光漂移补偿技术实现加载系统与复位系统的组合控制。系统总体精度为万分之三。

图11 BACS500自补偿静校台

天平各元最大校准载荷：X=±1 882.9 N，Y=±470.7 N，Z=±784.5 N，Mx=±64.7 Nm，My=±64.7 Nm，Mz=±94.1 Nm。

校准环境：温度18 ℃，湿度20%。

供桥电压：10.005 V。

以单分量加载为例，当天平升力收到Y=100 N的载荷作用时，此时升力测量通道将会有一个电压信号输出△U，考虑到供电电压U、应变计灵敏度系数K、桥路叠加系数N(一般情况N=1或2)，由此可以计算得到Y元应变计此时的应变值为：

(21)

根据材料力学[20]理论可以计算得到应变计粘贴位置应力值为：

σ=E·ε

(22)

于是将Y=100 N代入到天平载荷监测系数矩阵中，可以得到单元载荷Y作用下的六分量应力。此时需要注意的是，由应变计信号输出所计算得到的应力值是具有方向的，而不是此处的MISES应力值。例如，施加Y向载荷，一般Y元应变计粘贴方向为X向，则此时应该使用载荷监测系数法计算得到的σx与由应变计输出计算得到的应力σ进行对比验证，如果结果一致则说明本方法是有效的。

天平静态校准采集的是各个电桥的电压输出值，而载荷监控系数法计算得到的是监测点的应力值，因此需要确定二者之间的对应关系才能验证本文方法计算结果的准确性。

首先明确桥路与单个应变片输出的对应关系，而单个应变片可以近似为多个载荷监测点应变输出的平均值，由此可以找到桥路与监测点的对应关系。杆式天平桥路选择为全桥接线法[21]，如图12所示，其基本原理如下。

图12 全桥接线法

电桥的4个桥臂为四片粘贴在两个应变梁上测点处的电阻应变片R1、R2、R3、R4，即4个桥臂都是工作片，且电阻应变片的规格型号相同(应变片的灵敏系数、电阻值和丝栅尺寸均相同)。

当构件受力变形时，电阻应变片的电阻值变化分别为R1+ΔR1，R2+ΔR2，R3+ΔR3，R4+ΔR4。在小变形条件下，ΔR/R<<1，且R1=R2=R3=R4，根据两点压差与电阻的比例关系，则有：

(23)

(24)

UA-UC=Eu

(25)

整理上式，得：

(26)

上式分子分母除上R1，且R1=R2=R3=R4，并忽略二阶以上微量，可得：

(27)

(28)

式中，ΔUBD为电桥对角点BD输出电压变化量；Eu为电桥对角点AC接电源的电压。

设电阻应变片灵敏系数为K，则

(29)

由于4个桥臂的电阻应变片的灵敏系数K值相同，则上式可改写为：

(30)

其中：ε1、ε3、ε2、ε4分别为电阻应变片R1、R2、R3、R4产生的应变值，该值为代数值，其正负号由测点的变形决定。电阻应变仪上的示值为：

ε仪器=ε1+ε3-ε2-ε4

(31)

全桥连接法比单片测量灵敏度提高了4倍，应变仪的读数是测点实际应变的4倍。因此，测试点的实际应变值为：

ε实际=ε仪器

(32)

考虑到应变片粘贴位置位于应变梁对称处，各个应变片输出近似相等，因此可以得出：

ε实际≈ε单片

(33)

由于该天平实际校准数据过多，选择各元最大单元加载载荷作为对比数据，结合公式(21)，校准数据如表5所示。

表5 静态校准数据

根据各单元加载时天平梁受力变形情况，选取不同位置的载荷监控点。分别施加单元载荷，结合公式(22)，载荷监测系数法计算得到各应力监控区域的平均应力和平均应变，如表6所示。

表6 载荷监控系数法计算结果

在各元载荷单独作用下，本文方法与静态校准结果的对比如表7所示。

表7 载荷监测系数法与静态校准结果对比

由表7可得，六元单元加载情况下，载荷监控系数法计算得到的应力监测值与静态校准采集得到的数据差值均小于10%，其中最大差值出现在阻力元，差值为9.58%，其余各元误差值均小于4%。

计算结果与静态校准结果存在差异的主要原因在于：

1)理论计算模型是理想化的数学模型，其约束条件和加载方式不可避免地与真实校准存在差异。

2)有限元法的网格带有一定的厚度，其体积不为零。由于算法构造上的固有特点，仿真输出的是单元高斯积分点经外推后得到的结果，而非表面应变。因此，仿真解与应变片测出的结果在本质上存在微小的差别。

3)从校准角度分析，由于应变片尺寸无法做得很小，静态校准反映的只是受力后物体表面的宏观变形，无法像仿真那样描述各点的变形不均匀特性，这也是结果存在差异的客观原因之一。

由此可以认为载荷监控系数法与静态校准结果相当，可作为风洞天平载荷监控预警的理论依据。

5 结束语

本文针对风洞天平缺少有效的安全监控手段的现状，提出了一种基于应力叠加原理的风洞天平实时载荷监测方法，实现了对风洞天平全方位的实时载荷监测，大幅度降低风洞试验出现天平损坏的风险。通过本文研究，得到结论如下：

1)基于有限元法计算得到风洞天平危险区域载荷监测点的应力系数矩阵，为实现全域载荷监控奠定基础。

2)载荷监测系数法与有限元方法计算得到的载荷监测值误差小于0.02%，证明该方法具有较高的理论计算精度。

3)载荷监测系数法与静态校准得到的结果对比，阻力元误差小于10%，其余各元误差小于4%，证明该方法具有较高的实际应用价值。