把单元扭斜45°，刚度就涨了2.5倍 — 度量张量与本构张量的坐标变换

在标准正交基下，协变分量和逆变分量是同一组数字。这个巧合一旦破裂，直接沿用笛卡尔本构矩阵的代码就开始说谎。

同一个变形量两次，答案却对不上。单元没变，材料没变，实际发生的变形也没变。换掉的只是记录这个变形的坐标系。可是应变能却涨到了2.5倍。本文追查这个2.5倍从哪里来，用协变基、逆变基和度量张量把它讲清楚，再用Python验证四阶本构张量该怎么变换才能让数字回到原位。

变形没动，能量却涨了2.5倍#

壳单元的刚度矩阵通常这样组装：在高斯积分点上构造应变-位移矩阵 $\mathbf{B}$ ，乘以本构矩阵 $\mathbf{C}$ ，再积分 $\mathbf{B}^{\mathsf T}\mathbf{C}\mathbf{B}$ 。问题在于这两个矩阵诞生于不同的坐标系。

$\mathbf{C}$ 是材料属性，因此定义在材料试验所用的坐标系里，也就是贴在壳面上的局部正交笛卡尔坐标系。而 $\mathbf{B}$ 来自形函数的导数。形函数写在单元的自然坐标系 $(r^1, r^2, r^3)$ 中，单元一旦弯曲或畸变，这个坐标系既不正交也不归一。

有一种情况两个坐标系恰好重合：单元是矩形且曲面是平的。如果只用平面模型跑分片试验，那就只是在这个巧合上做验证。一旦把单元贴到曲面上，巧合就消失了。

基不正交，分量就有两套

自然坐标 $r^i$ 变化时物理位置 $\mathbf{x}$ 移动的方向，就是协变基。

\mathbf{g}_i = \frac{\partial \mathbf{x}}{\partial r^i}

$\mathbf{x}$ 是笛卡尔位置矢量， $r^i$ 是单元的自然坐标。这三个矢量彼此不正交，长度也不是1。

在不正交的基上写一个矢量的分量有两种做法。第一，沿基矢量按平行四边形法则分解，得到的系数是逆变分量 $v^i$ 。第二，向每个基矢量作垂线投影，那些投影是协变分量 $v_i$ 。

\mathbf{v} = v^i \mathbf{g}_i, \qquad v_i = \mathbf{v}\cdot\mathbf{g}_i

在标准正交基下这两种作图落在同一点。所以只用过笛卡尔坐标的人根本看不见这个区别。在下面的模拟里亲手操作一下。

skew of g_2 [deg]30| g_2 |1.40

v angle [deg]35

The white arrow is one physical vector. Blue dashes slide it along the green base vectors (contravariant vⁱ); yellow dashes drop it perpendicular onto them (covariant v_i). Push the skew slider to 0 and |g₂| to 1.00: the two rows of numbers merge and max | v^i - v_i | turns green. Anywhere else they differ, and only the metric row keeps the length right.

把 skew 拨到0、|g_2| 拨到1.00，蓝色平行四边形和黄色垂线交于同一点，max | v^i - v_i | 变成绿色。滑块只要稍微一动，两行数字就分开。操作时要盯住的关键是：白色箭头始终是同一个矢量。

度量张量把丢掉的长度还回来

连接这两套分量的就是度量张量。

g_{ij} = \mathbf{g}_i \cdot \mathbf{g}_j, \qquad v_i = g_{ij}\, v^j

$g_{ij}$ 把基矢量的长度和夹角装进一个矩阵。对角元 $g_{11}, g_{22}, g_{33}$ 是各基的长度平方，非对角元是夹角余弦乘以长度。

要量长度，就必须经过这个矩阵。

|\mathbf{v}|^2 = g_{ij}\, v^i v^j = v_i v^i

第二个等号才是要点。把协变分量和逆变分量配对相乘，度量会自动抵消。反过来，把逆变分量各自平方再相加( $v^i v^i$ )，那不是长度。可视化最底下一行用红色显示了这个值。

连续介质力学给应力用逆变、给应变用协变，原因就在这里。虚功必须是标量，只有把第二类Piola-Kirchhoff应力的逆变分量 $S^{ij}$ 与Green-Lagrange应变的协变分量 $E_{ij}$ 配对， $S^{ij}E_{ij}$ 才与坐标系无关。这和网格度量在流动求解器里干的事一样 — 曲线坐标变换与网格度量里讲的雅可比矩阵，在这里就是基矢量本身。

逆变基就是雅可比逆矩阵的行

要把协变分量转回逆变，需要 $g^{ij} = [g_{ij}]^{-1}$ ，用这个矩阵构造的基就是逆变基。

\mathbf{g}^i = g^{ij}\mathbf{g}_j, \qquad \mathbf{g}^i \cdot \mathbf{g}_j = \delta^i_j

$\delta^i_j$ 是克罗内克delta。也就是说 $\mathbf{g}^1$ 同时垂直于 $\mathbf{g}_2$ 和 $\mathbf{g}_3$ ，长度则调到与 $\mathbf{g}_1$ 的点积恰好为1。

实现时不需要求两次逆。把协变基按列堆成矩阵 $\mathbf{J}$ ，则

\mathbf{J} = [\,\mathbf{g}_1\;\mathbf{g}_2\;\mathbf{g}_3\,], \qquad \mathbf{g}^i = \big(\mathbf{J}^{-1}\big)_{i\bullet}

$\big(\mathbf{J}^{-1}\big)_{i\bullet}$ 表示 $\mathbf{J}^{-1}$ 的第 $i$ 行。

$\mathbf{J}$ 正是从自然坐标到物理坐标的雅可比矩阵。为了在每个积分点求 $\det\mathbf{J}$ ，这个矩阵早就算好了。逆变基是顺带得到的。

四阶张量要乘四次方向余弦

现在进入正题。把定义在局部笛卡尔基 $\mathbf{e}_p$ 上的本构张量 $C^{pqrs}$ 搬到自然坐标系。搬二阶张量要乘两次方向余弦，四阶张量就乘四次。

\tilde{C}^{ijkl} = (\mathbf{g}^i\!\cdot\mathbf{e}_p)(\mathbf{g}^j\!\cdot\mathbf{e}_q)(\mathbf{g}^k\!\cdot\mathbf{e}_r)(\mathbf{g}^l\!\cdot\mathbf{e}_s)\, C^{pqrs}

指标 $i,j,k,l$ 在自然坐标里， $p,q,r,s$ 在局部笛卡尔坐标里。应变那一侧方向相反。

\tilde{\varepsilon}_{ij} = (\mathbf{g}_i\!\cdot\mathbf{e}_p)(\mathbf{g}_j\!\cdot\mathbf{e}_q)\,\varepsilon_{pq}

本构张量挂上 $\mathbf{J}^{-1}$ ，应变挂上 $\mathbf{J}$ 。所以两者缩并时雅可比矩阵正好抵消，能量保持不变。反过来说，只变换应变而放着本构张量不管，就会剩下四个 $\mathbf{J}$ 。那四个就是下面要测的误差的真身。

skew of g_2 [deg]30| g_2 |1.00

The blue element on the left never changes its physics — only the coordinate lines you measure it with do. Drag the skew slider away from 0, or push |g₂| off 1.00, and watch the red bar climb while the green one stays put. Bring both back to the orthonormal setting and the two bars lock together: that is the only configuration where skipping the transform is free.

拖动 skew 和 |g_2|，左边单元的变形形状纹丝不动，右边只有红色柱子在长。切换 stretch/shear/mixed，看上方误差曲线的形状怎么变，这是观察要点。剪切模式下误差涨得最快。

用Python测出的各扭斜角度下的能量#

在平面应力各向同性材料上固定一个应变状态，只扭动坐标系，用两种方式计算应变能。材料常数与原始文档一致( $E = 2.1\times10^6$ ， $\nu = 0.3$ )。

import numpy as np
 
def natural_basis(skew_deg, stretch=1.0):
    """列为协变基 g_i = dx/dr^i 的雅可比矩阵 J"""
    a = np.deg2rad(skew_deg)
    g1 = np.array([1.0, 0.0])
    g2 = stretch * np.array([np.sin(a), np.cos(a)])
    return np.column_stack([g1, g2])
 
def plane_stress_tensor(E=2.1e6, nu=0.3):
    """局部笛卡尔坐标系下的各向同性平面应力四阶张量 C^{pqrs}"""
    lam = E * nu / (1.0 - nu**2)
    mu = E / (2.0 * (1.0 + nu))
    d = np.eye(2)
    return (lam * np.einsum('pq,rs->pqrs', d, d)
            + mu * (np.einsum('pr,qs->pqrs', d, d) + np.einsum('ps,qr->pqrs', d, d)))
 
def rotate_fourth_order(C, Jinv):
    """C^{ijkl} = (g^i.e_p)(g^j.e_q)(g^k.e_r)(g^l.e_s) C^{pqrs}, g^i = J^-1 的第i行"""
    return np.einsum('ip,jq,kr,ls,pqrs->ijkl', Jinv, Jinv, Jinv, Jinv, C)
 
def strain_energy(C, eps):
    return 0.5 * np.einsum('pqrs,pq,rs->', C, eps, eps)
 
def to_voigt2d(C):
    """2D Voigt: (00,11,01) -> 3x3"""
    idx = [(0, 0), (1, 1), (0, 1)]
    return np.array([[C[p, q, r, s] for (r, s) in idx] for (p, q) in idx])
 
# 一个物理应变状态(局部笛卡尔分量)。无论坐标怎么取，这个张量都不变。
eps_cart = np.array([[1.0e-3, 4.0e-4],
                     [4.0e-4, -6.0e-4]])
C_cart = plane_stress_tensor()
U_ref = strain_energy(C_cart, eps_cart)
 
print("skew   g11     g12     g22    | U_correct     U_naive      err%")
print("-" * 68)
for skew in [0, 5, 10, 15, 20, 30, 40, 45]:
    J = natural_basis(skew)
    g = J.T @ J                       # 度量张量 g_ij
    Jinv = np.linalg.inv(J)           # 行 = 逆变基 g^i
    eps_nat = J.T @ eps_cart @ J      # 协变应变分量
    C_nat = rotate_fourth_order(C_cart, Jinv)
    U_ok = strain_energy(C_nat, eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)   # 漏掉变换的代码
    err = 100.0 * (U_bad - U_ok) / U_ok
    print(f"{skew:3d}   {g[0,0]:.3f}  {g[0,1]:+.3f}  {g[1,1]:.3f} |"
          f" {U_ok:.6e}  {U_bad:.6e}  {err:+8.2f}")
 
print()
print("orthogonal (skew=0), only |g2| stretched")
for st in [1.0, 1.5, 2.0]:
    J = natural_basis(0, stretch=st)
    Jinv = np.linalg.inv(J)
    eps_nat = J.T @ eps_cart @ J
    U_ok = strain_energy(rotate_fourth_order(C_cart, Jinv), eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)
    print(f"  stretch={st:.1f}  g22={(J.T@J)[1,1]:.2f}  err% = {100*(U_bad-U_ok)/U_ok:+9.2f}")
 
print()
print(f"Cartesian reference U_ref   = {U_ref:.6e}")
J = natural_basis(30)
Jinv = np.linalg.inv(J)
C_nat = rotate_fourth_order(C_cart, Jinv)
eps_nat = J.T @ eps_cart @ J
print(f"skew=30 after transform     = {strain_energy(C_nat, eps_nat):.6e}  (invariant)")
 
# 折成Voigt形式后是否得到同样的值
Cv = to_voigt2d(C_nat)
ev = np.array([eps_nat[0, 0], eps_nat[1, 1], 2.0 * eps_nat[0, 1]])
print(f"skew=30 via Voigt 3x3       = {0.5 * ev @ Cv @ ev:.6e}")
 
# Voigt空间中的应变变换矩阵 A: e_v(nat) = A e_v(cart)
def voigt_map(J):
    cols = []
    for e in (np.array([[1, 0], [0, 0]]), np.array([[0, 0], [0, 1]]), np.array([[0, .5], [.5, 0]])):
        n = J.T @ e @ J
        cols.append([n[0, 0], n[1, 1], 2 * n[0, 1]])
    return np.array(cols).T
 
A = voigt_map(J)
Cv_cart = to_voigt2d(C_cart)
Ai = np.linalg.inv(A)
print("Voigt congruence C_nat = A^-T C_cart A^-1  residual =",
      f"{np.max(np.abs(Ai.T @ Cv_cart @ Ai - Cv)):.3e}")

skew   g11     g12     g22    | U_correct     U_naive      err%
--------------------------------------------------------------------
  0   1.000  +0.000  1.000 | 1.412308e+00  1.412308e+00     +0.00
  5   1.000  +0.087  1.000 | 1.412308e+00  1.486003e+00     +5.22
 10   1.000  +0.174  1.000 | 1.412308e+00  1.585621e+00    +12.27
 15   1.000  +0.259  1.000 | 1.412308e+00  1.722495e+00    +21.96
 20   1.000  +0.342  1.000 | 1.412308e+00  1.906550e+00    +35.00
 30   1.000  +0.500  1.000 | 1.412308e+00  2.437219e+00    +72.57
 40   1.000  +0.643  1.000 | 1.412308e+00  3.163176e+00   +123.97
 45   1.000  +0.707  1.000 | 1.412308e+00  3.567692e+00   +152.61
 
orthogonal (skew=0), only |g2| stretched
  stretch=1.0  g22=1.00  err% =     +0.00
  stretch=1.5  g22=2.25  err% =   +105.60
  stretch=2.0  g22=4.00  err% =   +407.84
 
Cartesian reference U_ref   = 1.412308e+00
skew=30 after transform     = 1.412308e+00  (invariant)
skew=30 via Voigt 3x3       = 1.412308e+00
Voigt congruence C_nat = A^-T C_cart A^-1  residual = 9.313e-10

有三点值得读。

第一，skew=0 那一行误差恰好为0。只用矩形单元搭起来的验证套件永远抓不到这个bug。第二，扭斜15°时误差已经是22%，这在真实曲面网格上是再普通不过的角度。第三，完全不扭斜，只把 $|\mathbf{g}_2|$ 拉到1.5倍，误差就是105%。这说明祸首不是畸变，而是度量不是单位矩阵。

折成Voigt形式后就是一个6×6矩阵#

没人会把四指标数组原样带进生产代码。应力和应变张量对称，独立分量只剩6个，四阶张量折成 $6\times6$ 矩阵(上面的代码是二维的，所以是 $3\times3$ )。

折叠之后变换依然在。如果应变的Voigt矢量按 $\tilde{\boldsymbol\varepsilon}_v = \mathbf{A}\,\boldsymbol\varepsilon_v$ 变换，那么能量不变就要求本构矩阵接受一次合同变换。

\tilde{\mathbf{C}}_v = \mathbf{A}^{-\mathsf T}\,\mathbf{C}_v\,\mathbf{A}^{-1}

$\mathbf{A}$ 的元素是方向余弦的乘积。上面输出最后一行的残差 $9.3\times10^{-10}$ 证实了这条路径与四阶张量缩并给出同一个答案。真实壳代码里常见的 T 矩阵就是它。剪切修正系数5/6和平面应力假设( $\sigma_{33}=0$ )要在构造 $\mathbf{C}_v$ 时先在局部笛卡尔坐标系里体现，顺序不能颠倒 — 平面应力条件只在定义了厚度方向的那个坐标系里才有意义。

有限体积法里同样的错误出现在哪

这不是结构代码独有的错误。曲线网格有限体积法计算粘性应力张量时会出现同样的结构。先用自然坐标的导数得到应变率张量，再把牛顿粘性定律 $\tau = 2\mu \mathbf{D}$ 按笛卡尔形式直接套上去，上表里的误差就会原封不动地复现。

确认三件事就够了。其一，手上的张量分量是物理分量，还是协变/逆变分量。其二，缩并时上指标和下指标是否成对。其三，验证算例里有没有哪怕一个畸变单元。

实践中第三点最要紧。正如非正交扩散通量修正里那样，非正交性造成的误差总是在正交网格测试100%通过之后才冒出来。如果说壳单元的剪切闭锁与MITC绑定讲的是修正单元的格式，那么本文讲的是在哪个坐标系里读这个格式。两者任错其一，平面分片试验照样通过。