1.       Bernoulli积分

前面讨论的充液系统动力学问题均限于液体在腔内完全充满情形。如不充满就有自由液面出现，在重力或惯性力作用下产生晃动 (sloshing)。对液体晃动的研究先假定刚体的运动已预先确定，使问题限制在流体力学范畴以内。然后再考虑液体晃动对刚体运动的影响。实际上刚体与液体的运动是相互影响同时发生的，与全充液情况相比，非全充的刚液耦合系统的动力学分析难度明显增加。

设 V  和 V′ 为液体和刚体占据的空间，Σ 和 Σ′ 为腔壁的湿润和非湿润表面，S 为自由液面（图 1）。刚体内 O 点的速度为 v₀，瞬时角速度为 ω，腔内流场在 P 点处的流速为 u，则 P 点处质点的绝对速度 v 为

                                                                                                                        (1)

图 1  带充液腔的刚体

博文 “充液系统动力学(一)”  中，曾基于分析力学的虚功率原理导出流体动力学方程。若忽略流体的粘性，称为欧拉流体动力学方程：

                                                                    (2)

将式 (1) 中解出的 u 代入上式，设无旋流体的速度势函数为 ϕ，作用力 f 为有势力，以 U 为力函数。令v =▽ϕ，f = ▽U，得到

              (3)

直接验算可证明以下关系式成立：

                                            (4)

利用上式将方程 (3) 化作

                                                         (5)

从而导出 Bernoulli 积分：

                                                                        (6)

上式右边的任意时间函数 C(t) 可并入左边的势函数 ϕ 而不影响流场的速度分布，可用积分常数 C 代替。如流速很小，忽略 v 的二次项，简化为

                                                                                                                                     (7)

式 (6) 或 (7) 也称为 Bernoulli 定理或 Bernoulli 方程。

2.       液面的表面张力

液体内部的分子引力在界面上表现为表面张力。在腔体体积足够大，或流体的体积力足够大的一般情况下可予以忽略。讨论失重或微重状态下的充液刚体运动时，液面的表面张力是不可忽略的因素。

自由液面 S 为液体与气体之间的界面，设液面上单位长度范围内的表面张力为 α，称为液体与气体之间的表面张力系数。在液面 S 上 的任意点 P 处，沿液面主曲率方向取边长为 ds₁ 和 ds₂ 的微元弧段组成微元液面 dσ =  ds₁ds₂。设液面沿主曲率方向的曲率半径为 R_j (j = 1,2)，则弧段 ds₁ 和 ds₂ 相对曲率中心 O_j (j = 1,2) 的张角为 φ_j = ds_j/R_j  (j = 1,2) （图 2）。

图 2   弧段 ds₁ 和 ds₂ 的张角

在 P 点处沿曲面法向取速度变分 δu_n = n·δu，使 P点处的微元液面 dσ 产生表面张力，与边线的长度变化成正比，导致虚功率增量 δP（图 3）：

                                                                   (8)

其中 K = (1/R₁)  + (1/R₂) 为 P 点处液面的曲率，为主曲率平均值的二倍。将上式在自由液面 S 范围内积分，积分值仍以 δP 表示：

                                                                                                                             (9)

图 3  法向虚速度产生的表面张力虚功率

        在自由液面与刚体腔壁相交的边缘 s_c 处，即刚体、液体、气体三种介质的界面 S,  Σ,  Σ′ 的公共交线处，除液体与气体界面 S 的虚功率增量 (9) 以外，还必须考虑因变分 δu_n 引起的刚体与液体的界面 Σ 和 刚体与气体的界面 Σ′ 的周边长度变化所产生的虚功率增量。分别以 α₁, α₂ 表示刚体与液体、刚体与气体之间的表面张力系数。设在界线 s_c 处，界面 S 与 Σ 的法线倾角为 θ_c，称为接触角。从图 4 可看出，虚速度  δu_n 在 s_c 处引起液体与气体、刚体与液体、刚体与气体之间的界线变化分别与 δu_n cotθ_c, δu_n cscθ_c, -δu_n cscθ_c 成正比，所导致的虚功率分别为 αδu_n cotθ_c, α₁δu_n cscθ_c, -α₂δu_n cscθ_c 。总的虚功率增量为

                                                                                     (10)

图 4   法向虚速度引起界面边界的长度变化

在博文 “充液系统动力学(一)” 中应用虚功率原理列写刚体与液体耦合系统的动力学方程时，若考虑液面的表面张力因素，应将式 (9), (10) 表示的虚功率增量加入此方程，化作

                       (11)

基于变分 δu_n 的任意性，导出考虑表面张力的自由液面边界条件：

                                                                                             (12)

3.       液体的静平衡自由液面

考虑液面的表面张力，将 Bernoulli 积分 (7) 中的压强 p 以式 (12) 代入， p₀ 并入积分常数 C，力函数 U  改以势能的负值 -V 表示，化作

                                                                                               (13)

设液面 S 处于静止状态，令上式中的速度项 ∂ϕ/∂t = 0 等于零。利用柱坐标 (r, θ, z) 表示流场中任意点  P  的位置（图 5），微重力场和离心力场的势能分别为

                             (14)

图 5   液体流场的柱坐标

其中 cg 为微重力（地球上 c = 1, 太空中 c < 1），ω₀ 为容器绕 z 轴的自旋角速度。将式 (14a) 代入式 (13)，选择某个特征长度 a，化作无量纲形式：

                                                                                                             (15)

各无量纲量定义为

                                                                                           (16)

其中 B 称为  Bond 数，其物理意义为体积力与表面张力之比。表面张力仅当  Bond  数足够大时方允许忽略，此时液面方程 (13) 简化为 V = const，自由液面为平面 z = z₀（重力场），或圆柱面 r = r₀（离心力场）。如重力与离心力共同作用，则为旋转抛物面  z = z₀+(ω₀²/2cg)r²（图 6）。

图 6   静平衡自由液面

4.       液体晃动的动力学方程

设刚体处于静止状态，势能 V 仅含微重力场 cgz 。设腔壁相对 z 轴对称，液体的平衡液面亦相对 z 轴对称。平衡液面和腔壁的曲面方程分别为

                                                                                                                 (17)

设晃动引起液面高度的摄动规律为 h(r, θ, t)，则受扰液面的曲面方程为（图 7）

                                                             (18)

图 7   腔体与液面的曲面形状

利用曲率 K 的计算公式，将式 (17) 代入计算，以下角标 f_r,  h_r,  h_θ 表示  f  和 h 对 r 和 θ 的偏导数。仅保留 h 及其导数的一次项，将导出的曲率 K 和式 (18) 代入 Bernoulli 方程 (13)，得到自由液面 S 处的动力学条件：

                                    (19)

令 h = 0，即得到无晃动时考虑表面张力的自由液面平衡条件：

                                                                          (20)

从方程 (19) 减去上式，化作扰动 h(r, θ, t) 的动力学方程：    

             (21)      

以 F(r, θ, z, t) = 0 表示腔壁 Σ 或液面 S 的曲面方程。此曲面法线的方向余弦为

                                                                                                   (22)

将上式代入腔壁  Σ  或液面 S 处无法向位移条件 u_n = u·n = 0 ，得到约束方程：

                                                                                           (23)

其中 u_r,  u_θ,  u_z 为相对流速 u 沿柱坐标的投影：

                                                                                                 (24)

对于腔壁 Σ 令 F = z - w(r)，对于液面 S 令 F = z - f(r) - h(r, θ, t)，代入式 (23) 且利用式 (24)，导出 Σ 和 S 处的运动学条件。

                                                                                                         (25)

综上所述，液面晃动问题归纳为以下 Neumann 边值问题：

 (26)

其中还应包括界面交线 s_c 处接触角 θ 的边界条件，为避免过多数学推导此处被省略。利用数值方法从上述边值问题可解出 ϕ(r, θ, t) 和 h(r, θ, t) 。若考虑刚体的牵连运动，则势函数V中应增加牵连加速度产生的惯性力势能，且应考虑牵连运动对腔壁Σ和液面S处运动学条件的影响。

5.       液体晃动的等效力学模型

以上对液体晃动的分析仅限于刚体静止或做独立运动条件下的流体力学计算。更符合实际的分析必须考虑刚体运动与液体晃动之间的耦合，理论分析的难度明显增大。1978 年俄罗斯学者 Mikishev 提出用一组单摆和阻尼器的简单模型，作为考虑刚体与液体耦合效应时液体晃动的等效力学模型^[1]。模型中的固定质量块模拟被刚体带动的底部液体质量，单摆模拟液面的晃动。单摆的固有频率可依据上述边值问题的计算结果。也可利用实验方法，在相同形状缩小的充液腔体的实验中测出液体的晃动频率及对腔壁的作用力和力矩，以确定图8表示的等效模型中固定质量块的惯性矩和位置、单摆的质量和摆长、悬挂点位置等参数。这种等效模型的处理方法适合于在复杂控制系统的工程设计中使用。

图 8   液体晃动的等效力学模型

关于充液刚体动力学的系列博文至此结束。有兴趣的读者可参阅系统性论述该领域的专著如[2],[3]，以及国内外近年来发表的丰富文献。

参考文献：

[1] Mikishev,G.N. Experimental methods of spacecraft dynamics. Moscow: Mechanical Manufacture, 1978. p.29

[2] Moiseyev,N.N., Rumyantsev,V.V. Dynamic stability of bodies containing fluid. Berlin: Springer-Verlag, 1968

 [3] Ibrahim,R.A. Liquid sloshing dynamics: theory and applications. Cambridge: Cambridge University Press, 2005

     （改写自：王照林，刘延柱. 充液系统动力学. 第4章. 北京：科学出版社，2002)