韩越洋， 朱翔,3， 郭文杰， 李天匀,3

(1.华中科技大学 船舶与海洋工程学院，湖北 武汉 430074；2.船舶与海洋水动力湖北省重点实验室，湖北 武汉 430074；3.高新船舶与深海开发装备协同创新中心， 上海 200240；4.华东交通大学 铁路环境振动与噪声教育部工程研究中心，江西 南昌 330013)

圆柱壳作为一种基本结构单元，被广泛应用于液体储存及运输方面。对于许多液体散装运输设备，如远洋运输船、公路铁路油罐车和太空运输飞船等，其方向稳定性、安全性以及结构的可靠性很容易受到内部部分填充的液体的影响。在接近共振的情况下，由液体晃荡产生的局部力和局部力矩可能会导致结构的不稳定，甚至圆柱壳结构和内部流体相互作用影响到结构设备性能，可能成为一些事故的原因[1-4]。学者通过数值方法等对壳体中的液体运动进行了研究[5-10]。其中Wang[5]利用几何比例边界有限元法，研究了水平圆柱体内部的液体晃动。Liu等[6]采用运动网格结合流体体积方法来捕捉晃动过程中液态氢罐自由表面的波动，研究了在正弦激励下晃荡的液体的水动力性能。文献[11-23]利用解析法、半解析法研究壳体中液体的晃荡。Evens等[11]利用伽辽金法构造矩阵求解特征值，对刚性边界下水平半充液圆柱壳的内部液体的晃荡频率进行了计算。McIver[12]采用二维双极坐标系并且将特征值问题转化积分方程求解问题，计算了水平二维充液圆柱壳的液体晃荡频率。Hasheminejad等[13]利用连续的共形坐标变换提出了精确的二维流体动力学分析，分析了内部有挡板的水平部分充液圆柱壳的液体横向晃荡特性。Papaspyrou[14-15]根据Evans的半解析法，先后分析了横向、纵向外部激励下，水平半充液圆柱壳的晃荡响应。文献[24-26]采用实验手段对水平充液壳体中液体的晃荡运动进行了研究。

针对二维充液壳体中液体晃荡的解析解或者半解析解大多只针对半充液圆柱壳或者矩形截面。由于建模和求解上的难度，针对充液深度变化的圆柱壳截面自由液面晃荡的解析方法，涉及文献较少。本文针对水平任意充液深度的二维圆柱壳中液体晃荡运动求解，提出一种半解析法。首先在圆柱壳和自由液面上分别建立2个极坐标系，根据2个坐标系的变换关系将圆柱壳和自由液面的2个坐标系联系起来。根据拉普拉斯方程得出自由液面的速度势函数，进而根据边界条件及坐标关系，得出自由液面的运动控制方程，并通过伽辽金法解出自由液面晃荡频率。通过改变充液深度，对二维圆柱壳自由液面晃荡频率的特性进行了分析。

1 二维任意深度充液刚性圆柱壳晃荡理论模型

如图1所示，二维圆柱壳上建立以壳体圆心O为坐标原点的极坐标系(r,φ)，且φ取值范围为[-α,α]。自由液面上则建立以液面中心Q为坐标原点的极坐标系(R,θ)，θ取值范围为[-π/2,π/2]。H为自由液面坐标系到圆柱壳结构坐标系的距离，H值为正则液面高度小于半径；H值为负则液面高度大于半径。2个坐标系在空间上任意一点的夹角定义为β=φ-θ。当H值为正时，∣β∣=θ-φ；当H值为负时，∣β∣=φ-θ。rs为圆柱壳半径。H值与充液深度h的关系为h=rs-H。

图1 二维刚性充液圆柱壳坐标系示意

假设液体为无粘性且不可压缩的理想流体，则对于小幅晃荡的理想流体，其速度势函数应满足拉普拉斯方程：

(1)

式中φ(r,θ)为液体速度势函数。

由于液面为小幅晃荡,还需要满足一维线性波的表面边界条件:

(2)

式中：θ=±π/2 ；K=ω2/g；ω为液体晃荡的圆频率。

假设壳体为刚性壳体，则速度势函数φ(r，θ)应满足流体与刚性壁面的边界条件，即：

(3)

其中r=rs，rs为圆柱壳半径。

通过上述约束条件，自由液面晃荡的速度势函数的解可利用分离变量法写为对称/反对称2种形式的解。首先，速度势函数满足拉普拉斯方程(1)，则速度势函数反对称与对称模态的解分别为：

(4)

(5)

对于自由液面反对称模态形式，为后续方便计算可将式(4)改写为：

(6)

将式(6)代入式(2)中，可得：

(7)

从式(7)中可得A2n-1与A2n之间的关系为A2n=-K/(2n)A2n-1，代入式(6)中，可得：

(8)

将式(3)从结构坐标系转换至液面坐标系上，即：

(9)

由图1可知结构坐标系和液面坐标系的关系：

R2=H2+r2-2Hrcosφ

(10)

Rsinθ=rsinφ

(11)

则可得液面坐标系中坐标R和θ对结构坐标系r求导关系式：

(12)

(13)

将式(8)、(12)以及(13)代入式(9)可将边界条件为：

2)θ+φ]-KR2n-1sin[(2n-1)θ+φ]]}=0

(14)

(D-KT){A2n-1}=0

(15)

晃荡频率的求解可利用齐次线性方程组系数矩阵行列式为零进行求解，即：

|D-KT|=0

(16)

利用搜根法即可解出晃荡频率。同理，对于自由液面对称模态，计算方法与反对称模态计算方法类似，首先将速度势函数改写：

(17)

速度势函数需要满足约束条件式(2)，可得；

(18)

则速度势函数可写为：

(19)

将式(12)、(13)以及(19)代入式(9)中可得：

KR2ncos[2nθ+φ]]}=0

(20)

(L-KY){B2n}=0

(21)

其中矩阵L=[lmn]和矩阵Y=[ymn]的元素分别为：

此时晃荡频率也可以通过令式(21)中系数矩阵行列式为零求解得到。

2 收敛性分析及正确性验证

通过本文的理论推导即可以得出自由液面晃荡频率。由于势函数是通过无穷级数表达，因此有必要针对级数的截断项数进行收敛性分析。给定圆柱壳半径rs=1 m, 液面高度h为0.5 m，分别计算m，n的截断数N=5,10,20,30,40时，壳体自由液面无量纲晃荡频率Ω=Krs，如表1所示。很明显在N=20时晃荡频率已经稳定。从表中得出N=40时可以得到更精确的值，但计算时间也会大大延长。所以本文采用截断数N=20，既能保证精度，又可以节省时间。

表1 二维刚性充液圆柱壳自由液面的无量纲晃荡频率Ω收敛性分析

为了验证本文理论模型和算法的正确性，计算二维刚性充液圆柱壳不同充液深度下自由液面的晃荡频率，并同文献结果进行对比分析。计算参数如下：圆柱壳半径为rs=1 m，H分别为0.5,0,-0.5 m，对应充液深度h分别为0.5,1,1.5 m。由方程(15)和(21)可以得出不同充液深度下二维刚性充液圆柱壳自由液面的晃荡频率，如表2所示。其中q为模态阶数。

表2 不同充液深度下二维充液圆柱壳自由液面的无量纲Ω晃荡频率对比

表3 二维刚性充液圆柱壳自由液面的晃荡频率f有限元和本文方法结果对比

进一步将本文计算的晃荡振型和有限元结果进行对比分析。图2为H=0.5 m时流体有限元单元模型。图3为本文方法和有限元法计算得到的自由液面前4阶模态振型对比，其中有限元计算出的振型根据模型中自由液面节点数据绘制得出。从图3可见本文计算得到的自由液面模态振型和有限元计算的模态振型也吻合很好。由此可见，和已有文献晃荡频率以及和有限元法晃荡频率和振型对比，都表明本文理论方法的准确性及可靠性。

图2 二维刚性充液圆柱壳中流体有限元网格模型

图3 二维刚性充液圆柱壳自由液面的晃荡振型

3 不同充液深度下自由液面晃荡频率分析

对壳体半径为rs=1 m，不同充液深度下(H=0.7 m，0.5 m，0 m，-0.5 m，-0.7 m)二维刚性充液圆柱壳自由液面晃荡频率进行计算并对比分析。

图4(a)给出了不同充液深度下晃荡频率随模态阶数q的变化。图4(b)则给出了晃荡频率随充液深度的变化曲线。由图4(a)可知，充液深度一定时，充液圆柱壳自由液面晃荡频率随着模态阶数q的增大而增加。根据文献[27]提出的晃动流体等效力学模型分析，随着模态阶数增加，内部液体晃荡的固定质量(即不参与晃荡的流体质量)占比逐渐增大，而晃荡质量占比则逐渐减小，晃荡固有频率逐渐增大。对于不同的充液深度，以上结论也都成立。

图4 二维刚性充液圆柱壳自由液面晃荡频率

结合图4可知，在第1阶模态的时候，晃荡频率随着深度的增大而增大。这是由于在1阶模态时，随着液深的增加，内部流体的固定质量占比逐渐增大，晃荡程度减小，晃荡质量占比逐渐降低，因此晃荡频率呈增大的规律。

但是在较高阶模态时，如第5阶模态时，晃荡频率则是随着充液深度的增大先减小后增大，且在半充液的时候出现转折。这可以解释为在高阶晃荡模态中，液面晃荡的节点增多，而在圆柱壳中液深关于壳体圆心上下的变化会导致液体自由表面区域的变小，相对刚度增大，使得激发高阶晃荡难度变大，则晃荡频率增大。

并且，在H=0.7 m时的晃荡频率分别在第2阶模态时超过H=0.5 m时的晃荡频率，第3阶模态的时候超过H=0 m的晃荡频率，第5阶模态的时候超过H=-0.5 m的晃荡频率。综上得出充液深度越小，晃荡频率随着模态阶数递增的速率也越快，图4(a)中曲线的斜率中也可直观反映出这一结论。

4 结论

1)特定深度下，自由液面晃荡频率随着模态阶数的增大而增大；在模态阶数q=1时，自由液面晃荡随充液深度的增大而增大；在高阶模态中，自由液面晃荡频率随着充液深度的增大先减小后增大，且在半充液状态时发生转折；不同充液深度下，自由液面晃荡频率增大速率随着充液深度的减小而增大。

2)本文的半解析法与有限元法相比建模简便，参数分析时无需重复建模，计算结果准确、可靠。与其他解析方法相比，具有良好的准确性和适应性，能适应不同充液深度下圆柱壳的晃荡频率求解，并且能进一步拓展到三维圆柱壳中的自由液面晃荡问题以及弹性圆柱壳的流固耦合振动分析。