任意转动激励下液体对矩形容器作用的精确解

An Exact Solution of Fluid Action on Rectangular Tank Subjected to Arbitrary Pitching Excitation

摘要: 矩形贮液容器承受转动运动激励时，液体会产生晃动，从而对容器及相关的结构产生附加的动力效应，对相关结构安全性、稳定性以及动力学行为产生重要影响。基于晃动动力学的一般方程，在小晃动假设下，本文采用线性势流理论，给出了矩形贮液容器在纯转动激励下的液体晃动线性化控制方程，通过求解液体运动方程得到液体速度势函数，由Bernoulli方程得到矩形容器内液体压力分布的精确解，通过压力积分进而得到液体对矩形容器的动态力与翻转力矩精确解，进行了液体共振激励的响应计算。精确解的适用范围为转动激励下的线性晃动问题，研究结果可为液体晃动计算及相关工程设计提供参考。

Abstract: When a rectangular liquid-storage tank is subjected to a pitching excitation, a slosh of the contained liquid is excited. The additional dynamic effects of sloshing will be exerted on the tank and the associated structure, which will lead to an important influence on the safety, stability and dynamic behavior of the associated structure. Based on the general equations of sloshing dynamics, under the assumption of small sloshing, this paper uses the linear potential flow theory to derive the linearized governing equations of liquid sloshing subjected to a pure rotating excitation. The velocity potential function is obtained through solving the governing equations of liquid motion. The exact solution of the liquid pressure distribution in the rectangular container is solved by using Bernoulli equation. The exact dynamic force and the overturning moment that act on the tank are further obtained through a pressure integral. A numerical example of liquid resonance response of is calculated. The application cope of present exact solution is limited to the linear sloshing problem subjected to a pitching excitation. The present result can be used as a reference for sloshing computation and the related engineering designs.

1. 引言

贮液箱内带有自由液面的液体在受外部运动激励下会产生晃动，晃动问题广泛存在于船舶、土木水利工程等领域。工程设计中需要知道液动压力分布以及作用在贮液箱上的液动力、力矩，这些液体的动力学效应会对贮液箱及其相关结构安全性、稳定性以及动力学行为产生重要影响。

大型运输LNG (Liquefied Natural Gas)液体的船舶在波浪作用下，由于船体的运动将引发舱内巨量LNG液体的晃动，晃动液体将会对船体产生附加的力和翻转力矩，可能影响船舶结构及航行安全，晃动引起的载荷效应已成为载液船舶安全性评估的重要内容之一。土木水利工程中的架空水箱、渡槽、升船机等结构物在地震作用下，其结构内的水体也会发生晃动，这些结构的抗震设计需要计算水体晃动对结构产生的附加液动荷载。

矩形容器是工程中使用最广泛的一类贮液箱，当矩形容器承受水平加速度运动时(例如：水平地震作用)，其液体的晃动问题得到了充分的研究 [1] [2] [3] [4]。当矩形容器承受转动运动激励时，Ibrahim [2] 在其专著中并未给出一般转动激励下的液体晃动解答；Faltinsen与Timokha [3] 在其著作中给出了一般激励(包括转动激励)下液体的晃动控制方程，提出了求解液体的晃动方程的一般步骤，然而文献 [3] 并未给出矩形容器承受转动激励时的精确解。

本文根据转动激励下的晃动控制方程，在小晃动假设下，采用线性势流理论，推导矩形容器内液体的动水压力分布，进而得到液体对矩形容器的动态力与翻转力矩，供相关工程设计参考。

2. 转动激励下液体晃动的控制方程及解答

图1显示了一个转动的刚性矩形容器，其中充有部分液体，液体静止时的宽度为2a，水深为H，设容器内壁厚度(垂直于纸面)为B。坐标 $oxyz$ 固定在容器上，坐标 $oxyz$ 平面与静止液面重合，坐标原点位于静止液面几何中心，其中符号 $\Omega$ 、 $\partial S_f$ 及 $\partial S_w$ 分别表示液体区域、自由液面及湿边界。设坐标系 $o^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ 为固定在地球上的坐标系，假定坐标系 $oxyz$ 在初始时与坐标系 $o^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ 重合，容器发生了绕 $o^{\prime }{y}^{\prime }$ 轴的转动，设转动位移、转动速度，转动加速度分别为 $\theta$， $\stackrel{\dot{}}{\theta }$， $\ddot{\theta }$，设顺时针转动为正。假定液体为无粘性、不可压缩且无旋的液体，液体在纯转动作用下的二维线性晃动方程可以写为 [3] ：

Figure 1. Rotating rigid rectangular tank filled with fluid

图1. 转动的充液刚性矩形容器

(1)

${\frac{\partial \Phi }{\partial n}|}_{\partial S_w}=\left(\omega \times r\right)\cdot n$ (2)

${\frac{\partial \Phi }{\partial n}|}_{\partial S_f}={\left[\left(\omega \times r\right)\cdot n+\frac{\partial h}{\partial t}\right]}_{\partial S_f}$ (3)

${\left[\frac{\partial \Phi }{\partial t}-\left(\omega \times r\right)\cdot \nabla \Phi \right]}_{\partial S_f}+gh=0$ (4)

式中：t为时间， $\Phi =\Phi \left(x,z,t\right)$ 为总的速度势函数， $h=h\left(x,t\right)$ 为自由液面上的竖向位移，g为重力加速度， $n$ 为湿边界的法向量， $\partial /\partial n$ 表示法向求导。设 $e_i$ 与 ${e^{\prime }}_i$ ( $i=1,2,3$ )分别为坐标系 $oxyz$ 及 $o^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ 下的单位基矢，则 $r=x{e}_1+z{e}_3$ 为任一液体质点的径向矢量， $\omega =\omega \left(t\right)=\stackrel{\dot{}}{\theta }\left(t\right)e_2$ 为容器绕y轴转动的角速度矢量，其中 $\stackrel{\dot{}}{\theta }\left(t\right)$ 为小的激励。液体压力可以由Bernoulli方程表达为：

$p=p\left(x,z,t\right)=-\rho \frac{\partial \Phi }{\partial t}-\rho U_g,\left(x,z\right)\in \Omega$ (5)

式中：ρ为液体质量密度， $U_g=-g\cdot r$ 为重力势，其中 $g=-g{e^{\prime }}_3$，将 $r=x{e}_1+z{e}_3$， $\omega =\stackrel{\dot{}}{\theta }\left(t\right)e_2$ 及 $g=-g{e^{\prime }}_3$ 代入方程(1)~(5)，注意到： ${e^{\prime }}_3=e_3\cos \theta -e_1\sin \theta$，在小转动位移下： ${e^{\prime }}_3\approx e_3-\theta e_1$，所以以上方程(1)~(5)变为：

${\nabla }^2\Phi =0,\left(x,z\right)\in \Omega$ (6a)

${\frac{\partial \Phi }{\partial x}|}_{x=\pm a}=\stackrel{\dot{}}{\theta }z$ (6b)

${\frac{\partial \Phi }{\partial z}|}_{z=-H}=-\stackrel{\dot{}}{\theta }x$ (6c)

${\frac{\partial \Phi }{\partial z}|}_{z=0}=-\stackrel{\dot{}}{\theta }x+\frac{\partial h}{\partial t}$ (6d)

${\frac{\partial \Phi }{\partial t}|}_{z=0}+gh-gx\theta =0$ (6e)

$p=-\rho \frac{\partial \Phi }{\partial t}-\rho gz+\rho xg\theta ,\left(x,z\right)\in \Omega$ (6f)

由傅里叶级数及级数收敛的必要性条件，设速度势函数的解为：

$\Phi =\stackrel{\dot{}}{\theta }\phi \left(x,z\right)+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{{\stackrel{\dot{}}{q}}_n\left(t\right)\cdot \cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \frac{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}\left(z+H\right)\right]}{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}\right\}}$ (7)

将式(7)代入方程(6a)~(6c)，求解方程，得速度势函数为：

$\begin{array}{c}\Phi =\stackrel{\dot{}}{\theta }xz+\stackrel{\dot{}}{\theta }{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{\frac{32a^2}{{\left[\left(2n-1\right)\pi \right]}^3\cosh \left[\frac{\left(2n-1\right)\pi H}{4a}\right]}\cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \sinh \left[\frac{\left(2n-1\right)\pi }{2a}\left(z+\frac{H}{2}\right)\right]\right\}}\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{{\stackrel{\dot{}}{q}}_n\left(t\right)\cdot \cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \frac{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}\left(z+H\right)\right]}{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}\right\}}\end{array}$ (8)

将式(8)代入方程(6d)和(6e)，然后合并两个方程，得：

$\begin{array}{l}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \left\{{\stackrel{⃛}{q}}_n\left(t\right)+g\left[\frac{\left(2n-1\right)\pi }{2a}\right]\tanh \left[\frac{\left(2n-1\right)\pi }{2a}\right]{\stackrel{\dot{}}{q}}_n\left(t\right)\right\}}\\ =-gx\stackrel{\dot{}}{\theta }-g\stackrel{\dot{}}{\theta }{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{\frac{16a}{{\left[\left(2n-1\right)\pi \right]}^2}\cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\right\}}\\ -\stackrel{⃛}{\theta }{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{\frac{32a^2}{{\left[\left(2n-1\right)\pi \right]}^3}\tanh \left[\frac{\left(2n-1\right)\pi H}{4a}\right]\cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\right\}}\end{array}$ (9)

在零初始条件下( $q_n\left(0\right)={\ddot{q}}_n\left(0\right)=\theta \left(0\right)=\ddot{\theta }\left(0\right)=0$ )，对方程(9)两边对时间t积一次分，再将方程(9)两边乘以 $\cos \left[\frac{\left(2m-1\right)\pi }{2a}\left(x+a\right)\right]$ ( $m=1,2,3,\cdots$ )，方程两边对x积分得：

${\ddot{q}}_n\left(t\right)+{\omega }_n^2{q}_n\left(t\right)=A_{1n}\theta \left(t\right)+A_{2n}\ddot{\theta }\left(t\right),n=1,2,3,\cdots$ (10)

式中：

${\omega }_n^2=g\left[\frac{\left(2n-1\right)\pi }{2a}\right]\tanh \left[\frac{\left(2n-1\right)\pi H}{2a}\right]$ (11)

$\{\begin{array}{l}{A}_{1n}=-\frac{8ga}{{\left[\left(2n-1\right)\pi \right]}^2}\\ A_{2n}=-\frac{32a^2}{{\left[\left(2n-1\right)\pi \right]}^3}\cdot \tanh \left[\frac{\left(2n-1\right)\pi H}{4a}\right]\end{array}$ (12)

根据Duhamel积分，方程(10)的解为：

$q_n\left(t\right)=\frac{A_{1n}}{{\omega }_n}{\displaystyle \underset{0}{\overset{t}{\int }}\theta \left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }+\frac{A_{2n}}{{\omega }_n}{\displaystyle \underset{0}{\overset{t}{\int }}\ddot{\theta }\left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }$ (13)

${\ddot{q}}_n\left(t\right)=-A_{1n}{\omega }_n{\displaystyle \underset{0}{\overset{t}{\int }}\theta \left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }-A_{2n}{\omega }_n{\displaystyle \underset{0}{\overset{t}{\int }}\ddot{\theta }\left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }+A_{1n}\theta \left(t\right)+A_{2n}\ddot{\theta }\left(t\right)$ (14)

所以液体压力为：

$\begin{array}{l}p\left(x,z,t\right)=-\rho \frac{\partial \Phi }{\partial t}-\rho gz+\rho xg\theta \\ =-\rho \ddot{\theta }xz-\rho \ddot{\theta }{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{\frac{32a^2}{{\left[\left(2n-1\right)\pi \right]}^3\cosh \left[\frac{\left(2n-1\right)\pi }{\text{4}a}H\right]}\cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \sinh \left[\frac{\left(2n-1\right)\pi }{2a}\left(z+\frac{H}{2}\right)\right]\right\}}\\ -\rho {\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{{\ddot{q}}_n\left(t\right)\cdot \cos \left[\frac{\left(2n-1\right)\pi }{2a}\left(x+a\right)\right]\cdot \frac{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}\left(z+H\right)\right]}{\cosh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}\right\}}-\rho gz+\rho xg\theta \end{array}$ (15)

于是液体对容器x方向的合力 $F_x\left(t\right)$ 及对容器底部中点 $A^{\prime }$ 的翻转力矩 $M_{A^{\prime }}$ 可以由下式积分得到：

$\begin{array}{c}{F}_x\left(t\right)={\displaystyle \underset{-H}{\overset{0}{\int }}\left[{p\left(x,z,t\right)|}_{x=a}-{p\left(x,z,t\right)|}_{x=-a}\right]B\text{d}z}=\langle \rho a{H}^2+4\rho a{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{A_{2n}\tanh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}{\left(2n-1\right)\pi }}\rangle B\ddot{\theta }\left(t\right)\\ +\left\{2\rho gaH+4\rho a{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{A_{1n}\tanh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}{\left(2n-1\right)\pi }}\right\}B\theta \left(t\right)\\ -4\rho aB{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{A_{2n}{\omega }_n\tanh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}{\left(2n-1\right)\pi }}{\displaystyle \underset{0}{\overset{t}{\int }}\ddot{\theta }\left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }\\ -4\rho aB{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{A_{1n}{\omega }_n\tanh \left[\frac{\left(2n-1\right)\pi }{2a}H\right]}{\left(2n-1\right)\pi }}{\displaystyle \underset{0}{\overset{t}{\int }}\theta \left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }\end{array}$ (16)

$\begin{array}{c}{M}_{A^{\prime }}\left(t\right)={\displaystyle \underset{-H}{\overset{0}{\int }}\left[{p\left(x,z,t\right)|}_{x=a}-{p\left(x,z,t\right)|}_{x=-a}\right]\left(z+H\right)B\text{d}z}+{\displaystyle \underset{-a}{\overset{a}{\int }}{p\left(x,z,t\right)|}_{z=-H}xB\text{d}x}\\ =B\langle \frac{1}{3}\rho a{H}^3+\frac{2}{3}\rho a^{3}H+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{128\rho a^{3}H}{{\left[\left(2n-1\right)\pi \right]}^4}\left\{1-\frac{6}{\left[\left(2n-1\right)\pi \right]}\left(\frac{a}{H}\right)\tanh \left[\frac{\left(2n-1\right)\pi H}{4a}\right]\right\}\begin{array}{c}\\ \\ \\ \end{array}}\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{8\rho a^2{A}_{2n}}{{\left[\left(2n-1\right)\pi \right]}^2}\frac{1}{\text{cosh}\left[\frac{\left(2n-1\right)\pi H}{2a}\right]}}+\text{2}\rho {\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{A}_{2n}{B}_n}\rangle \ddot{\theta }(t)\end{array}$

$\begin{array}{l}+B\left(\rho ga{H}^2+\frac{2}{3}\rho g{a}^3+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{8\rho a^2{A}_{1n}}{{\left[\left(2n-1\right)\pi \right]}^2}\frac{1}{\text{cosh}\left[\frac{\left(2n-1\right)\pi H}{2a}\right]}}+2\rho {\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{A}_{1n}{B}_n}\right)\theta \left(t\right)\\ -B{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{2\rho A_{1n}{B}_n{\omega }_n+\frac{8\rho a^2{A}_{1n}{\omega }_n}{{\left[\left(2n-1\right)\pi \right]}^2}\frac{1}{\text{cosh}\left[\frac{\left(2n-1\right)\pi H}{2a}\right]}\right\}{\displaystyle \underset{0}{\overset{t}{\int }}\theta \left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }}\\ -B{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left\{2\rho A_{2n}{B}_n{\omega }_n+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{8\rho a^2{A}_{2n}{\omega }_n}{{\left[\left(2n-1\right)\pi \right]}^2}\frac{1}{\text{cosh}\left[\frac{\left(2n-1\right)\pi H}{2a}\right]}}\right\}{\displaystyle \underset{0}{\overset{t}{\int }}\ddot{\theta }\left(\tau \right)\sin \left[{\omega }_n\left(t-\tau \right)\right]\text{d}\tau }}\end{array}$ (17)

其中：

$B_n=\frac{2aH}{\left(2n-1\right)\pi }\tanh \left[\frac{\left(2n-1\right)\pi H}{2a}\right]-\frac{4a^2}{{\left[\left(2n-1\right)\pi \right]}^2}\left\{1-\frac{1}{\cosh \left[\frac{\left(2n-1\right)\pi H}{2a}\right]}\right\}$ (18)

3. 数值算例

由1) 方程的正确，2) 得到的解答满足方程，3) 级数解答收敛可知其精确解的可靠性。现设有一矩形水箱，静止水截面尺寸2a = 8 (m)，H = 6 (m)，容器内壁(垂直于)厚度为B = 1 (m)，流体质量密度为 $\rho =1000\left(\text{kg}/{\text{m}}^{\text{3}}\right)$，液体一阶晃动自然圆频率为1.944 (rad/s)，设容器的转动(共振)位移激励为谐波函数 $\theta \left(t\right)=0.1\sin \left(1.944t\right)\left(\text{rad}\right)$，则角加速度激励为： $\ddot{\theta }\left(t\right)=-0.3764\sin \left(1.944t\right)\left(\text{rad}/{\text{s}}^{\text{2}}\right)$。取前三阶模态叠加计算得到作用在容器上液体力 $F_x\left(t\right)$ 及翻转力矩 $M_{A^{\prime }}\left(t\right)$ 如图2与图3所示，表明容器内的液体发生了共振相应。

Figure 2. The force $F_x\left(t\right)$ acting on the tank by fluid

图2. 液体作用在容器上的力 $F_x(t)$

Figure 3. The overturning moment $M_{A^{\prime }}\left(t\right)$ acting on the tank by fluid

图3. 液体作用在容器上的翻转力矩 $M_{A^{\prime }}(t)$

4. 结语

在小晃动假设下，导得了纯转动激励下的液体晃动线性化的控制方程，得到了矩形容器内液体压力分布的精确解，通过积分得到液体对矩形容器的动态力与翻转力矩的精确解，本文的精确解的适用范围为转动激励下的线性晃动问题，本文解答可供相关工程设计计算参考。

基金项目

国家自然科学基金面上项目(51879191)。

NOTES

^*通讯作者。