圆柱贮箱内液体非线性稳态晃动实验及动力学特性分析

吴文军，高超南，岳宝增，邓明乐

(1. 广西科技大学机械与交通工程学院，柳州 545006；2. 广西汽车零部件与整车技术重点实验室，柳州 545006；3. 北京理工大学宇航学院，北京 100081；4. 中国空间技术研究院，北京 100094)

0 引 言

液体晃动问题广泛存在于大型充液航天器、液体火箭、带贮液罐类重型车辆和船舶等运载工具中[1-7]，贮箱内液体晃动产生的动态晃动力和力矩会严重影响运载工具操纵系统或控制系统的任务执行效率、控制精度和整体稳定性，长久以来对于液体晃动问题的研究均备受国内外大量学者们的青睐。

现有的液体晃动分析方法主要有理论解析法、数值仿真法和实验研究法。其中，理论解析法需要依据壁面边界条件推导具备正交性的模态坐标基后，利用液面处的运动边界条件和动力学边界条件，建立模态坐标函数系数的状态微分方程和求解体系。理论解析法具有计算效率高，便于研究不同重力环境和荷载工况下的液体晃动本质特征问题，分析结果可作为参考依据，用于指导实验研究和实际工程应用，但理论解析法仅能适用于矩形箱、平底圆柱箱、球箱和卡西妮贮箱等形状规则贮箱[8-12]。数值仿真法指利用计算流体力学(Computational fluid dynamics, CFD)理论研究液体晃动问题，目前已开发出了多个相应的CFD商业软件，例如：Fluent、Flow3D、Adina等。数值仿真法具有良好通用性和适用性，可研究不同贮箱尺寸、不同液体属性和贮箱内设置防晃板等更难现实的晃动问题[13-19]。理论解析法和数值仿真法可统称为理论分析方法，但由于液体晃动状态和行为特征受到液体属性、外界扰动、贮箱的几何形状、防晃板类型和计算算法的不足等多种复杂因素的影响，致使目前针对贮箱内的液体晃动动力学特性的所有理论分析方法中，尚无精确、可靠和高效的计算方法能得到普及应用，液体晃动问题已普遍成为相关工程领域技术发展亟待解决的重要热点之一。

实验研究法作为分析液体晃动问题最直观可信的研究方法，不仅可以有效避免理论分析和数值模拟存在的不足，而且可以利用准确、有效的实验数据进一步对理论分析和数值模拟进行验证与修订，因此通过实验方法对研究液体晃动动力学特性具有较为重要的现实意义。在液体晃动问题的实验研究方面，文献[20]最早通过理论结合实验，研究了柱形和球形贮箱的液体晃动问题，得到了液体晃动对贮箱结构的影响以及晃动频率与隔板参数之间的关系。文献[21]研究了运载火箭中椭球形贮箱液体晃动动力学特性，并叙述了晃动试验中观测到的非线性效应。文献[22]将有限元仿真计算结果与实验测定结果进行对比，研究了储箱内无隔板和有隔板情况下黏性流体晃动的阻尼机理。文献[23]对半球形容器中液体自由晃动的非线性行为进行了实验研究，得到了频率与阻尼比的时间历程关系。文献[24]对贮箱内水平带孔挡板进行了实验研究，通过改变实验工况，详细观察了边壁处的自由液面波动以及晃动的共振频率，讨论了挡板的晃动抑制效果。文献[25]利用振动台试验测试了水箱中液体的动水压和波高衰减数据，计算求得液体的晃动频率和阻尼比，并将试验结果与数值分析、等效模型理论计算结果进行比较。文献[26]基于透射散斑相关性测量方法并根据斯涅尔定律计算出液体表面变形，结合多向牛顿迭代算法计算出液面的高度，并验证了该方法的可行性和有效性。文献[27]利用参数振动的原理，激发各阶稳态晃动模态后，通过激光测量液体表面波高的自由衰减曲线，得到不同容器内液体各阶自由晃动频率与阻尼比系数。文献[28]研究了圆柱形储液罐在水平地震激励下晃动模态的动力学特性，并通过实验观测得到液体晃动的占优模态。文献[29]通过振动台对充液容器施加横向简谐激励，观察自由液面的运动形式，并采用压力传感器测量容器侧壁某点处的液体压力时间历程。文献[30]采用两种有限元仿真计算方法分别与试验测量结果进行了对比分析，校验了两种仿真方法的有效性，有效指导了工程试验。文献[31]通过实验测定了旋转晃动对储罐不同位置应变和加速度发展的影响。

上述研究文献在分析液体晃动模态坐标、晃动阻尼和动力学特性等方面已取得了大量丰富的重要研究成果，但大多数文献的研究往往忽视了影响运载工具整体性能的实际主导因素——液体晃动时产生的动态晃动力和力矩。鉴于此，本工作将以晃动力和晃动力矩的高效、准确测量为基本目标，依据液体晃动动力学基础理论，基于前期工作基础和现有实验条件，利用力传感器、加速度传感器、信号发生器、功率放大器、激振器和动态信号采集仪等仪器设备，设计一种液体晃动力和晃动力矩实时采集与分析的实验平台，同时开展大量的圆柱贮箱内液体晃动实验研究。将实验测量结果与CFD软件的仿真结果进行对比分析，并依据研究结论较客观的评价了两种分析手段的适用性、可行性、有效性和准确性。以期为进一步实现液体晃动动力学特性的在线实时测量、荷载反馈、晃动状态识别、晃动抑制和运载工具稳定性控制等相关动力学与控制问题的研究提供一定的参考数据和实验研究基础。

1 液体晃动动力学实验装置和实验原理

基于作者们的前期工作基础[32]，文中所设计的液体晃动动力学实验平台装置由信号发生器、功率放大器、激振器、应变放大器、动态信号采集仪和自主设计的贮箱夹紧装置等仪器设备构成，所需设备的型号和数量具体见表1。

表1 设备型号和数量规格表Table 1 Table of device models and quantity specification

为了最大化的采集和分析贮箱内液体非线性稳态晃动行为的表现特征和动力学特性，激发液体晃动的各阶模态，所设计的液体晃动实验平台装置既能实现充液贮箱结构固定安装，又能实现最大化的测量贮箱内液体晃动产生的非平衡晃动力和晃动力矩，其三维示意图、采集原理示意图和实际照片分别如图1(a)、(b)和(c)所示。其中，图1(b)中的坐标原点为力传感器安装位置的参考中心点，激振器的激励方向为X轴；d1和d2分别为坐标原点至安装在X轴和Y轴力传感器的水平距离；F1，F2和F3分别表示测量晃动力的三轴力传感器沿三个坐标轴方向采集的数据；F4，F5和F6分别表示测量晃动力矩的三个单轴力传感器采集的数据。

图1 液体晃动实验平台装置的安装方式Fig.1 Installation method of liquid sloshing experimental platform

由于实验过程中传感器采集的数据同时包含有液体晃动产生的动态力和实验装置固件(主要包括贮箱、晃动支撑板、夹具、球铰等)运动产生的惯性力，为了得到真实的液体晃动动态力和力矩，在实验数据的处理中需要将其消除。图2给出了固件质量与质心距离示意图，其中，Fx和Fs分别为液体晃动等效力和实验中采集的固件惯性力；h1和h0分别为贮箱底面和固件质心距直角坐标原点的垂直距离。为了得到真实的晃动力及晃动力矩，需首先识别出固件的总质量m0和质心距离h0。当贮箱内未充液时，依据动力学平衡关系，分别有:

图2 充液贮箱固件质量与质心距离示意图Fig.2 Schematic diagram of mass and centroid height of liquid filled tank

(1)

式中:ax为实验中采集的加速度；e和ε分别代表不同的综合测量误差，与传感器动态性能有关，可设其平均值为零。根据最小二乘法原理，若令:

(2)

将式(2)求导并取零值后，最终可得固件总质量m0的无偏估计量为:

(3)

同理，根据最小二乘原理，可得质心距离h0的无偏估计量为:

(4)

根据上式识别出的mo和h0，可得贮箱内液体晃动时沿水平坐标轴产生的晃动力及晃动力矩分别为:

Fx=F1-m0ax

(5)

Fy=F2

(6)

Mx=F4d1

(7)

My=F5d2-m0axh0-Fxh1

(8)

2 液体晃动实验分析与研究

2.1 动态测量系统可靠性检验实验

为了准确稳定的测量液体晃动产生的动态力和力矩信号，实验中需尽量避免实验设备工作状态的不稳定和人为操作不当影响实验结果，因此实验过程严格按以下标准步骤进行：实验设备预热、测量通道参数设置、滤波器设置、采集系统平衡清零和数据后处理及分析。

首先，为了校核测量系统及相应实验数据采集、分析与处理方法的准确性和可靠性，利用激振器、加速度传感器和铝块等设计出如图3所示的实验装置，通过对激振器输入不同频率的正弦信号，利用测得的加速度和轴力值，进一步运用最小二乘法识别出铝块的质量，从而实现动态测量系统可靠性的检验。

图3 动态校核实验装置示意图Fig.3 Schematic diagram of dynamic verification

依据牛顿第二定律，轴力fi与铝块质量m和加速度ai的线性关系为:

fi=mai+ei

(9)

式中:ei代表平均值为零的综合测量误差。同理，依据最小二乘原理，易推导出铝块质量m的识别公式为:

(10)

以一个质量为1.003 kg铝块为例，开展质量识别的实验，表2给出了动态测量系统可靠性检验的实验数据。由实验结果可知，该动态实验均能较准确的识别出铝块质量，且识别误差均在2%以内，表明该测量系统可以得到较精确的测量结果，具有可靠的测量性能，为后续实验的精确测量提供了有力保证。

表2 动态测量系统可靠性检验的实验结果Table 2 The experimental results of reliability testing for the dynamic measurement system

2.2 液体小幅晃动固有频率实验验证

以一个内径R=0.09 m的平底圆柱玻璃贮箱为例，通过设置信号放大器输出电流最大幅值保持不变，输入信号的频率范围从0.5 Hz～5.0 Hz进行扫频分析，分别得到当静液面高度h=1.0R、1.5R和2.0R时的共振频率和自由液面晃动模态。经多次重复实验验证，取均值后确定液体晃动的一阶固有频率。表3列出了本文实验结果与文献[2]中理论公式、CFD软件仿真的对比结果。由于在理论公式和CFD软件仿真中，均未能真实考虑贮箱尺寸加工精度、材料表征性能和液体属性等参数的不确定性，因此理论计算值和CFD软件仿真结果均相比实验值明显偏高。

表3 理论值与实验值相对误差表Table 3 Relative error table of theoretical value and experimental value

2.3 液体晃动实验与仿真对比分析

以液体静液面高度h=1.5R时为例，开展实验结果与CFD软件仿真结果的对比研究。参考贮箱内液体晃动的模态分析结果，设置频率变化范围为0.5～4.0 Hz，以贮箱内的液体产生非线性稳态晃动(液面不产生破碎或无规则大幅晃动)为准则，设置相应的激振器激励幅值，通过多组独立反复的定频驻留实验分析，得到多种类液体非线性晃动行为特征。同时，以液体晃动实验过程中采集到的贮箱加速度作为输入，利用CFD软件对相应的晃动问题进行仿真计算，并将仿真结果与实验数据进行对比分析，用于验证实验数据测量的准确性和CFD软件在分析贮箱内液体晃动问题时的有效性和精确性。

因篇幅有限，文中仅给出了激励频率在非旋转晃动和旋转晃动两种状态下液体晃动动态晃动力和力矩的时程曲线对比，具体结果分别如下图4～6所示。

由图4可知，因实验装置各部件间不可避免的存在细微的安装间隙，以及实验平台运动过程中摩擦力的影响，实验采集的数据相比CFD仿真结果始终存在高阶非线性项，但在此类非旋转晃动状态下，实验测量值与CFD仿真结果基本一致，有效验证了实验数据测量、分析和处理方法的准确性和可行性。

图4 当激励频率为1.92 Hz时的非旋转晃动Fig.4 Non-rotational sloshing at excitation frequency of 1.92 Hz

然而，通过对比图5、图6中所给出的旋转晃动状态下晃动力及晃动力矩时程曲线，却可明显看出CFD的仿真结果与实验值相差较大。为了避免因计算误差而影响结论的正确评价，经作者们多次实验验证和CFD仿真计算后，依然发现：在所有可能激发贮箱内液体产生旋转晃动的频率范围内，实验结果中均未出现CFD仿真中始终存在的“理想节拍旋转晃动”现象，基于动力学特性分析原理，可知其主要原因是CFD软件对液体晃动阻尼相关属性的描述与真实情况存在偏差，液体旋转晃动时的整体晃动系统阻尼比小于真实值从而导致仿真计算结果偏差。由此可见，虽然CFD软件可以研究一般情况下的液体晃动行为，但当液体晃动行为表现为旋转状态时，由于自身算法误差和液体属性描述的不完全真实性，导致了CFD仿真无法准确的分析贮箱内液体旋转时的液体晃动力学特性，甚至会得到错误的计算结论。因此，若单纯独立依靠CFD仿真研究贮箱内的液体问题，将不能准确揭示其真实的动力学特性，必须针对实际情况，综合利用理论解析法、数值仿真法和实验研究法等开展相关问题的研究。

图5 当激励频率为2.14 Hz时的旋转晃动Fig.5 Rotational sloshing at excitation frequency of 2.14 Hz

图6 当激励频率为2.23 Hz时的旋转晃动Fig.6 Rotational sloshing at excitation frequency of 2.23 Hz

2.4 液体非线性晃动行为表现特征研究

鉴于CFD数值仿真未能真实分析液体旋转晃动状态下的晃动力和晃动力矩，所以本小节将着重通过实验手段进一步细致研究液体非线性稳态晃动行为特征随着激励频率的变化规律。

为了能够准确分析和捕捉晃动行为表现特征及相应的晃动力、晃动力矩实验数据，在实验过程中，当频率不处于临界频率附近时，实验采取的激励频率以0.1 Hz的增量改变；而当频率处于临界频率附近时，激励频率则以0.01 Hz的增量改变。同样，经多次重复实验，通过数据分析与概率统计，最终归纳总结得到液体非线性晃动行为特征随着激励频率的变化规律，具体见表4。

由表4可知，随着激励频率的变化，贮箱内的液体晃动行为表现特征按占优晃动模态可大致分为一阶反对称模态占优的低阶状态、复杂旋转状态、一阶或二阶反对称模态占优的混合状态和高阶模态占优的高阶状态。而液体晃动行为表现特征随着外激励频率的变化，将分别呈现如下几类晃动行为:

表4 液体晃动行为特征变化规律统计Table 4 Statistics of liquid sloshing behavior characteristics

1)小幅或有限幅驻波晃动行为

当激励频率较明显的小于2.07 Hz时，在具有一定频率的外激励作用下，贮箱内的液体仅表现为小幅晃动，图7给出了当外激励频率等于1.80 Hz时，液体晃动产生的晃动力、晃动力矩。

图7 当外激励频率为1.80 Hz时的晃动力和晃动力矩Fig.7 Sloshing forces and moments when the external excitation frequency is 1.80 Hz

随着频率的增加，液体晃动波高、晃动力和晃动力矩均逐渐增大；而当激励频率接近液体一阶固有频率时，液体在开始阶段就发生明显的共振现象，图8给出了共振时液体晃动力及晃动力矩时程曲线，随着时间的推移，晃动力和晃动力矩最终均将趋于稳定的有限幅驻波晃动。

图8 共振时的晃动力和晃动力矩Fig.8 Sloshing forces and moments at resonance equipment

2)多类非线性行为交替变换的旋转晃动

当外激励略大于2.07 Hz时，自由液面的运动初期就展现为明显的非线性驻波晃动行为，随着晃动幅值的变化，逐渐呈现出非线性驻波晃动、有限幅的顺时针旋转与逆时针旋转交替的晃动行为；当频率进一步增大时，顺时旋转运动的持续时间变长，逆时针旋转时长逐渐减少，直到逆时针旋转完全消失，典型的多类非线性行为交替变换的旋转晃动产生的晃动力、晃动力矩曲线如图9所示(注：图中“非”代表非线性驻波晃动；“顺”代表顺时针旋转；“逆”代表逆时针旋转，文中以下图中均依此方法进行标注)，由图中晃动力Fx、Fy曲线特性可知，液体晃动随着时间的推移将逐渐趋于稳态旋转晃动。

图9 当外激励频率为2.09 Hz时的晃动力和晃动力矩Fig.9 Sloshing forces and moments when the external excitation frequency is 2.09 Hz of liquid filled tank

3)非线性驻波与顺时针旋转交替变换的旋转晃动

当外激励大于2.15 Hz时，逆时针旋转完全消失，液体自由液面表现为有限幅非线性驻波晃动与顺时针旋转交替的晃动行为，但随着时间的推移，非线性驻波晃动逐渐消失，晃动行为逐渐表现为稳态顺时针旋转运动，相应的晃动力和晃动力矩的曲线特性如图10所示。

图10 当外激励频率为2.16 Hz时的晃动力和晃动力矩Fig.10 Sloshing forces and moments when the external excitation frequency is 2.16 Hz

4)一致顺时针稳态旋转晃动

外激励频率在2.19～2.28 Hz范围内液体自由液面表现为顺时针旋转共振的晃动行为，开始阶段晃动幅值和旋转速度逐渐增大，但随着时间的推移，两者均将趋于稳定状态，而稳定状态下的晃动波高随着外激励频率的增加而减小，但旋转速度随着外激励频率的增加会先增大后减小。典型的一致顺时针稳态旋转晃动力曲线特性如图11所示。

图11 当外激励频率为2.23 Hz时的晃动力和晃动力矩Fig.11 Sloshing forces and moments when the external excitation frequency is 2.23 Hz

5)小幅或有限幅驻波晃动及有限幅拍振晃动行为特征

液体晃动实验过程中发现：当外激励频率大于2.28 Hz时，液体晃动仅出现小幅晃动，晃动幅度逐渐趋于平缓。且大多数情况下，晃动力和力矩在开始阶段均表现出拍振行为特性，但是由于存在边界阻尼耗散和液体黏性阻尼耗散，致使液体晃动出现的拍振现象只能维持一段时间，并最终趋于稳定的非线性驻波晃动，图12给出了当外激励频率为2.50 Hz时的晃动力和晃动力矩。

图12 当外激励频率为2.50 Hz时的晃动力和晃动力矩Fig.12 Sloshing forces and moments when the external excitation frequency is 2.50 Hz

6)高阶模态占优的小幅晃动行为特征

当激励频率大于3.5 Hz时，液体自由液面表现为对称的高阶晃动模态，且晃动幅值均较小。随着频率的继续增加，峰值会发生多点破碎，自由液面内出现多个波峰和波谷的高阶模态占优的小幅晃动行为。

2.5 液体非线性晃动动力学特性研究

取静液面深度h=1.5R时进行实验研究，由于晃动力及晃动力矩在0.5～1.5 Hz范围内变化微弱，所以仅给出激励频率在1.5～4.0 Hz范围内变化时液体晃动产生的最大稳态晃动力和晃动力矩的晃动动力学特性，结果如图13所示。

图13 不同频率外激励产生的最大晃动力和晃动力矩Fig.13 Max sloshing forces and moments under different frequencies

研究表明：1)当液体晃动处于低阶状态时，最大晃动力和晃动力矩均随着外激励频率的增加而逐渐增大，对运载工具的控制系统产生一定的影响;2)当液体晃动处于复杂旋转状态时，自由液面内液体旋转方向和晃动方式均与激励幅值大小无关，幅值很小的持续激励也能引发非线性旋转晃动。由于液体旋转晃动状态剧烈且复杂，贮箱内有更多的液体参与晃动，自由液面表现为剧烈的非线性大幅晃动，致使该状态下晃动力和晃动力矩稳态峰值明显增大；晃动力和晃动力矩的最大值均不随旋转状态的改变而发生变化，与文献[21]得到的实验现象和液体晃动特性高度一致，复杂的旋转晃动将严重影响运载工具控制系统的任务执行效率和控制精度，甚至会导致整体系统的失稳;3)当液体晃动处于混合状态时，液体晃动波高幅值逐渐趋于平缓，随着外激励频率的增加，晃动响应先减小后又逐渐增大，对运载工具控制系统的影响较小;4)当液体晃动处于高阶状态时，随着外激励频率增加，整体平台的加速度值将增大，此时虽然液体晃动高波的幅值较小，但贮箱整体近似于刚体运动，由于惯性力的存在，相应的力和力矩依然会增大，几乎不会影响运载工具控制系统的执行效率和控制精度。

3 结论与展望

本文设计了一种用于测量贮箱内液体动态晃动力及力矩的实验平台装置，实现了贮箱内液体三维晃动力和晃动力矩动态信号的实时、准确测量，并通过大量实验，研究了不同工况下液体的晃动动力学特性。实验结果表明：1)由于贮箱尺寸的加工精度、结构材料表征性能、液体属性等参数不确定性的影响，理论公式和CFD仿真计算得到的液体固有频率均偏高;2)虽然CFD软件可以研究一般情况下的液体晃动行为，但当液体晃动行为表现为旋转时，由于自身算法误差和液体阻尼描述的不完全真实性，导致了CFD仿真无法准确的分析贮箱内液体旋转时的液体晃动力学特性，甚至会得到错误的计算结果，不利于揭示贮箱内液体真实的晃动特性;3)即使在单一横向平动外激励作用下,当外激励频率在一定范围内变化时，液体自由液面表现出多种类型的非线性晃动行为特征，其中液体非线性旋转晃动时的行为特征最为显著，相应的晃动幅值和旋转方式与激励幅值无关，幅值很小的持续激励也能引发非线性旋转晃动，实际工程中，这将更容易影响运载工具中动力学控制系统的任务执行效率和控制精度。

基于本文所设计的液体晃动实验平台和研究方法，可进一步深入开展如下更具现实意义的液体晃动研究：1)在研究经费、实验场地和设备性能等条件允许的前提下，可以较完善的模拟和测量具有多种类几何外形贮箱、液体属性以及防晃板类型等不同情况下的液体晃动力及力矩，便于更真实的研究液体晃动动力学特性；2)能够有效避免理论分析存在的不足，更有利于揭示液体晃动的本质特征，为研究液体晃动系统计算模型的建立、修正，以及等效力学模型的建立和旋转晃动的抑制方法研究，降低液体晃动对运载工具的影响研究提供重要参考和技术支撑；3)为进一步实现充液贮箱内液体晃动动力学特征的在线实时反馈、液体晃动系统参数辨识、液体晃动状态识别、液体晃动抑制和运载工具整体系统的动力学与控制等问题的研究提供可靠的测量手段和实验基础。