重心在上却不会翻 — 稳心高度 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。按单位长度把 $b^3/12$ 除以 $\nabla = b T = 64$ ，又得到这 5.333 m。船宽加大 20%， $b^2$ 变成 1.44 倍，正文里的 44% 就是这个平方。于是 $0.44 \times 5.333$ m 直接加进 $GM$ 。 $KG = 6.6$ m 不动，0.733 m 的余量就按这笔加法拉开。薄舷侧如果用壳来算，这份余量会先被锁住。若剪切闭锁把挠度压死，读 $GZ$ 之前先怀疑刚度矩阵。

$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 的损失是一样的。

损失按 $1/n^2$ 下降，原因在下一式的计数里。

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

把宽度为 $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° 处 $GZ$ 为 0.0657， $GM\sin\phi$ 为 0.0639，误差 −2.7%。25° 处是 0.5550 对 0.3099，误差拉到 −44.2%。甲板入水角 26.6° 刚过的 30°，真实 $GZ$ 是 0.7207，线性式是 0.3667，大约一半。甲板入水之前水线面一直变宽，真实复原力超过线性预测。最大复原力臂 0.727 m 在 31°。45° 处误差符号翻成 +83.3%，50° 的真实 $GZ$ 是 0.0001，线性式仍留下 0.5618。消失角 51° 之后，55° 与 60° 的真实值是 −0.3116、−0.6406。线性式那时仍是 +0.6007、+0.6351。过了消失角还承诺正的复原力臂，就是往危险一侧偏。

自由液面损失不用切多边形也能核对。甲板入水之前舷侧是垂直的，复原力臂有闭式解。损失是虚拟上升，把 $\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

把舱壁数 n 从 1 拨到 2，红色损失条降到四分之一，躺着的驳船立起来。

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)。加一道舱壁，红色损失条立刻降到四分之一，船就立起来了。

一道舱壁把 $i$ 从 144 降到 36，正好是四分之一。 $\rho_f = 1000$ 、 $\rho = 1025$ 、 $i = 144$ 、 $\nabla = 64$ 时， $\delta GM = 1000 \times 144 / (1025 \times 64) = 2.195$ m。 $n = 4$ 时 $i = 9$ ， $\delta GM = 0.137$ m，20° 复原力臂是 0.325 m。还没到没有损失时的 0.3716 m。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 计算里把船放进水之前#

把初始吃水送进求解器之前，先把 $GM = 0.733$ m 留在纸上。自由横摇周期 $T_\phi = 2\pi k/\sqrt{g\,GM}$ 和这个手算值对不上，就先怀疑网格或转动惯量。 $k$ 是横摇回转半径。 $BM$ 是水线面上按 $x^2$ 加权的积分，舷侧网格贡献最大。自由液面在两三个网格上抹开，有效宽度就缩，复原力矩变小。

满舱计算里，没有舱壁的 $\delta GM = 2.195$ m 是准静态基线。那一行的有效 $GM$ 是 $-1.462$ m。液体与船体相位错开时，倾斜可以超过这份损失。隐式处理自由水面的 Casulli–Zanolli 嵌套 Newton 法要的就是这一步。手算不是答案。它是怀疑那段横摇时捏在手里的数字。