只转了30°，应变却读出-13% — 大变形壳单元的应力·应变共轭对

大变形下不能把任意应力和任意应变相乘。只有 $S:\dot{E}$、$P:\dot{F}$、$J\sigma:d$ 给出同一个值。

只是转了一下，应变片却指向 -13%#

把一个壳单元在平面内转了30°。没有拉伸，也没有扭转。只有刚体转动。

可是计算工程应变 $\varepsilon_{xx} = \partial u_x / \partial x$ ，结果是 $-0.134$ 。也就是13.4%的压缩。单元哪里都没有变形，应变片却读出了压缩。

这个数不是笔误，也不是离散化误差。 $\cos 30° - 1 = -0.134$ ，正好就是这个值。小应变的定义本身，就把转动误读成了变形。

本文要讲的是：在哪里切断这种误读，切断之后又该把应力换成什么。读完可以带走三样东西。对转动免疫的应变度量、应力与应变配错对时能量偏差的实测值，以及闭合壳单元厚度方向未知量的那个方程。

在下面的模拟中亲手操作一下。

theta30°lam1.00gam0.00

Press rigid rotation only and let it spin: the patch never changes shape, yet the red eps bars swing all the way across while the green E bars stay pinned at zero. At theta = 30° the small-strain gauge reads 0.0000 against 0.0000 — a gap of 0.0000 invented by the rotation alone. Now add lam and gam: E moves, and it keeps the same value at every theta.

点击 rigid rotation only 让转动跑起来。单元形状纹丝不动，红色的 $\varepsilon$ 条却在左右大幅摆动。绿色的 $E$ 条始终贴在0上。把 lam 调大， $E$ 才开始动，而且这个值无论转角怎么变都不改变。

滤掉转动的不是 $F$ ，而是 $F^{T}F$ #

变形的故事从变形梯度(deformation gradient，从参考构形到当前构形的局部映射)开始。

F_{iJ} = \frac{\partial x_i}{\partial X_J}

$x$ 是变形后的坐标， $X$ 是变形前的坐标。用极分解(polar decomposition)可以把 $F$ 拆成 $F = RU$ 。 $R$ 是转动， $U$ 是纯拉伸。问题在于 $F$ 本身还带着 $R$ 。

$\varepsilon = \tfrac{1}{2}(F + F^{T}) - I$ 直接用了 $F$ ，于是 $R$ 就漏了进来。只有转动时 $\varepsilon_{xx} = \cos\theta - 1$ ，原因就在这里。

但把 $F$ 平方一下，转动就消失了。

C = F^{T}F = U^{T}R^{T}RU = U^{T}U

因为 $R^{T}R = I$ ， $C$ 里只剩下 $U$ 。这个 $C$ 就是右柯西–格林张量。减去单位张量再除以二，得到的就是格林–拉格朗日应变。

E = \tfrac{1}{2}\left(F^{T}F - I\right)

刚体转动下 $F^{T}F = I$ ，所以 $E$ 精确为0。上面模拟中绿色条不动的原因，就是这两行。 "换个坐标系物理不该改变"这个要求，在本构张量的坐标变换里是靠基的变换解决的；在这里，则是靠应变度量的定义本身解决。

应力张量是除以哪个面积得到的

应变已经拉回参考构形，应力也得跟着一起拉回来。应力是"力除以面积"，但在大变形里，这个面积是变形前的还是变形后的，答案会分岔。一张表就能理清。

张量	力作用的面	相除的面积	对称性	共轭的应变率
Cauchy $\sigma$	变形后	变形后	对称	$d$ （但要写成 $J\sigma:d$ ）
第一 PK $P$	变形后	变形前	非对称	$\dot{F}$
第二 PK $S$	拉回参考构形	变形前	对称	$\dot{E}$
工程 $\varepsilon$ · $\sigma$	不区分	不区分	对称	仅在小应变极限下

它们之间的关系如下。

P = FS, \qquad \sigma = J^{-1} F S F^{T}, \qquad J = \det F

$J$ 是体积比。第一皮奥拉–基尔霍夫张量 $P$ 之所以非对称，是因为它的两条腿踩在不同的构形上。一个指标指向变形后，另一个指向变形前。所以有限元程序里要存 $P$ ，就得把9个分量全部带着。 $S$ 的两条腿都放在参考构形上，6个分量就够了。

共轭配对的含义是功率相等

"共轭对(conjugate pair)"不是口味问题，而是一个等式：单位参考体积的内部功率必须给出同一个值。

\dot{W} = S : \dot{E} = P : \dot{F} = J\,\sigma : d

其中 $d = \operatorname{sym}(\dot{F}F^{-1})$ 是变形率(rate of deformation)张量。这三种写法是同一个物理量在三个构形上的表达，数值必须严格一致。

反过来， $\sigma : \dot{E}$ 或 $S : d$ 不对应任何物理量。量纲对得上，也算得出来，但那个数不是功率。用这种组合去搭有限元残差，刚度矩阵就不再是能量泛函的海森矩阵。这和伽辽金对称性讲的是同一件事。没有可最小化的能量，牛顿迭代就失去二阶收敛。

用 Python 在同一时刻测了三种配对的功率#

构造一条把拉伸、剪切和转动混在一起的变形路径，在 $t=0.7$ 处测三种组合的数值。材料用圣维南–基尔霍夫模型， $S = \lambda\,\mathrm{tr}(E)I + 2\mu E$ 。

import math
 
I3 = [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]]
 
def mul(A, B):
    return [[sum(A[i][k]*B[k][j] for k in range(3)) for j in range(3)] for i in range(3)]
 
def tr(A):
    return [[A[j][i] for j in range(3)] for i in range(3)]
 
def add(A, B, s=1.0):
    return [[A[i][j] + s*B[i][j] for j in range(3)] for i in range(3)]
 
def scale(A, s):
    return [[s*A[i][j] for j in range(3)] for i in range(3)]
 
def ddot(A, B):
    return sum(A[i][j]*B[i][j] for i in range(3) for j in range(3))
 
def trace(A):
    return A[0][0] + A[1][1] + A[2][2]
 
def det(A):
    return (A[0][0]*(A[1][1]*A[2][2] - A[1][2]*A[2][1])
          - A[0][1]*(A[1][0]*A[2][2] - A[1][2]*A[2][0])
          + A[0][2]*(A[1][0]*A[2][1] - A[1][1]*A[2][0]))
 
def inv(A):
    d = det(A)
    C = [[0.0]*3 for _ in range(3)]
    for i in range(3):
        for j in range(3):
            m = [[A[r][c] for c in range(3) if c != j] for r in range(3) if r != i]
            C[j][i] = ((-1)**(i+j))*(m[0][0]*m[1][1] - m[0][1]*m[1][0])/d
    return C
 
def sym(A):
    return scale(add(A, tr(A)), 0.5)
 
def green_lagrange(F):
    return scale(add(mul(tr(F), F), I3, -1.0), 0.5)
 
def linear_strain(F):
    return sym(add(F, I3, -1.0))
 
LAM, MU = 100.0, 60.0                      # 圣维南–基尔霍夫常数
 
def pk2_stress(E):
    return add(scale(I3, LAM*trace(E)), E, 2*MU)
 
def cauchy_stress(F, S):
    return scale(mul(mul(F, S), tr(F)), 1.0/det(F))
 
def rot_z(th):
    c, s = math.cos(th), math.sin(th)
    return [[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]
 
print("--- 1. pure rotation, no stretch ---")
print(" theta   eps_xx(linear)   E_xx(Green-Lagrange)")
for deg in (0, 5, 10, 30, 60, 90):
    F = rot_z(math.radians(deg))
    print("%5.0f   %14.5f   %20.2e" % (deg, linear_strain(F)[0][0], green_lagrange(F)[0][0]))
 
def defo_path(t):
    """先施加拉伸和剪切，再做 40度 * t 的刚体转动"""
    U = [[1 + 0.20*t, 0.15*t,     0.0],
         [0.15*t,     1 - 0.05*t, 0.0],
         [0.0,        0.0,        1 - 0.08*t]]
    return mul(rot_z(math.radians(40.0)*t), U)
 
def rate(f, t, h=1e-6):
    A, B = f(t + h), f(t - h)
    return [[(A[i][j] - B[i][j])/(2*h) for j in range(3)] for i in range(3)]
 
print()
print("--- 2. work rate at t=0.7, three pairings ---")
t = 0.7
F = defo_path(t)
Fd = rate(defo_path, t)
E = green_lagrange(F)
Ed = rate(lambda s: green_lagrange(defo_path(s)), t)
S = pk2_stress(E)
P = mul(F, S)
J = det(F)
sig = cauchy_stress(F, S)
d = sym(mul(Fd, inv(F)))                   # 变形率张量
 
print("  J = det F              = %.6f" % J)
print("  S : Edot   (2nd PK  x GL rate)   = %12.6f" % ddot(S, Ed))
print("  P : Fdot   (1st PK  x F rate)    = %12.6f" % ddot(P, Fd))
print("  J sigma: d (Cauchy  x stretching)= %12.6f" % (J*ddot(sig, d)))
print("  sigma : Edot   <- wrong pair     = %12.6f" % ddot(sig, Ed))
print("  S : d          <- wrong pair     = %12.6f" % ddot(S, d))

--- 1. pure rotation, no stretch ---
 theta   eps_xx(linear)   E_xx(Green-Lagrange)
    0          0.00000               0.00e+00
    5         -0.00381               0.00e+00
   10         -0.01519              -5.55e-17
   30         -0.13397               0.00e+00
   60         -0.50000               0.00e+00
   90         -1.00000               0.00e+00
 
--- 2. work rate at t=0.7, three pairings ---
  J = det F              = 1.028087
  S : Edot   (2nd PK  x GL rate)   =    10.522306
  P : Fdot   (1st PK  x F rate)    =    10.522306
  J sigma: d (Cauchy  x stretching)=    10.522306
  sigma : Edot   <- wrong pair     =     9.961790
  S : d          <- wrong pair     =     4.823514

5% 和 54% — 配错共轭对会产生的两类误差#

三种正确组合小数点后六位完全一致，都是 $10.522306$ 。同一个量只是分写在三个构形上，这一点由数值直接证实。

两种错误组合的错法却不一样。 $\sigma : \dot{E}$ 得到 $9.961790$ ，低了5.3%。 $\sigma$ 和 $S$ 相差一个 $J^{-1}F(\cdot)F^{T}$ ，而这个变换本身不大。 $J = 1.028$ ，拉伸也就20%出头，所以误差止步于这个量级。

$S : d$ 得到 $4.823514$ ，低了54%。这不是尺度问题，而是另一类错误。它无视了 $\dot{E} = F^{T} d F$ 这个关系，直接把 $d$ 塞了进去，于是转动分量也混了进来。转角越大，这个误差越大。

工程上更危险的是那个5%。54%在第一个载荷步就发散，马上会被抓住。5%却能收敛，只是收敛到一个略微错误的答案。加密网格也消不掉。

厚度方向由体积闭合，而不是靠第三个方程

到这里为止都是一般连续体的内容。壳还要再加一条。

三维壳单元把厚度方向的伸长当作未知量。在 Sussman 和 Bathe 的三维壳格式中，厚度方向需要3个未知量，而平面应力条件只能给出2个方程。少了一个。

补上这一个的，是不可压缩条件。

J = \lambda_1 \lambda_2 \lambda_3 = 1 \quad\Longrightarrow\quad \lambda_3 = \frac{1}{\lambda_1 \lambda_2}

$\lambda_i$ 是主伸长比。在橡胶或金属塑性这类体积几乎守恒的材料里，厚度不是独立未知量，而成了面内伸长的因变量。面内拉伸 $\lambda_1 = 1.20$ 、 $\lambda_2 = 1.05$ 时，厚度变为 $0.794$ 倍，减薄20.6%。板料成形中计算减薄率，用的正是这个关系。

这个条件和 MITC 绑定有交汇之处。如果在单元内部直接插值厚度方向伸长，薄单元上会出现体积锁死(volumetric locking)。就像用绑定点解决剪切锁死那样，厚度伸长也只在每个单元的少数几个点上独立取值，其余用插值绑住。

把转动截断到一阶，director 会长 16%#

壳的第二条是转动。壳单元随身带着中面法向量，也就是 director。牛顿迭代给出转动增量 $\Delta\boldsymbol{\theta}$ 后，就要把 director 转过相应的角度。转动矩阵由罗德里格斯(Rodrigues)公式生成。

R(\boldsymbol{\theta}) = I + \frac{\sin\theta}{\theta}\,\Theta + \frac{1 - \cos\theta}{\theta^{2}}\,\Theta^{2}

$\Theta$ 是 $\boldsymbol{\theta}$ 的反对称矩阵， $\theta = |\boldsymbol{\theta}|$ 。不少代码觉得增量很小，就截断成 $R \approx I + \Theta$ 。这个矩阵并不正交。 $\det(I + \Theta) = 1 + \theta^{2} \neq 1$ ，所以每一步 director 都会稍微变长一点。

在下面的模拟中亲手操作一下。

d.theta0.100

Leave d.theta at 0.10 rad and let it march: after 0 increments the red 1st-order director is 0.0 % too long and has climbed off the dashed unit circle, while the green closed-form arrow sits on it. Drag d.theta down — the drift shrinks in proportion, so halving the load step only halves the error. The yellow 2nd-order curve (0.00 %) shows what one more term buys.

把 d.theta 设为0.10再让它推进，红色的一阶截断箭头会螺旋着甩出虚线单位圆之外。绿色的闭式解则精确地留在圆上。把 d.theta 减半，漂移也减半。也就是说，这个误差对增量大小是一阶的。

import math
 
def matvec(A, v):
    return [sum(A[i][k]*v[k] for k in range(3)) for i in range(3)]
 
def matmul(A, B):
    return [[sum(A[i][k]*B[k][j] for k in range(3)) for j in range(3)] for i in range(3)]
 
def skew(w):
    return [[0.0, -w[2], w[1]],
            [w[2], 0.0, -w[0]],
            [-w[1], w[0], 0.0]]
 
def rodrigues(w, order):
    """order = 1, 2 为级数截断，0 为闭式"""
    th = math.sqrt(sum(c*c for c in w))
    W = skew(w)
    W2 = matmul(W, W)
    if order == 1:
        a, b = 1.0, 0.0
    elif order == 2:
        a, b = 1.0, 0.5
    else:
        a = math.sin(th)/th
        b = (1.0 - math.cos(th))/(th*th)
    return [[(1.0 if i == j else 0.0) + a*W[i][j] + b*W2[i][j]
             for j in range(3)] for i in range(3)]
 
def spin_director(dth, steps, order):
    """把中面法向绕 y 轴每步转动 dth 弧度"""
    d = [0.0, 0.0, 1.0]
    for _ in range(steps):
        d = matvec(rodrigues([0.0, dth, 0.0], order), d)
    return d
 
print("--- 3. director after 30 increments of 0.10 rad (exact total 171.89 deg) ---")
print(" order        |d|      length err %   angle(deg)   angle err(deg)")
for order, name in ((1, "1st"), (2, "2nd"), (0, "closed")):
    d = spin_director(0.10, 30, order)
    n = math.sqrt(sum(c*c for c in d))
    ang = math.degrees(math.atan2(d[0], d[2]))
    if ang < 0:
        ang += 360.0
    print(" %-6s  %10.5f   %11.2f   %10.3f   %12.3f"
          % (name, n, 100*(n - 1.0), ang, ang - math.degrees(3.0)))
 
print()
print("--- 4. same total rotation, smaller increments (1st order) ---")
print(" steps   dtheta      |d|     length err %")
for steps in (30, 60, 150, 300, 3000):
    dth = 3.0/steps
    d = spin_director(dth, steps, 1)
    n = math.sqrt(sum(c*c for c in d))
    print(" %5d   %6.4f  %8.5f   %11.3f" % (steps, dth, n, 100*(n - 1.0)))
 
print()
print("--- 5. thickness closed by J = 1, not by a 3rd equation ---")
print(" lam1   lam2    lam3=1/(lam1 lam2)   thickness change %")
for l1, l2 in ((1.20, 1.05), (1.20, 1.00), (1.10, 1.10), (1.30, 0.95)):
    l3 = 1.0/(l1*l2)
    print(" %4.2f   %4.2f   %16.5f   %16.1f" % (l1, l2, l3, 100*(l3 - 1.0)))

--- 3. director after 30 increments of 0.10 rad (exact total 171.89 deg) ---
 order        |d|      length err %   angle(deg)   angle err(deg)
 1st        1.16097         16.10      171.318         -0.570
 2nd        1.00038          0.04      172.173          0.286
 closed     1.00000          0.00      171.887          0.000
 
--- 4. same total rotation, smaller increments (1st order) ---
 steps   dtheta      |d|     length err %
    30   0.1000   1.16097        16.097
    60   0.0500   1.07778         7.778
   150   0.0200   1.03045         3.045
   300   0.0100   1.01511         1.511
  3000   0.0010   1.00150         0.150
 
--- 5. thickness closed by J = 1, not by a 3rd equation ---
 lam1   lam2    lam3=1/(lam1 lam2)   thickness change %
 1.20   1.05            0.79365              -20.6
 1.20   1.00            0.83333              -16.7
 1.10   1.10            0.82645              -17.4
 1.30   0.95            0.80972              -19.0

一阶截断只用30步就把 director 拉长了16.1%。多加一项的二阶降到0.04%。长度精度提高了400倍，角度误差反而是二阶的更大。这两种误差彼此独立。

第4张表更重要。把增量缩到十分之一，误差也只降到十分之一。把载荷步切得更细，解决不了这个问题。要么用闭式解，要么每步对 director 做归一化，要么把转动改用四元数来携带。

如何从代码里分辨是 $S$ 还是 $\sigma$ #

打开别人的大变形代码时，光看变量名认不出应力的身份。看三个地方就能分辨。

看应力进入刚度矩阵的位置。如果 $B$ 矩阵是按 $\partial E / \partial u$ 构造的，与之相乘的应力就是 $S$ 。如果 $B$ 是 $\partial \varepsilon / \partial u$ ，那就是 $\sigma$ 。正如本构矩阵的变换中所见， $B$ 的定义决定了应力的身份。

看积分时的雅可比。如果像 $\int_{V_0} \cdots \, dV_0$ 那样在参考构形体积上积分，而且没有乘 $J$ ，那就是 $S$ 。如果乘了 $J$ ，说明正在把 $\sigma$ 拉回参考构形。

看输出例程。如果后处理在算 von Mises 之前有 $\sigma = J^{-1}FSF^{T}$ 这一步变换，说明内部一直是用 $S$ 在跑。也有代码省掉这步变换，直接把 $S$ 的分量塞进 von Mises 去画。小变形下看不出破绽，一旦伸长超过20%，两者就分道扬镳。

三处都查不出来，还剩下一个测试。让一个单元只做刚体转动，看应力是否保持为0。如果30°时冒出13%，就说明某个地方没有把 $F$ 平方。