要素を45°ゆがめたら剛性が2.5倍になった — 計量テンソルと構成テンソルの座標変換

正規直交基底では共変成分と反変成分が同じ数字になります。その偶然が崩れた瞬間、構成行列をそのまま使ったコードは嘘をつき始めます。

同じ変形を二度測って答えが割れます。要素はそのままです。材料もそのままです。実際に起きた変形もそのままです。変えたのは、その変形を書き留めた座標系だけです。それでもひずみエネルギーは2.5倍になります。2.5倍が共変・反変基底と計量のどこで生まれるのか、4階の構成テンソルをどう回せばエネルギーが1.412308に固定されるのかを、表の角度で合わせます。

skew を0から45°へ押すか、|g_2| を1から上げると、白い矢印はそのままなのに共変の数字と反変の数字が離れます。

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 | は緑です。角度を少しずらすだけで二行が離れます。矢印自体は動きません。ゆがんだ要素の剛性が狂う場所は、この割れ目です。

変形は動いていないのにエネルギーが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$ は要素の自然座標です。この3本のベクトルは互いに直交せず、長さも1ではありません。

直交しない基底の上で1本のベクトルを成分で書く方法は二つあります。基底に沿って平行四辺形で分解すると、反変成分 $v^i$ が出ます。各基底ベクトルの上に垂直に落とすと、共変成分 $v_i$ です。

\mathbf{v} = v^i \mathbf{g}_i, \qquad v_i = \mathbf{v}\cdot\mathbf{g}_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$ )、それは長さではありません。可視化の最下行がその値を赤で見せています。

連続体力学が応力に反変、ひずみに共変を使う理由がここにあります。仮想仕事はスカラーでなければならず、第2Piola-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$ はクロネッカーのデルタです。つまり $\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}$ を求めるために既に計算してある行列です。反変基底はおまけで出てきます。

4階テンソルは方向余弦を4回掛ける#

ここからが本題です。局所デカルト基底 $\mathbf{e}_p$ で定義された構成テンソル $C^{pqrs}$ を自然座標系へ移します。2階テンソルを移すとき方向余弦を2回掛けたなら、4階テンソルは4回掛けます。

\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}$ が4つ残ります。その4つが表の誤差列です。skew と |g_2| でその列がどう膨らむかは、出力の直後に棒で見ます。

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):
    """局所デカルト座標系での等方平面応力4階テンソル 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から45°へ押すか |g_2| を1から上げると、左の要素の形はそのままで、右の赤いバーだけが伸びます。shear に切り替えると、バーがいちばん急です。

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.

45°の+152.61%と stretch 1.5 の+105.60%が、その二本のバーです。stretch は角度が0でも |g_2| だけで誤差曲線を持ち上げます。

表のエネルギーは角度ごとの同じ引き算です。正しい側の $U$ は1.412308で固定です。座標を45°までひねっても、Voigtに畳んでも、この数字は動きません。変換を省いた側は育ちます。5°で1.486003、誤差は+5.22%です。15°で1.722495、+21.96%です。1.722495を1.412308で割ると約1.22で、表の百分率と揃います。曲面格子で15°はありふれた角度です。長方形だけを回した検証は skew=0 の行の0.00だけ見て通り過ぎます。

45°の行が題の2.5倍です。誤ったエネルギー3.567692を1.412308で割ると約2.53です。差だけでも2.155あります。そのとき $g_{12}$ は+0.707、つまり $\sin 45^\circ$ です。30°では $g_{12} = +0.500 = \sin 30^\circ$ で、誤差は+72.57%、エネルギーは2.437219です。角度が0でも長さが違えば同じ列が膨らみます。 $|\mathbf{g}_2| = 1.5$ なら $g_{22} = 2.25$ で誤差は+105.60%です。2.0なら $g_{22} = 4.00$ 、誤差は+407.84%です。原因はゆがんだ角ではなく、計量が単位行列でないことです。最終行の残差 $9.3\times 10^{-10}$ は、4階の縮約とVoigtの合同変換が同じ1.412308を返す確認にすぎず、変換を省いた列は救いません。

Voigtに畳めば6×6行列一つになる#

添字4つの配列を実務コードにそのまま持ち歩くことはありません。応力・ひずみテンソルは対称なので独立成分は6個だけになり、4階テンソルは $6\times6$ 行列に畳まれます(上のコードは2Dなので $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}$ が、この式が4階テンソルの縮約と同じ答えを出すことの確認です。実際のシェルコードでよく見る T 行列がこれで、せん断補正係数5/6と平面応力仮定( $\sigma_{33}=0$ )は $\mathbf{C}_v$ を作るときに局所デカルト座標系で先に反映します。順序を入れ替えてはいけません。平面応力条件は、板厚方向が定義されているその座標系でしか意味を持たないからです。

有限体積法で同じ間違いが出る場所

これは構造コードだけの間違いではありません。曲線格子の有限体積法で粘性応力テンソルを計算するときも同じ構造が出ます。ひずみ速度を自然座標の微分から得て、ニュートンの粘性 $\tau = 2\mu \mathbf{D}$ をデカルト形のまま掛けると、上の表の誤差列がそのまま戻ります。45°で+152.61%だったあの列です。手元の成分が物理成分か共変か、縮約で上添字と下添字が組になっているかが先です。直交格子だけを通した検証は、この列を見ません。

曲線セルの面積を対角の外積で測る場所でも、Green–Gaussは同じ計量を見ます。シェルのせん断ロッキングとMITCタイイングが要素の定式化を直す場所なら、ここではその定式化をどの座標系で読むかが問題です。定式化を直しても、座標系を外しても、平面パッチテストは通ります。45°が一枚、あるいは $|\mathbf{g}_2| = 1.5$ が一枚、表を割ります。