|
摘要：在理想流体及小幅度晃动的假设条件下，表征液体晃动的压力函数的控制方程可以线性化，经有限元方法离散后，便能导出结构振动与液体晃动耦合问题的二阶线性微分方程组，其中系数矩阵是非对称的。采用静态缩聚技术，方程中的液体非自由面上的节点压力可被凝聚掉。进一步用模态坐标表示液体自由面的压力向量，得到了具有对称系数矩阵的耦合方程组。最终，液体是以超单元形式出现在方程中，从而可以用NASTRAN 程序求解结构与液体的耦合动力学问题。
关键词：流固耦合；模态分析；有限元
1 引言
结构振动与液体晃动是航天器结构动力学的一个重要问题。分析结构与液体的耦合模态特性不外乎解析和数值两种方法[1][2]。解析方法适应于构型简单的结构，其中NASA 技术报告[3][4]是处理液体运载火箭贮箱内液体的经典建模方法。当结构模型较为复杂时，一般只能用有限元方法进行液固晃动耦合模态分析。
考虑到工程应用的实际需要，在用有限元方法分析结构-液体耦合模态时，不仅要求理论上是正确的，还要求方法能够在通用程序的框架下实现，如通用结构分析程序NASTRAN。满足这样要求的方法才便于应用、更具有可靠性，适合工程实际使用。
对于小幅度晃动的结构-液体耦合问题，表征液体运动的方程可以线性化，王文亮在其著作文献中已给出了不同形式的控制方程[5]，其中用液体晃动主模态表示的液固耦合控制方程具有对称性。在这组方程中，液体如同附加在结构上的超单元(也可以认为是子结构)，因此能将液体表示为超单元，在NASTRAN 程序中实施结构-液体的耦合分析。
2 小幅晃动条件下结构—液体耦合方程的基本形式
假设液体是无粘、无旋的不可压缩理想流体，液体所占区域Ω 内的压力p 满足Laplace方程 （图片） 在边界条件(2)和(3)中，n 表示液体表面的法向单位矢量(指向液体外部)， nu 为结构沿n 方向的加速度，ρ为液体密度， g 为重力加速度。
用有限元方法对液体压力进行插值，将压力函数p 表示成（图片） 其中p 是液体有限元模型的节点压力向量， N 为插值型函数。根据插值公式(4)式可以推导出液体压力∇p梯度的表达式（图片） 于是，方程(1)和(2),(3)边界条件可以写成离散形式（图片） 其中（图片） 结构有限元模型的控制方程为（图片） 式中u 为结构的节点位移向量， S M 和S K 分别是结构的质量和刚度矩阵， S Q 为结构除液体压力之外所受的外力向量。
将(5)和(7)写成联立形式，就得到结构振动与液体晃动的耦合方程（图片） 3 液体节点压力自由度缩聚
将液体的节点压力向量划分成自由面压力f p 和非自由面压力r p 两部分，即（图片），液体质量矩阵 MF 和刚度矩阵KF 也相应地写成分块形式。从公式(6)可知，计算矩阵 MF 的面积分仅限于自由面，因此矩阵 MF 只有与自由面相关的元素非零，因此（图片） 将分块形式的 MF 和 KF 代入方程(5)，即有（图片） 其中L=[Lf Lr]。从方程(9)的第二式可以确定 pr 与 pf 约束关系（图片） 将上式代回到从方程(9)的第一式，得到仅由自由面节点压力表示的液体控制方程（图片） 其中（图片） 再将(10)代入方程(7)，则可得到结构的控制方程（图片） 其中（图片） 将(11)和(12)写成联立形式后得到（图片） 4 液体压力主坐标表示的耦合方程
求解广义特征值问题（图片） （图片） 实际上广义特征值问题(14)就是当结构为刚性时液体自由晃动的特征方程，特征值ω1,...ωnf 就是刚壁容器内液体的自由晃动频率，且ω1=0。主座标ξ1表示液体自由面的均匀压力，即自由面平直升降。特征向量满足正交性条件（图片） 其中（图片） 将主坐标表示的液体自由面压力向量(15)代入到方程(11)中，利用正交性条件可以得到（图片） 式中ξE= {ξ2....ξnf}T。
在零初始条件下，从方程(17)的第一式可解出（图片） 再将(18)和(19)代入到方程(12)中，又可得到（图片） 其中（图片） 用（图片）同乘以方程(17)第二式两端，并于方程(20)联立后最终得到（图片） 这是一组对称形式的方程组，便于用常规多自由度动力学方法求解。若将液体看作由界面坐标u 和内部坐标 ξE 表示的超单元，其质量和刚度矩阵分别是（图片） 将此超单元综合到结构上，耦合以后方程就是(21)。于是，只要构造出矩阵μ 和κ ，就能用MSC.Nastran 程序进行后续的液固耦合动力学分析。
5 算例与结果
(a) 圆柱贮液器液固耦合模态
圆柱贮液器模态分析算例的参数取自文献[6]，模型如图所示。贮液器腔体为铝合金材质，弹性模量E= 60GPa，泊桑比μ=0.35，密度 ρ= 2000kg/m3 ，高度L=625mm，直径D=200mm，厚度1mm。底板和加筋条的材质与腔体相同，厚度均为3mm、宽度为28mm。贮液器注水高度0.7L = 438mm。
生成矩阵μ 和κ 以后，用Nastran 的DMIG 格式写入文件Fluid.pch，而圆柱贮液器结构模型用Nastran 常规的bdf 格式文件Structure.bdf 表示。为了计算液体晃动与结构振动的耦合模态，写出Nastran 的模态分析求解代码如下。（图片）
图1 圆柱贮液器结构示意图 SOL 103
CEND
ECHO = NONE
SUBCASE 1
METHOD = 1
K2GG = KAAX
M2GG = MAAX
VECTOR(SORT1,REAL)=ALL
BEGIN BULK
PARAM POST 0
PARAM AUTOSPC YES
EIGRL 1 256 0
include 'Structure.bdf'
include 'Fluid.pch'
ENDDATA
运行Nastran 以后即可得到耦合模态和频率。部分结果如表1。计算和实测结果基本吻合。表1 圆柱贮液器液固耦合模态的计算式实测结果
（图片）(b)带端壳圆柱贮箱的液固耦合模态
某贮箱结构如图2 所示，贮箱中段为圆柱壳，两端是椭球壳。贮箱内部分充液。采用算例1 相同的方法，在生成液体质量和刚度矩阵以后，用Nastran 分析耦合模态。部分频率计算结果如表2。（图片） 根据形状，模态可以分为四种类型。第1 种如图3，结构呈刚体运动(上下平动、左右平动，以及绕z 轴转动)，液体表面几乎没有晃动，频率近似为零；第2 种如图4，结构呈刚体运动(绕x, y轴摆动)，液体随结构晃动；第3 组如图5，结构几乎不动，液体在“固壁容器”内晃动；第4 组如图6~8，结构与液体存在比较明显的耦合振动。（图片） 在第四组模态中，液体表面周向波数为0 属于液固耦合纵向模态，波数为1 属于横向模态，其余是局部模态。在液体运载火箭的动力学分析中，重点需要关注周向波数为0 和1 的模态。
6 结论
将液体视为超单元，就可以利用通过结构分析程序对充液贮箱进行液固耦合模态分析，极大地简化了此类问题的分析。计算与实测结果基本吻合，验证了方法的正确性。文中的算例虽然只给出了模态分析，但是采用同样的建模方法，也能够进行频率响应和瞬态响应分析，并且也适用于多个贮箱问题。
参考文献
[1] Miao Guo-ping, Analytical Solutions for the Sloshing loading on Circular Cylindrical Liquid Tanks with Interior Semi-porous Barriers, Journal of Hydrodynamics , Ser. B ,2(2001), 32-39.
[2] Alfredo Bermudez, et. al., Finite element computation of sloshing modes in containers with elastic baffle plates, Int. J. Numer. Meth. Engng., 2003; 56:447–467.
[3] G.B. Paddock, Dynamic Stability of Space Vehicles Volume I – Lateral Vibration Modes, NASA CR-935.
[4] J.A. Staley, Dynamic Stability of Space Vehicles Volume II – Determination of Longitudinal Vibration Modes, NASA CR-936.
[5] 王文亮等，结构动力学与动态子结构方法：复旦大学出版社，第一版(1985)。
[6] 万水等，圆柱贮液器固液耦合模态分析，工程力学，第17 卷第3 期(2000)，87~92。
5/22/2011
|
