王佐 吴文军 姜丙超 高超南
摘 要：为解决传统等效力学模型描述贮箱内液体晃动特性时在准确性和完整性等方面存在不足的问题，首先，运用势流理论推导出部分充液圆柱贮箱的动力学状态方程，依据方程中液体晃动速度势函数与晃动波高函数中傅里叶-贝塞尔级数展开系数的正交性条件，将球形贮箱内液体晃动的求解区域向与之外切的圆柱贮箱求解区域内进行扩展，建立用于描述球形贮箱内液体晃动问题的参数化模型;然后，通过CFD数值模拟实验得到球形贮箱受横向简谐激励时小幅晃动情况下的液体晃动力和力矩;最后，利用Matlab灰箱辨识函数对所建立的参数化模型中的待定系数进行准确辨识，实现了球形贮箱内液体横向晃动问题的高精度等效建模。结果表明：所建立的参数化等效力学模型能较准确和完整地描述球形贮箱内的液体晃动固有频率、晃动力和晃动力矩等动力学特性，辨识精度均能达到92%以上，为进一步在线实时辨识与实际工程应用提供参考。
关键词：球形贮箱;液体晃动;参数化模型;CFD数值模拟;灰箱辨识;系统辨识
中图分类号：O353.1 DOI：10.16375/j.cnki.cn45-1395/t.2022.03.001
0 引言
液体晃动现象普遍存在于人们的生活与生产中，液体晃动导致的安全和稳定性问题长期影响着各充液系统应用行业技术的发展。对于受外部干扰作用的贮箱內液体晃动产生的动态非平衡晃动力和晃动力矩的研究，在交通运输[1-4]、液体能源储存[5-6]和航空航天[7-13]等工程领域受到学者们的普遍关注。
建立液体晃动系统的力学模型是研究液体晃动特性的重要手段。传统单摆模型研究中，包光伟[14]针对平放式贮箱内的液体晃动特性建立单摆模型来对其进行描述;苗楠等[15]对单摆模型各个参数插值建立航天液体燃料晃动模型，并进行了变充液比工况下的输出响应仿真验证。质量-弹簧模型研究中，刘嘉一等[16]利用建立的三维质量-弹簧模型计算了水平载荷时的液体作用力;岳宝增等[17]在解析带柔性附件充液航天器耦合特性时将液体晃动等效为球摆模型。此类传统等效力学模型具有计算量小和效率高的优点，但是简化假设较多，制约了传统力学模型描述液体高阶晃动模态时的完整性，且可控、可调参数的数目较少，使其准确性也受到了限制。近年来涌现出的各类新型模型有复合模型[18]、运动脉动球模型（moving pulsating ball model，MPBM）[19]、深度学习预测模型[20]、幅度组合模型[21]和参数化模型[22-23]，以上模型对液体晃动系统的特性表达精度较传统等效力学模型有了较大提高，其中参数化模型不仅可控、可变参数多，而且描述高阶晃动模态时精准度高。文献[22-23]中的参数化模型均是在传统等效单摆模型的基础上对模型进行参数确定方法的优化，尽管比传统等效模型有所提升，但受限于传统力学模型框架结构单一的特点而无法对液体复杂工况下的晃动行为进行描述。动力学系统建模需对研究对象进行系统辨识和参数估计，测量实验和CFD模拟实验均可获得系统的输入、输出响应，但实验测量法[24-26]往往存在实验误差，且相似比选取不恰当时模拟实际工况程度较低或成本高，而CFD数值模拟方法[27-30]成本低、适用性强和准确性高，且对液体晃动系统进行激励输入时准确无延迟。
结合以上研究现状，本研究拟通过势流理论推导出部分充液圆柱贮箱液体晃动的线性化状态空间表达式，将其扩展应用到充液球形贮箱内液体晃动的动力学系统，利用CFD数值模拟对充液球形贮箱进行自振特性和不同频率下的同系统晃动模拟仿真实验，并在Matlab中利用线性灰度模型（linear grey-box models）进行参数估计，以期建立高精准度的充液球形贮箱内液体晃动的参数化模型。
1 液体晃动的动力学模型
假设圆柱贮箱内的液体为理想流体，贮箱壁面为刚体，模态坐标如图1所示。
图1中全局坐标为直角坐标系[o-xyz]，局部柱坐标[o-rθz]，以贮箱底部圆心为坐标原点，随贮箱的移动而移动，[h0]为初始液面与贮箱底部的高度差，[η（r， θ， t）]为液面波高函数。根据文献[31]中的势流理论，且不考虑环向晃动模态，当贮箱重力加速度受到[xt]的横向水平激励作用时，液体相对壁面的速度势函数和贮箱对液体的牵连速度势函数简化为：
[ϕrr， θ， z， t=n=1∞ant（cosθ）J1λnrR0coshλnzcoshλnh0]， （1）[ϕe=v（t）rcosθ]. （2）
式中：[ϕ]为液体晃动速度势函数，分为2个部分，相对速度势函数[ϕr]和牵连速度势函数[ϕe];[r]为径向坐标，取值范围是[0<r<R0];[J1]为一阶Bessel函数;[λn]为方程[J1′（λn）=0]的正根;[ant]是关于时间的任意函数;[R0]为柱体半径;[v（t）]是贮箱的横向运动速度。
波高函数为：
[ηr， θ， t=-n=1∞bnt（cosθ）J1λnrR0] . （3）
式中：[bnt]是关于时间的任意函数。
将式（1）和式（3）代入如下的自由液面处动力学等压条件和运动学等速边界条件，可得：
[∂ϕ∂t+gηSt=0] ， （4）
[∂η∂t+∂ϕ∂zSt=0] . （5）
式中：[g]为垂直于理想静止液面向下的重力加速度。为得到动力学控制方程，利用傅里叶-贝塞尔函数展开系数的正交性，综合考虑小幅晃动情况下液体晃动的一阶反对称模态占优[32]，经过代数运算，最终可将系统的状态方程转换为如下一阶微分状态方程：
[ddtb1ta1（t）=01-ω2-2ωζb1（t）a1（t）+0Aax（t）]. （6）
式中：[ω]為贮箱内液体周期晃动的一阶圆频率，[ζ]为阻尼比，[A]为输入矩阵系数。
系统输出方程为：
[FxMy=B11B12B21B22b1（t）a1（t）+C1C2ax（t）]. （7）
式中：[Fx]为x轴方向上液体对贮箱壁面的横向晃动力;[My]为y轴方向上液体对贮箱整体的晃动力矩;[B11]、[B12]、[B21]、[B22]、[C1]、[C2]均为输出矩阵系数;[ax（t）]为随时间变化的横向加速度激励。液体压强在贮箱内壁面上进行积分即可算得晃动力和晃动力矩，结合势函数和波高函数展开系数来对状态向量进行唯一确定，式（6）和式（7）构成圆柱贮箱内液体小幅晃动的线性参数化模型。在充液圆柱贮箱液体晃动建模的推导中，液体晃动的速度势函数和波高函数的傅里叶-贝塞尔级数展开系数可以构成一个完备的正交系，球形贮箱空间内切于此正交系所在空间，如图2所示。
当球形贮箱内液体小幅晃动时，其内部不产生气泡，自由液面也不出现破碎的现象，则可假设液体晃动的速度势函数和自由液面处的波高函数均满足连续无间断点条件。利用傅里叶-贝塞尔级数展开法，可将未知的球形贮箱液体晃动速度势函数和波高函数在与球形贮箱外切的圆柱贮箱区域内进行级数展开，最终也可以分别建立类似于式（6）的动力学状态方程和式（7）系统输出方程的线性参数化模型，从而实现球形贮箱内液体晃动问题的参数化建模。实际上，虽然传统的等效单摆模型与质量弹簧阻尼模型的最终动力学状态方程也具有和式（6）、式（7）相似的动力学状态方程和系统输出方程，但由于前者存在几何物理性质的限制，相应的方程形式也受到约束，实际应用中并不能实现液体晃动动态特性的完整描述。而本文提出的参数化模型却不存在该类限制，将更具有灵活性和适应性，可依据实际的晃动特性进行合理的修订和优化。这一点在有限幅非线性晃动问题上具有明显的优势，可参照非线性度相关性分析结果，对模型的输入、输出形式进行相应调整后建立适当的模型。
2 液体晃动数值模拟实验
利用CFD数值仿真软件对部分充液球形贮箱内液体小幅晃动进行模拟，设定球形贮箱的内径为300 mm，壁厚4 mm，罐内壁面弹性属性为刚体，贮箱内液体为水且体积占球形贮箱体积的一半，重力加速度与图1全局直角坐标系中[z]方向相反且大小为9.8 m/s2。
2.1 自振特性实验
工程上常用瞬态激振法进行自振特性实验，其原理是当稳定的系统结构受到瞬态激励后不再继续受到任何外力时，系统会以某一频率进行周期振动，此频率为系统的自振频率（固有频率）。
为后续建立液体晃动的动力学系统模型时对初始频率进行精准设定，寻找系统一阶固有频率十分必要。本次采用瞬态激振法对充液球形贮箱内液体晃动系统进行自振特性实验研究。在数值模拟实验中，设定外部激励为单次短时横向加速度激励，加速度幅值使液体晃动幅度在小幅晃动范围内;然后采集液体对贮箱壁面的晃动力信号和液体晃动对贮箱整体的晃动力矩信号，并对其进行傅里叶变换。晃动力和晃动力矩性质相似，仅以晃动力为例绘制晃动力信号自由变减时域图与频域图如图3、图4所示。结果表明，半充液球形贮箱内液体小幅晃动系统的一阶固有频率约为1.58 Hz，以此作为参数化等效灰度模型系统辨识时的初始频率值。
2.2 强迫晃动实验
将复杂工况下充液贮箱受到的振动简化为稳定横向简谐振动。为防止液体晃动发生共振现象而导致晃动行为超出小幅晃动范围，仿真实验中的激励频率范围设置为1.43～1.52 Hz，且加速度幅值不宜过大，实现液体晃动时径向模态为一阶占优。仿真实验完成后导出液体对贮箱的晃动力和晃动力矩信号值。以1.43 Hz激励为例将数值模拟实验数据进行时域图和频域图的绘制，图5—图7分别为加速度、晃动力和晃动力矩信号的时域图，图8为晃动力信号的频域转换且经过低通滤波滤除的频域图。
图6—图7中液体晃动的拍振现象在仿真晃动实验开始25 s左右基本消失，晃动响应逐渐由拍振过渡到标准的简谐稳态振动，基频约为加速度激励的频率值1.43 Hz。由图8可知，当横向加速度外激励频率为1.43 Hz时，充液贮箱液体晃动系统的晃动力输出信号的一阶固有频率与外激励频率同步。
3 参数估计和系统辨识
数值模拟仿真实验中的液体晃动系统输入信号和输出信号已知，而部分系统参数未知，此类系统辨识问题属于灰度模型参数辨识问题。
3.1 参数估计
依据推导得到的状态空间方程式（6）和输出方程式（7），在Matlab中建立含有未知参数的初始动力学参数化模型，基于经验设定初始化的频率、阻尼比、输入系数、晃动力输出系数、晃动力矩输出系数、直接输出系数，导入模拟实验数据后对灰度模型进行参数估计，系统辨识的流程如图9所示。
当贮液容器内充液比固定时，该系统的特征属性唯一确定，虽然不同激励下的系统输出响应不同，但是输入信号与输出信号的转换关系必然相同，因此，进行同一系统的不同频率晃动模拟实验。以1.43 Hz和1.52 Hz加速度激励时数值仿真实验与辨识所得模型的输出响应结果的对比为例，如图10—图13所示。根据2-范数计算曲线间的距离并换算辨识度后可得：激励频率为1.43 Hz时，参数化模型对比数值仿真实验的晃动力和晃动力矩输出响应的辨识度为96%以上;激励频率为1.52 Hz时，辨识度为92%以上。
3.2 系统辨识结果
在同一激励频率下进行多次数值模拟实验。为了使系统参数尽可能接近真实值，系统辨识后需要对所得同类参数进行最小二乘法的统计学处理，不同激励频率处理后的结果如表1所示。由表1可知，不同激励频率下系统的一阶固有频率[（f）]辨识為1.57 Hz左右，与自振特性实验所得基频接近程度达到98%以上，阻尼比[ζ=0.006～0.010]，输入、输出系数较为稳定且可信度高，本参数化模型对数值模拟仿真实验晃动力和晃动力矩的输出响应的辨识度均在92%以上。
4 结论
本文扩展充液圆柱贮箱求解区域后建立含有未知参数的充液球形贮箱晃动系统参数化动力学模型，通过CFD仿真软件对充液球形贮箱进行自振特性实验和不同激励频率下的强迫晃动实验，依据液体压强等式和数值仿真实验数据在Matlab中对系统状态向量和未知参数进行计算和辨识，研究参数化模型参数估计结果并对比分析数值模拟实验与参数化模型的输出响应，得出以下结论：
1）依据圆柱贮箱横向激励速度势函数和波高函数中傅里叶-贝塞尔函数的正交性且圆柱空间与球形空间相切的特性，得到横向激励下充液球形贮箱系统的状态空间表达式具有可行性，可将此原理类推到与圆柱贮箱相切的任意旋转对称贮箱。
2）自振特性实验频域分析能够精准计算分析出系统的基频，是确定系统辨识频率初值的可靠 方法。
3）小幅简谐横向加速度激励下的液体受迫晃动的数值仿真实验结果表明：实验前期液体晃动的“拍振”现象较为显著，原因在于紧靠贮箱内壁面的液体和远离壁面的液体晃动相位不同步;之后液体晃动的“拍振”现象逐渐减小直到消失，这是由于在外部激励的持续作用下，液体晃动系统的各个子系统的相位均同步于外激励相位。
4）本研究对不同频率横向简谐加速度频率下的充液球形贮箱内液体晃动系统进行辨识后，所得的参数化模型中参数值之间略有差别，系统基频的辨识极佳且晃动力和晃动力矩输出响应的辨识准确度均在92%以上，总体实现了充液球形贮箱液体晃动系统的高精度参数化等效建模的目标。
参考文献
[1] 茅海剑， 李波， 贝绍轶， 等. 基于流固耦合的液罐车侧倾稳定性研究[J].江苏理工学院学报，2021，27（4）：56-67.
[2] 万滢. 液罐车辆车-液耦合动力学特性与防侧翻控制方法研究[D].长春：吉林大学，2018.
[3] 黄中烈. 带贮罐类重型车辆动力学建模与刚-液耦合特性研究[D].柳州：广西科技大学，2019.
[4] 李超，吴文军，周凤霞，等. 重型罐车内液体晃动参数化数学模型辨识[J].广西科技大学学报，2021，32（2）：65-71.
[5] 邓传奇. 基于ANSYS的船用双层C型LNG储罐结构设计研究[D].淮南：安徽理工大学，2018.
[6] 赵鹤锦. 大型储液罐在地震作用下的响应分析[D].大连：大连理工大学，2021.
[7] 于瑾，李宝海，王求生，等.贮箱防晃挡板结构形式和晃动阻尼研究[J].导弹与航天运载技术，2021（5）：147-149，154.
[8] 刘峰， 岳宝增， 唐勇. 多充液贮腔航天器耦合动力学与姿态控制[J].宇航学报，2020，41（1）：19-26.
[9] 岳宝增，于嘉瑞，吴文军. 多储液腔航天器刚液耦合动力学与复合控制[J].力学学报，2017，49（2）：390-396.
[10] 孙梓煜， 岳宝增， 刘峰，等. 充液柔性航天器刚-液-柔耦合动力学研究的凯恩方法[J]. 宇航学报，2021， 42（5）：552-561.
[11] 张海涛， 孙蓓蓓. 匀加速激励下的液体晃动力解析计算[J].振动与冲击，2020，39（12）：158-163.
[12] 杨旦旦， 岳宝增. 低重环境下旋转轴对称贮箱内液体晃动研究[J].宇航学报，2013，34（7）：917-925.
[13] 陈青全. 卫星贮箱液体晃动等效力学模型参数辨识[D].长沙：国防科学技术大学，2014.
[14] 包光伟. 充液卫星平放式贮箱内液体晃动的等效力学模型[J].宇航学报，1996，17（1）：67-70，106.
[15] 苗楠，刘战合，王晓璐，等. 考虑燃料消耗的飞机副油箱燃油晃动等效建模[J].科学技术与工程，2020，20（29）：11851-11857.
[16] 刘嘉一， 白欣然， 李晓轩. 质量-弹簧模型在储液容器抗震分析中的应用[J].原子能科学技术，2013，47（6）：947-951.
[17] 岳宝增，祝乐梅. 携带晃动燃料柔性航天器姿态机动中的同宿环分叉研究[J].宇航学报，2011，32（5）：991-997.
[18] 刘峰， 岳宝增， 马伯乐， 等. 燃料消耗下充液航天器等效动力学建模与分析[J].力学学报，2020，52（5）：1454-1464.
[19] 邓明乐， 岳宝增， 黄华. 液体大幅晃动类等效力学模型研究[J].宇航学报，2016，37（6）：631-638.
[20] 金晓威， 赖马树金， 李惠. 物理增强的流场深度学习建模与模拟方法[J]. 力学学报，2021，53（10）：2616-2629.
[21] 苗楠， 李俊峰， 王天舒. 横向激励下液体大幅晃动建模分析[J].宇航学报，2016，37（3）：268-274.
[22] 李超， 吴文军， 周凤霞， 等. 重型罐车内液体晃动参数化数学模型辨识[J].广西科技大学学报，2021，32（2）：65-71.
[23] 周凤霞， 吴文军， 王佐， 等. 带贮罐类重型车辆建模与操纵稳定性研究[J].广西科技大学学报，2020，31（2）：68-79.
[24] 王坤. 水平激励下矩形贮箱内液体晃动力和倾覆力矩实验研究[D].泰安：山东农业大学，2018.
[25] 吴文军， 高超南， 岳宝增， 等. 圆柱贮箱内液体非线性稳态晃动实验及动力学特性分析[J].宇航学报，2021，42（9）：1078-1089.
[26] 王琼瑶. 部分充液罐车罐体内液体晃动的动力学特性研究[D].广州：华南理工大学，2017.
[27] WANG W Y， GUO Z J， PENG Y， et al. A numerical study of the effects of the T-shaped baffles on liquid sloshing in horizontal elliptical tanks[J]. Ocean Engineering， 2016， 111：543-568.
[28] 李骏，龚思惠，张鹏飞，等. 基于FLUENT的货车燃油箱液体晃动数值模拟[J].科学技术与工程，2020，20（34）：14152-14156.
[29] 张展博，李胜强. 圆柱水箱中水平多孔挡板对液面晃动影响的数值模拟研究[J].清华大學学报（自然科学版），2018，58（10）：934-940.
[30] 刘汉武， 张华， 胡震宇， 等. 基于SPH方法航天器含贮箱液体晃动分离动力学研究[J].中国科学：技术科学，2021，51（8）：938-947.
[31] 吴文军. 带多充液贮箱及机动柔性附件航天器耦合动力学研究[D].北京：北京理工大学，2015.
[32] 李遇春. 液体晃动动力学基础[M].北京：科学出版社，2017.
System identification and parametric equivalent modeling of liquid lateral sloshing in spherical tanks
WANG Zuo， WU Wenjun*， JIANG Bingchao， GAO Chaonan
（School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology，
Liuzhou 545616， China）
Abstract：
The study aims at inaccuracy and incompleteness of the traditional equivalent mechanical model in describing the sloshing characteristics of the liquid in the tank. Firstly， we derive the dynamic state equation of partially liquid-filled cylindrical tank by using potential flow theory and expand the solution area of liquid sloshing in the spherical tank to the solution area of the cylindrical tank tangent to the outside based on the orthogonality condition of the Fourier-Bessel series expansion coefficient in the liquid sloshing velocity potential function and the sloshing wave height function in the equation， thus establishing a parameterized model to describe the liquid sloshing in the spherical tank. Then， we obtain the liquid sloshing force and torque of the spherical tank under lateral harmonic excitation with small sloshing through the CFD simulation experiment. Finally， we use the gray box identification function of Matlab to identify the undetermined coefficients in the established parameterized model， realizing the high-precision equivalent modeling of the liquid lateral sloshing in the spherical tank. The results show that the parameterized equivalent mechanical model can describe the natural frequency， sloshing force and torque of the liquid in the spherical tank more accurately and completely， and the identification accuracy can reach over 92%， which can give reference for the online real-time identification and practical engineering.
Key words：
spherical tank; liquid sloshing; parametric modeling; CFD simulation; gray box identification; system identification
（责任编辑：罗小芬）