重心在上却不会翻 — 稳心高度 GM 与自由液面效应

GM 是水线面惯性矩挣来的余量。舱内一个自由液面就可能把这份余量整个拿走。

集装箱船的重心比水下部分的中心高出好几米。重的一端压在上面，船却能自己立起来。反过来，重心低于浮心的驳船也会翻。可见仅凭上下关系什么都定不了。本文要讲的是真正做出判定的那个量 — 稳心高度 $GM$ — 它从何而来，以及舱内液体为什么能把这个值整个抹掉。

重心位于浮心之上

先统一符号。 $K$ 是龙骨(keel，船底中心)， $G$ 是重心， $B$ 是浮心。浮心是水下体积的形心(centroid)。阿基米德给出的只是大小：浮力等于排开水的重量，作用线通过排开体积的形心。

正浮时 $G$ 与 $B$ 在同一条铅垂线上。重力与浮力大小相等方向相反，力矩为零。此时 $G$ 在上还是在下都不起作用。稳定性不是由静止状态定义的，而是由稍微倾斜之后出现了什么定义的。

一旦倾斜情况就变了。 $B$ 是水下形状的形心，形状变了它就跟着移动。 $G$ 只随船体一起转动，在船内位置不变。两条铅垂线错开后就出现了力臂。这个力臂称为复原力臂(righting arm) $GZ$ ，对排水量 $\Delta$ 而言复原力矩为 $M = \Delta \cdot GZ$ 。

一倾斜，水下形状就变了

在下面的模拟中亲手倾斜一下。改变船宽、吃水和重心高度，再按 release at 15 deg，看它能不能回正。

beam B [m]16.0draft T [m]4.0KG [m]6.6

phi = 15.0 deg

Raise KG past 7.33 m and the barge no longer comes back — it rolls out to an angle of loll and stays there. Hold the heel and drag it past 27 deg: the deck edge goes under, the waterplane stops widening, and the green GZ curve peels away from the dashed straight line GM sin phi. That gap is why a single number cannot certify a hull.

$B$ (绿色)被拉向低舷一侧，它与通过 $G$ (黄色)的铅垂线之间的间距就是 $GZ$ 。紫色的 $M$ 是浮力作用线与船体中心线的交点。把 KG 滑块推到 7.33 m 以上， $M$ 就落到 $G$ 之下，船不再回正，而是歪向一侧停住。

$BM = I/\nabla$ — 全在水线面的惯性矩#

在小倾角 $\phi$ 下，水下形状的变化可以归结为两个楔形：一舷露出水面的楔形，另一舷新浸入的楔形。排水量不变，因此两者体积相等，而左右分开。正是这次体积转移把 $B$ 挪走了。

设水线面(waterplane，水面切过船体的截面)上距中心线 $x$ 处有微元条带 $dA$ ，它带来的体积变化是 $x\phi\,dA$ ，力矩是 $x^2\phi\,dA$ 。全部相加就得到截面惯性矩。

\overline{BB'} = \frac{\phi \int_A x^{2}\,dA}{\nabla} = \frac{I\,\phi}{\nabla}, \qquad BM = \frac{I}{\nabla}

$I$ 是水线面对中心线的截面惯性矩， $\nabla$ 是排水体积。对宽度为 $b$ 的矩形剖面，单位长度上 $I = b^3/12$ ， $\nabla = b\,T$ ，于是

BM = \frac{b^{2}}{12\,T}

关键在于：宽度以三次方进入，又被分母的 $\nabla$ 抵消一次，最后留下平方。船宽加大 20%， $BM$ 就增大 44%。独木舟不稳而木筏稳，原因就在这里。判定式本身只有一行。

GM = KB + BM - KG

$KB$ 是浮心高度， $KG$ 是重心高度。 $GM > 0$ 则回正， $GM < 0$ 则倾覆。宽 $b = 16$ m、吃水 $T = 4$ m 的驳船， $KB = 2.0$ m， $BM = 5.333$ m，故 $KM = 7.333$ m。若 $KG$ 为 6.6 m，余量只有 0.733 m。

$GM \sin\phi$ 失效的角度#

$GM$ 是 $\phi \to 0$ 处的斜率。小角度下 $GZ \approx GM \sin\phi$ 吻合得很好。问题是这个近似何时失效。

失效的位置由几何决定，因为只有水线面不变时 $I$ 才是常数。上面那条驳船干舷(freeboard，水面以上的船体高度)为 4 m，因此 $\tan\phi = 2 \times 4/16$ ，即 26.6° 时甲板边缘入水。从那一刻起水线面不再变宽，反而收窄。 $BM$ 随之崩塌， $GZ$ 曲线开始下折。

在上面的模拟中按 hold heel，把角度推过 30°，就能看到绿色实线(精确值)从紫色虚线( $GM\sin\phi$ )上剥离。两条曲线符号分道扬镳的角度就是消失角。

舱内液体把 GM 整个拿走#

现在往船内放一个液舱。压载舱、燃油舱、货舱，乃至被消防水淹没的汽车甲板，都是同一个问题。

液体的重量不会因倾斜而改变，它只是移动。可这一移动会把 $G$ 拽向低舷一侧。计算结构与 $BM$ 完全相同：舱内自由液面同样是一个水线面，一倾斜就有两个楔形互换位置。

\delta GM = \frac{\rho_f\, i}{\rho\, \nabla}

$i$ 是舱内自由液面的截面惯性矩， $\rho_f$ 是舱内液体密度， $\rho$ 是舷外水的密度。这个量并不是真的从船上减去，而是当作 $G$ 上升了同样的高度来处理(虚拟上升，virtual rise)。**舱里装了多少液体并不出现在公式中。**只要不是灌满到没有自由液面、也不是完全排空，装 20 cm 和装 2 m 的损失是一样的。

拖动下面的舱壁滑块试试。

tank width b [m]12.0bulkheads n1fill [m]2.0KG solid [m]7.2

phi = 10.0 deg

Start at n = 1 and press release: the liquid slides to the low side and the barge never comes back. Now drag the bulkhead slider. The liquid volume never changes, the centre of gravity of the ship at rest never changes — only the width of each free surface does, and at n = 2 the red loss bar has already dropped to a quarter.

黄色是舱内液体。无论船体怎么倾，液面始终保持水平。在 n = 1 时松手，船不会回正，而是歪在一侧保持平衡(固定横倾，angle of loll)。加一道舱壁，红色损失条立刻降到四分之一，船就立起来了。

一道舱壁把损失降到四分之一

把宽度为 $b_t$ 的液舱分成 $n$ 格，每格宽 $b_t/n$ ，共 $n$ 格。

i = n \cdot \frac{(b_t/n)^{3}}{12} = \frac{b_t^{3}}{12\,n^{2}}

损失按 $1/n^2$ 而非 $1/n$ 下降：宽度的三次方压过了格数的一次方。上面那条驳船若有一个宽 12 m 的液舱，则 $i = 144$ m⁴， $\nabla = 64$ m³/m，于是 $\delta GM = 2.195$ m。0.733 m 的 $GM$ 早就没了。只要在中线加一道舱壁，损失降到 0.549 m，船就重新站住。

油轮中线上那道纵舱壁就是为此而设。载货量不变，重心高度也不变，变的只是每个自由液面的宽度减半。

用 Python 画出的 GZ 曲线与消失角#

$GM$ 不过是一条切线的斜率。大角度下的 $GZ$ ，直接切出水下多边形再求形心反而更快。把船体剖面写成多边形，用二分法找出使水下面积等于排水面积的水线高度，再计算剩下那块的形心。

import math
 
def poly_area_centroid(poly):
    a = cx = cy = 0.0
    n = len(poly)
    for i in range(n):
        x0, y0 = poly[i]
        x1, y1 = poly[(i + 1) % n]
        cr = x0 * y1 - x1 * y0
        a += cr
        cx += (x0 + x1) * cr
        cy += (y0 + y1) * cr
    a *= 0.5
    if abs(a) < 1e-14:
        return 0.0, 0.0, 0.0
    return a, cx / (6 * a), cy / (6 * a)
 
def clip_below(poly, phi, h):
    """在倾角 phi 下，只保留大地坐标高度不超过 h 的部分"""
    def f(p):
        return p[1] * math.cos(phi) - p[0] * math.sin(phi) - h
    out = []
    n = len(poly)
    for i in range(n):
        p, q = poly[i], poly[(i + 1) % n]
        fp, fq = f(p), f(q)
        if fp <= 0.0:
            out.append(p)
        if (fp < 0.0) != (fq < 0.0):
            t = fp / (fp - fq)
            out.append((p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1])))
    return out
 
def waterline_offset(poly, phi, area0):
    """用二分法求使水下面积等于 area0 的水线高度"""
    lo, hi = -60.0, 60.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        a, _, _ = poly_area_centroid(clip_below(poly, phi, mid))
        if a < area0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)
 
def righting_arm(beam, depth, draft, kg, phi):
    hull = [(-beam / 2, 0.0), (beam / 2, 0.0), (beam / 2, depth), (-beam / 2, depth)]
    h = waterline_offset(hull, phi, beam * draft)
    _, xb, yb = poly_area_centroid(clip_below(hull, phi, h))
    return xb * math.cos(phi) + (yb - kg) * math.sin(phi), xb, yb
 
def gz_curve(beam, depth, draft, kg, deg_max=70, step=1):
    return [(d, righting_arm(beam, depth, draft, kg, math.radians(d))[0])
            for d in range(0, deg_max + 1, step)]
 
BEAM, DEPTH, DRAFT = 16.0, 8.0, 4.0
KG = 6.6
KB = DRAFT / 2
BM = BEAM ** 2 / (12 * DRAFT)
GM = KB + BM - KG
print(f"KB={KB:.3f} m  BM=B^2/12T={BM:.3f} m  KM={KB+BM:.3f} m  KG={KG:.3f} m  GM={GM:.3f} m")
print(f"deck edge immerses at phi = {math.degrees(math.atan(2*(DEPTH-DRAFT)/BEAM)):.1f} deg")
print()
print(" phi[deg]   GZ[m]   GM*sin(phi)[m]   error[%]")
for d, gz in gz_curve(BEAM, DEPTH, DRAFT, KG, 60, 5):
    lin = GM * math.sin(math.radians(d))
    err = f"{100 * (lin - gz) / gz:8.1f}" if abs(gz) > 0.01 else "       -"
    print(f"  {d:5.1f}   {gz:7.4f}   {lin:11.4f}   {err}")
 
curve = gz_curve(BEAM, DEPTH, DRAFT, KG, 70, 1)
dmax, gzmax = max(curve, key=lambda t: t[1])
vanish = next((d for d, g in curve if d > 5 and g <= 0.0), None)
print()
print(f"max GZ = {gzmax:.3f} m at {dmax} deg,  vanishing stability at {vanish} deg")

KB=2.000 m  BM=B^2/12T=5.333 m  KM=7.333 m  KG=6.600 m  GM=0.733 m
deck edge immerses at phi = 26.6 deg
 
 phi[deg]   GZ[m]   GM*sin(phi)[m]   error[%]
    0.0    0.0000        0.0000          -
    5.0    0.0657        0.0639       -2.7
   10.0    0.1417        0.1273      -10.2
   15.0    0.2394        0.1898      -20.7
   20.0    0.3716        0.2508      -32.5
   25.0    0.5550        0.3099      -44.2
   30.0    0.7207        0.3667      -49.1
   35.0    0.6823        0.4206      -38.4
   40.0    0.5196        0.4714       -9.3
   45.0    0.2828        0.5185       83.3
   50.0    0.0001        0.5618          -
   55.0   -0.3116        0.6007     -292.8
   60.0   -0.6406        0.6351     -199.1
 
max GZ = 0.727 m at 31 deg,  vanishing stability at 51 deg

5° 处误差为 2.7%。到 30°， $GM\sin\phi$ 只有真实 $GZ$ 的一半，因为在甲板入水之前水线面一直变宽，真实复原力超过了线性预测。可越过 45° 之后符号反转，这时线性近似就往危险的一侧说谎：50° 时真实 $GZ$ 已经为零， $GM\sin\phi$ 却仍然承诺 0.56 m。

自由液面损失不用切多边形也能核对。甲板入水之前舷侧是垂直的，复原力臂有闭式解。损失是虚拟上升，把 $\delta GM$ 从 $GM$ 里减掉再代入即可。

import math
 
BEAM, DRAFT, KG = 16.0, 4.0, 6.6
NABLA = BEAM * DRAFT
BM = BEAM ** 2 / (12 * DRAFT)
GM = DRAFT / 2 + BM - KG
 
def free_surface_loss(rho_f, rho, b, n, nabla):
    """宽 b 的液舱被 n 道舱壁分隔后的自由液面损失 [m]"""
    i_free = n * (b / n) ** 3 / 12
    return rho_f * i_free / (rho * nabla)
 
def wall_sided_gz(gm, bm, phi):
    """甲板入水之前、舷侧仍为垂直时的复原力臂"""
    return (gm + 0.5 * bm * math.tan(phi) ** 2) * math.sin(phi)
 
print(" tanks   i[m^4]   dGM[m]   GM_eff[m]   GZ(20deg)[m]")
for n in (1, 2, 3, 4):
    d_gm = free_surface_loss(1000.0, 1025.0, 12.0, n, NABLA)
    gz20 = wall_sided_gz(GM - d_gm, BM, math.radians(20))
    print(f"  {n:4d}  {n*(12.0/n)**3/12:7.1f}  {d_gm:7.3f}   {GM-d_gm:8.3f}   {gz20:11.3f}")

 tanks   i[m^4]   dGM[m]   GM_eff[m]   GZ(20deg)[m]
     1    144.0    2.195     -1.462        -0.379
     2     36.0    0.549      0.185         0.184
     3     16.0    0.244      0.489         0.288
     4      9.0    0.137      0.596         0.325

一道舱壁把 $i$ 从 144 降到 36，正好是四分之一。20° 处的复原力臂由 −0.379 m 变为 +0.184 m，符号翻转。这里用的闭式解 $GZ = (GM + \tfrac{1}{2}BM\tan^2\phi)\sin\phi$ ，在甲板入水之前与多边形计算一致：没有损失时 20° 给出 0.3716 m，与上面的表格相同。

在 6-DOF VOF 计算里把船放进水之前#

一旦用 CFD 求解浮体，上面这些数字都会变成校核基准。

第一，在设定初始吃水之前先手算 $GM$ 。用 VOF 捕捉自由液面(界面捕捉方法对比)，再挂上 6-DOF 刚体运动(基于四元数的刚体 FEM)，船自己就会摇。若这段自由衰减振荡的周期 $T_\phi = 2\pi k/\sqrt{g\,GM}$ 与手算值对不上，就该怀疑网格或者输入的惯性矩。 $k$ 是横摇回转半径。

第二，网格分辨率会直接吃掉 $BM$ 。 $BM$ 是水线面上按 $x^2$ 加权的积分，因此靠近舷侧的网格贡献最大。自由液面若在两三个网格上抹开，有效水线面宽度就模糊了，复原力矩被低估。在水线附近做局部加密是更便宜的办法。

第三，如果带着满舱计算，自由液面损失会在计算中自动出现 — 同时晃荡(sloshing)也一起出现。液体与船体运动相位错开时，结果可能比准静态的 $\delta GM$ 更糟，也可能更好。隐式处理自由水面的方法(Casulli–Zanolli 的嵌套 Newton 法)正是在这里派上用场。手算在那时依然有效，只是不再作为答案，而是作为怀疑计算结果的基线。