越薄反而越硬 — 壳单元的剪切闭锁与MITC绑定

闭锁不是单元太弱，而是因为剪切应变由求导得到。MITC 改为在绑定点上读取该应变再重新插值。

为流固耦合分析写了一段壳单元代码，用悬臂梁做验证。厚度 0.2 m 时，挠度与理论值吻合到小数点后三位。把厚度降到 2 mm，挠度变成理论值的 0.5%。载荷、网格、材料参数都没有动。200 倍的差距从哪来，MITC（Mixed Interpolation of Tensorial Components，张量分量混合插值）改了代码的哪一行才把它抹掉，记在这里。

单元变薄之后，答案僵住了

把细长比滑块往右推，红色完全积分梁的挠度会灭掉，蓝色绑定梁的挠度还在。

slenderness L/t100elements4

Push the slenderness slider right. The blue beam keeps the same shape all the way to L/t = 1000; the red one flattens against the dashed exact curve and its energy bar turns almost entirely red — that red is transverse shear energy that a thin beam is not supposed to have. Adding elements buys the red beam a little back, but the gap reopens as soon as you make it thinner again.

红色梁是完全积分（full integration）单元。把 $L/t$ 提到 100，它在肉眼尺度上几乎不弯曲。蓝色梁只在单元中央一点读取剪切应变，无论细长比多大，形状都不变。右下角能量条中的红色区域就是剪切能量。

这种现象叫剪切闭锁（shear locking）。问题不是单元太弱，恰恰相反：出现了本不该存在的刚度，把结构拽住了。加密网格能缓解，但不能消除。把单元从 4 个增加到 32 个，在 $L/t = 500$ 时仍然只有正确答案的 1.3%。

两种能量之比按细长比的平方拉开

Reissner–Mindlin 壳理论假设厚度方向的纤维（director，法向纤维）变形后仍保持直线，但不必垂直于中面。于是应变能中留下两项。

\Pi = \frac{1}{2}\int_A \left( D\,\kappa^{2} + k\,G\,t\,\gamma^{2} \right) dA, \qquad D = \frac{E\,t^{3}}{12(1-\nu^{2})}

$\kappa$ 是弯曲曲率， $\gamma$ 是横向剪切应变， $D$ 是弯曲刚度， $k = 5/6$ 是剪切修正系数。弯曲刚度与厚度的三次方成正比，剪切刚度与厚度的一次方成正比。看两者之比。

\frac{k\,G\,t}{D} = \frac{6\,k\,(1-\nu)}{t^{2}}

取 $\nu = 0.3$ 、 $k = 5/6$ ，系数是 $3.5/t^{2}$ 。用单元长度无量纲化后，这个值按 $(L/t)^2$ 增长。也就是说，剪切项是压在弯曲项前面的一个巨大罚项。

在连续理论里这不成问题。厚度变小时 $\gamma$ 以同样的速度趋于零，乘积保持有限。这就是 Kirchhoff 极限。问题在于离散单元能否表示 $\gamma = 0$ 。

扫描代码把 $L = 1$ 固定，再代入 $t = 1/(L/t)$ 。 $L/t = 100$ 时 $t = 0.01$ ，取 $\nu = 0.3$ 、 $k = 5/6$ ，则 $k G t / D = 3.5/t^2 = 35000$ 。 $L/t = 500$ 时 $t = 0.002$ ， $t^2$ 缩小到二十五分之一，罚项升到 875000。 $L/t = 2000$ 时 $t = 5\times 10^{-4}$ ，罚项到 $1.4\times 10^{7}$ 。连续理论里 $\gamma$ 随厚度一起变小， $\gamma^2$ 把这个 $1/t^2$ 抵消掉。 $\gamma_h = -\kappa r$ 不是零，抵消就断了。假剪切原样乘上这个系数。

线性单元做不出零剪切应变

降到一维，原因一行就能看清。考虑一个用同样线性形函数插值挠度 $w$ 和截面转角 $\phi$ 的两节点单元。横向剪切应变的定义是

\gamma(r) = \frac{\partial w}{\partial r} - \phi

把纯弯曲状态直接加到节点上。曲率为 $\kappa$ 的弯曲意味着 $w = \tfrac{1}{2}\kappa r^{2}$ 、 $\phi = \kappa r$ 。单元节点在 $r = \pm 1$ ，所以两个节点挠度都等于 $\kappa/2$ 。对两个相等值做线性插值，结果是常数。

w_h(r) = \frac{\kappa}{2}, \qquad \phi_h(r) = \kappa r \;\;\Longrightarrow\;\; \gamma_h(r) = 0 - \kappa r = -\kappa r

转角是一次函数，能被精确再现；挠度是二次的，再现不了。这个不匹配全部流进 $\gamma_h$ 。真实的 $\gamma$ 是零，单元却带着 $-\kappa r$ ，而这个值还要乘上 $3.5/t^2$ 的罚项。在桁架有限元里只有一个轴向应变，不会出现这种配对。闭锁出现在应变有多个、其中一个必须消失的时候。

单元里面有一个点值得盯住。 $\gamma_h = -\kappa r$ 在 $r = 0$ 处恰好为零。错误的值并非均匀地铺满整个单元，而是有一个点是对的。两个高斯点 $\pm 1/\sqrt{3}$ 上，假剪切还剩 $\mp \kappa/\sqrt{3}$ ，代码里的 GAUSS 就是这两点。完全积分在那里攒刚度。绑定只取 $r = 0$ ，因为这道式子只有那里是零。

壳有三套坐标系 — 该在哪里动手

一维梁只要说"在单元中央读取"就够了。曲面壳则要先确定这句话写在哪套坐标系里。壳单元有三套坐标系。自然坐标 $(r, s, t)$ 把单元展平成 $[-1,1]^3$ 立方体，供计算使用。局部坐标是与中面相切的板坐标系，整体坐标 $(x, y, z)$ 则是装配和加载所在的地方。

动应变的地方是自然坐标。用自然坐标的基向量 $\mathbf{g}_i = \partial \mathbf{x}/\partial \xi^i$ 度量的协变（covariant）分量，与单元几何无关地保持同一个物理含义。

\varepsilon_{ij} = \frac{1}{2}\left( \mathbf{g}_i \cdot \frac{\partial \mathbf{u}}{\partial \xi^{j}} + \mathbf{g}_j \cdot \frac{\partial \mathbf{u}}{\partial \xi^{i}} \right)

即使单元被压扁， $\varepsilon_{rt}$ 仍然是" $r$ 方向的线与厚度方向纤维之间的夹角变化"。在整体坐标里做同样的绑定，绑定的对象会随单元形状而改变。

反过来，本构关系和 $B$ 矩阵需要在整体坐标下使用。于是又多一次变换，而这里有个常见陷阱。按工程剪应变的约定，Voigt 记法的第 4、5、6 个分量上乘了 2。把 $B$ 从自然坐标搬到整体坐标时，必须先把 Voigt 向量还原成 3×3 对称张量（剪切项减半），旋转之后再折回去。直接乘一个 6×6 旋转矩阵，剪切行会差 2 倍或 4 倍。它不像闭锁那样显眼，所以能活得更久。

MITC — 在点上读取应变，再重新插值#

解法只有一句话。不要用位移场求导来得到横向剪切应变，而是在指定的绑定点（tying point）上读取，再对这些值做插值。 对 MITC4， $\gamma_r$ 在 $A(0,-1)$ 和 $C(0,+1)$ 两点读取， $\gamma_s$ 在 $B(1,0)$ 和 $D(-1,0)$ 两点读取。

\gamma_r^{\text{MITC}}(r,s) = \frac{1-s}{2}\,\gamma_r(0,-1) + \frac{1+s}{2}\,\gamma_r(0,+1)

右端根本没有出现 $r$ 。前面那个 $-\kappa r$ 关于 $r$ 是一次的，而这个插值里没有地方装这样的项。

它容易和减缩积分（reduced integration）混淆，但两者结果一致只发生在一维这个特例上。减缩积分靠减少积分点让刚度矩阵变软，因此容易招来零能模式（沙漏模式）。MITC 保持积分不变，改的是应变插值空间本身。刚度仍然被精确积分，也不会出现秩亏。

MITC3+ 和 MITC4 把点放在哪里#

四边形的 MITC4 到这两对就结束了。三角形要麻烦得多。既要对称处理三条边又要保持各向同性的绑定布置并不显然，早期的 MITC3 在畸变网格上收敛很差。

MITC3+ 在转角场上增加一个位于单元中心的气泡（bubble）自由度，并把绑定点从三条边的中点向内收。代价是多一个自由度，换来的是在畸变三角网格上也能得到与厚度无关的（uniformly optimal）收敛。在必须用三角形覆盖任意曲面的工程实践中，这个差别很大。

用 Python 数出的闭锁系数#

把两种格式放进同一段代码，把悬臂梁自由端挠度与理论值对比。区别只有剪切规则那一行。

import numpy as np
 
E, NU, KS, B = 210e9, 0.3, 5.0 / 6.0, 1.0        # 弹性模量、泊松比、剪切修正系数、宽度
G = E / (2 * (1 + NU))
GAUSS = (-3 ** -0.5, 3 ** -0.5)
 
 
def strain_operators(h, tied):
    """两节点线性单元的B矩阵。tied=True 时只在 xi=0 一点读取剪切。"""
    Bb = np.array([0.0, -1 / h, 0.0, 1 / h])                      # phi'
    if tied:                                                      # MITC绑定
        rules = [(np.array([-1 / h, -0.5, 1 / h, -0.5]), h)]
    else:                                                         # 两点高斯：精确积分
        rules = [(np.array([-1 / h, -(1 - x) / 2, 1 / h, -(1 + x) / 2]), h / 2)
                 for x in GAUSS]
    return Bb, rules
 
 
def solve_tip(L, t, nel, tied, P=1.0):
    EI, GA = E * B * t ** 3 / 12, KS * G * B * t
    h, ndof = L / nel, 2 * (nel + 1)
    Bb, rules = strain_operators(h, tied)
    Ke = EI * h * np.outer(Bb, Bb) + sum(GA * w * np.outer(Bs, Bs) for Bs, w in rules)
    K = np.zeros((ndof, ndof))
    for e in range(nel):
        idx = [2 * e, 2 * e + 1, 2 * e + 2, 2 * e + 3]
        K[np.ix_(idx, idx)] += Ke
    f = np.zeros(ndof)
    f[-2] = P                                                     # 自由端横向载荷
    u = np.zeros(ndof)
    free = np.arange(2, ndof)                                     # 固定端 w0 = phi0 = 0
    u[free] = np.linalg.solve(K[np.ix_(free, free)], f[free])
    Ub = sum(0.5 * EI * h * (Bb @ u[2 * e:2 * e + 4]) ** 2 for e in range(nel))
    Us = sum(0.5 * GA * w * (Bs @ u[2 * e:2 * e + 4]) ** 2
             for e in range(nel) for Bs, w in rules)
    return u[-2], Us / (Ub + Us)
 
 
def exact_tip(L, t, P=1.0):
    EI, GA = E * B * t ** 3 / 12, KS * G * B * t
    return P * L ** 3 / (3 * EI) + P * L / GA                     # 弯曲 + 剪切
 
 
def sweep_slenderness(nel):
    print(f"nel = {nel:2d}    w_fem / w_exact          shear energy fraction")
    print("  L/t       full        tied            full      tied")
    for ratio in (5, 20, 100, 500, 2000):
        L, t = 1.0, 1.0 / ratio
        ex = exact_tip(L, t)
        wf, ff = solve_tip(L, t, nel, tied=False)
        wm, fm = solve_tip(L, t, nel, tied=True)
        print(f"{ratio:5d}   {wf / ex:10.5f}   {wm / ex:10.5f}"
              f"      {ff:8.4f}   {fm:.4f}")
 
 
sweep_slenderness(4)
print()
sweep_slenderness(32)

nel =  4    w_fem / w_exact          shear energy fraction
  L/t       full        tied            full      tied
    5      0.66631      0.98485        0.3639   0.0307
   20      0.11095      0.98441        0.8910   0.0020
  100      0.00497      0.98438        0.9951   0.0001
  500      0.00020      0.98438        0.9998   0.0000
 2000      0.00001      0.98438        1.0000   0.0000
 
nel = 32    w_fem / w_exact          shear energy fraction
  L/t       full        tied            full      tied
    5      0.99224      0.99976        0.0380   0.0303
   20      0.88873      0.99976        0.1132   0.0019
  100      0.24213      0.99976        0.7579   0.0001
  500      0.01262      0.99976        0.9874   0.0000
 2000      0.00080      0.99976        0.9992   0.0000

把模式从 pure bending 拨到 true shear，红色楔形消失，蓝色曲线贴上真实剪切。

section s0.50half-length / thickness20

Start in pure bending: the mid-surface stays flat while the directors fan out, and every red wedge is shear strain the element invented. Only at r = 0 does the wedge close — that is why the tying point sits there. Switch to true shear and watch the blue MITC curve sit exactly on the red one: tying removes the spurious strain without touching the real one. Then drag the aspect-ratio slider and read the energy numbers.

在 pure bending 模式下，中面是平的，法向纤维却像扇子一样张开。它们之间红色楔形的角度就是单元凭空造出的剪切应变，只有在 $r = 0$ 处楔形闭合。切到 true shear，蓝色 MITC 曲线正好压在红色曲线上。绑定只抹掉假的应变，不碰真的。

nel = 4 那张表横着读，完全积分的挠度比一格一格塌下去。 $L/t = 5$ 时的 0.66631，到 20 是 0.11095，100 是 0.00497，500 是 0.00020，2000 是 0.00001。细长比变成 5 倍的两段里，20 到 100 大约缩到二十二分之一，100 到 500 大约缩到二十五分之一，罚项 $(L/t)^2$ 在同样的跨度上也放大 25 倍。只有 5 到 20 停在大约 6 倍。理论挠度 $P L^3/(3 E I) + P L/G A$ 里的剪切项还活着，闭锁单独的平方规律还没有露出来。

同一行的剪切能量占比是 0.3639、0.8910、0.9951、0.9998、1.0000。 $L/t = 100$ 时能量里已经有 0.9951 是剪切。 $L/t = 2000$ 时是 1.0000。弯曲问题，全靠剪切在撑。绑定那一列从 0.98485 开始，贴在 0.98438 上。缺的 1.5% 把厚度推到 2000 也不会再拉开。那是四个单元的离散误差。

加到 32 个单元，绑定列是 0.99976。完全积分仍然是 $L/t = 100$ 的 0.24213， $L/t = 500$ 的 0.01262， $L/t = 2000$ 的 0.00080。网格加密八倍，细长比 500 时仍只有正确答案的 1.3%。闭锁不是把网格切细就能抹掉的误差。

把壳装进流固耦合代码之前

在 CFD 这一侧，壳通常出现在耦合分析里。薄板或薄膜随流动摆动，其位移再传回网格或浸没边界标记点。这时闭锁会安静地给出错误答案。结构过硬会把固有频率推高，颤振起始速度和附加质量效应也会随之一起偏移。残差照样下降，迭代照样收敛，出错的是刚度矩阵。

所以在接入壳单元之前有三件事要确认。细长比提高十倍时，归一化挠度能否保持？故意把单元弄畸变后，这个值还站得住吗？如果需要大变形，初始应力（几何刚度）项是否也遵循同一套坐标变换规则？前两件用一根悬臂梁半小时就能做完。跳过它们，之后就会到流体那一侧去找原因。