带填充结构圆柱形容器内液体横向晃荡特征研究

王肖飞，刘 俊，叶文斌

(1.大连理工大学 海岸与近海工程国家重点试验室， 辽宁 大连 116024；2.大连理工大学 工程抗震研究所 建设工程学部水利工程学院， 辽宁 大连 116024)

液体晃荡，广泛存在于运输、航天、土木、核电领域等诸多实际工程中，也是储液罐设计主要考虑因素之一。对液体晃荡的研究由来已久，基于小波振幅的假设，1954年Senda等[1]率先求得了圆柱形容器液体晃荡的解析解。近几年，圆柱形容器液体晃荡在实验及数值研究方面取得了很多的成果，并且FLURNT等流体力学仿真软件也得到了广泛的应用[2-3]。忽略流体黏性并假设流体不可压缩，Tetsuya等[4]研究了地震激励下带浮顶圆柱形储桶的晃荡响应。Sygulski[5]在不可压缩势流假设的基础上，采用边界元方法研究了三维液体晃荡问题。采用变分原理，Takahara等[6]分析了在俯仰激励下圆柱形容器非线性晃荡特性，并通过实验验证了非线性分析结果的有效性。Lee等[7]通过假设流动为小振幅理想晃动，分析了底部有间隙的环形圆槽内液体晃荡的固有频率。Zhen等[8]应用变量分离和拉普拉斯变换方法研究了圆柱形储罐内分层液体的三维晃荡运动。Chen等[9]采用正则化边界积分法研究了二维和三维容器(包括圆柱形容器)中的非线性液体晃荡问题。

比例边界有限元方法(Scaled Boundary Finite Element Method，SBFEM)，由Wolf等[10]在20世纪90年代提出，该方法结合了有限元方法和边界元方法的优点，仅需离散计算域边界，使计算维度减一。相对于边界元方法，该方法无需基本解；相对于有限元方法，它具有在比例坐标下径向解析的特性，从而提高了数值计算精度。许多学者已经成功将该方法应用于液体晃荡领域。例如，Lin等[11]将比例边界有限元方法应用在水平外部激励下圆柱形或者轴对称水箱液体晃荡的分析。Liu等[12]基于线性势流理论，研究了带T形挡板的卧式圆柱罐内液体晃荡问题。应用比例边界有限元方法，Wang[13]研究了二维椭圆形容器内T形挡板对液体晃荡响应的影响。本文基于线性势流理论，应用变分原理推导出带填充结构圆柱形容器内SBFEM控制方程并引入相应边界条件进行求解。同时开展了在横向外部激励作用下相对波数、填充介质系数和填充半径等参数对液体晃荡特征影响的研究。

1 比例边界有限元方程推导

如图1所示，假设液体无黏且不可压缩，圆柱形容器半径为b，同轴圆柱形填充结构半径为a，容器内液体深度为H。假设容器底部固定，容器壁为刚性。在x轴方向上作用外部水平位移激励x=Ae-iωt，其中A为位移的幅值，ω为角频率，t为时间，同时假设A足够小，液体自由表面满足线性势流理论。

图1带填充结构圆柱形容器示意图

将整个计算域划分为填充子域Ω1和非填充子域Ω2。其中，填充结构内可采用“固体泡沫球”[14]等填充物，使得该填充结构满足Sollitt等[15]在1972年所建立的流体与大孔隙多孔介质作用数学模型。

基于线性势流理论，分离出时间因子，将速度势表示为

Φ(x,y,z,t)=φ(x,y,z)e-iωt

(1)

复速度势函数φ(x,y,z)在整个计算域内满足拉普拉斯方程

(2)

及线性自由液面条件

(3)

容器边界条件

(4)

其中g为重力加速度。

引入边界条件，将变量z分离出来，速度势函数φ(x,y,z)可表示为

(5)

上式右边第一项代表传播模态，第二项代表非传播模态。k0和km(m=1,2,…,∞)分别为传播模态和非传播模态的波数，且满足下述色散方程

ω2=gk0tanhk0H

(6)

ω2=-gkmtankmH，m=1,2,…,∞

(7)

计算域中φ0(x,y)和φm(x,y)分别满足二维Helmholtz方程及修正的Helmholtz方程

(8)

(9)

用φⅠ和φⅡ分别代表计算子域Ω1、Ω2的速度总势，在r=a处，两个子域之间的耦合边界条件可表示为

εφⅠ,n=-φⅡ,n

(10)

(s+if)φⅠ=φⅡ

(11)

其中：ε为多孔介质孔隙度系数；f为线性阻力系数；s为惯性系数。

容器边壁处边界条件表示为

φⅡ,n=-iωAcosθ，r=b

(12)

其中：θ为域内任一点与 轴正向所成的角度。

采用变分原理推导比例边界有限元方程。φ0(x,y)和φm(x,y)在边界上的法向导数分别满足

(13)

(14)

方程(8)、式(9)的等效泛函形式可以表示为

(15)

(16)

图2比例边界有限元计算方法示意图

(17)

(18)

在比例坐标系统中拉普拉斯算子表示为：

(19)

根据等参变换概念，速度势函数可采用插值函数N(s)进行离散：

φ(ξ,s)=N(s)φ(ξ)=φT(ξ)NT(s)

(20)

将式(19)、式(20)代入式(15)，可得传播模态SBFEM控制方程及边界条件：

(21)

(22)

(23)

式中系数矩阵

(24)

E1=0

(25)

(26)

同理得非传播模态SBFEM控制方程及边界条件

(27)

(28)

(29)

2 比例边界有限元方程求解

域Ω1内速度势可以重新表示为：

(30)

式(21)、式(27)为齐次贝塞尔微分方程，其解的形式可表达为：

(31)

(32)

其中：Tj为n阶向量；c0j和cmj为系数；Jrj(ξ)为第一类贝塞尔函数；Irj(ξ)为第一类修正的贝塞尔函数。考虑贝塞尔函数和第一类修正贝塞尔函数的性质并将式(31)代入式(21)，将式(32)代入式(27)，可得

(33)

(34)

同时考虑c0jJrj(ξ)、cmjIrj(ξ)的任意性可得特征方程

(35)

类似地，Ω2内速度势可以表示为

(36)

其中φⅡ0(ξ)和φⅡm(ξ)分别满足比例边界有限元控制方程式(21)和式(27)，以及边界条件式(22)、式(23)、式(28)和式(29)。方程式(21)和式(27)的解的形式可表达为：

=T(J(ξ)Α0+Y(ξ)B0)

(37)

=T(I(ξ)Αm+K(ξ)Bm)

(38)

其中：a0j、b0j、amj和bmj为系数；Yrj(ξ)为第二类贝塞尔函数；Krj(ξ)为第二类修正的贝塞尔函数。

(39)

(40)

当r=a，由式(9)得：

(41)

将式(31)、式(32)代入式(23)，将式(37)、式(38)代入式(28)，联立后式(11)、式(41)表示为：

(s+if)J0aC0=J0aA0+Y0aB0

(42)

(s+if)ImaCm=ImaAm+KmaBm

(43)

(44)

(45)

联立式(31)、式(32)、式(37)、式(38)、式(39)、式(40)、式(42)、式(43)、式(44)、式(45)，可求得所有系数矩阵，代入式(30)、式(36)可求得Ω1、Ω2的速度势。

水平向晃荡力Fx以通过对z进行积分得到

(46)

3 计算结果及分析

为了验证SBFEM方法在填充结构圆柱形容器内液体晃荡问题的精确性和高效性，将本方法计算结果与解析解[16]进行对比。选取计算参数如下：填充结构半径a=0(相当于无填充结构)，容器半径b=1.0 m，重力加速度g=9.81 m/s2，液体密度ρ=1.0 g/cm3，容器内液面深度H=2.0 m，A=0.1 m，波数k=0.6。采用4个、12个、20个二阶3节点单元对容器边壁有限元离散。图3给出了SBFEM计算结果与解析解归一化后容器边壁处自由液面液面高程(η/A)的对比。从图3中可以看出，在较少的离散单元下，即可得到与解析解相吻合的计算结果，从而验证了本方法具有很高的精度及计算效率。

图3容器边壁液面高程与解析解的对比

用本文提出的方法，对四组不同填充结构系数情况下容器内晃荡问题进行了研究。系数ε、f取值分别为：(1)ε=0.20、f=5.0；(2)ε=0.45、f=2.0；(3)ε=0.80、f=2.0；(4)ε=1.00、f=0.0。其中惯性系数取定值s=1.0。对于流体而言，即流域内没有填充物时，孔隙率系数ε=1.0，阻力系数f=0，惯性力系数s=1.0。在此基础上研究了相对波数kb、填充结构半径a、孔隙度系数ε和线性阻力系数f等参数对容器内晃荡力及自由液面高程的影响，以下为数值计算的结果。

图4描述了填充半径a=0.3、a=0.5时相对波数kb变化对系统晃荡力Fx的影响。从图中可以看出，随着相对波数的增加，Fx呈现出先增大后减小之后再缓慢增大的趋势；四组参数中，随着ε的增加及f的减小，Fx峰值先减小后增大，峰值对应的相对波数随之增大；且当ε=0.45、f=2.0时所对应的Fx的峰值最小。

图4系统晃荡力随相对波数的变化

图5描述了当填充系数为ε=0.45、f=2.0时，考虑四种填充半径a(取值分别为0.1、0.3、0.5、0.7)条件下系统晃荡力Fx及x正方向上容器边壁处液面高程η/A随相对波数的变化。由图中可以看出，填充半径a增大时，Fx峰值及η/A峰值均减小，Fx及η/A峰值对应相对波数也随之减小；当填充半径a由0.1增加0.5时，Fx峰值及η/A峰值减小幅度较大，当填充半径a由0.5增加到0.7时，Fx峰值及η/A峰值减小趋势变缓。以上情况表明，填充结构可以有效减缓液体晃荡的程度，减小晃荡产生的冲击压力对容器的破坏。同时，相对于改变填充结构参数，改变容器内填充区域的范围能取得更好的减晃效果，在实际工程中也易于实现。

图6描述了填充系数为ε=0.8、f=0.2和两种波数(kb=1.0，2.0)条件下，填充半径a变化对晃荡力的影响。图中可以看出，随着填充半径的增加，系统晃荡力与容器壁处晃荡力变化趋势一致，数值差别很小；填充半径a增加时，填充结构边壁处晃荡力变化较为平稳，且始终远小于系统晃荡力与容器壁处晃荡力。同时，由图6(a)可知，在k=1.0时，随着晃荡半径的增加，系统晃荡力、容器壁处晃荡力及填充结构边壁处晃荡力增加较为平稳；由图6(b)可知，在k=2.0时，随着填充半径a的增加，系统晃荡

图5 不同填充半径下，相对波数对晃荡的影响

图6晃荡力随填充半径a的变化

力与容器边壁晃荡力先减小后增大，并且在a=0.6左右时晃荡力值取得最小。

4 结 论

(1) 与解析解结果对比表明，比例边界有限元方法在很少的离散单元即可得到很高的计算精度，验证了本文采用方法在液体晃荡问题研究的精确性及高效性。

(2) 填充半径a一定时，孔隙度系数ε和线性阻力系数f的变化，对无量纲晃荡力峰值及峰值对应相对波数有显著的影响；所研究四组填充参数情况中，采用ε=0.45、f=2.0时，晃荡力峰值取得最小值。

(3) 随着填充半径a的增加，无量纲晃荡力峰值及峰值对应相对波数都随之减小；当填充半径a大于0.5并继续增大时，峰值减小趋势变缓，即此时增大填充区域后减晃的效果不明显。