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

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

同じ変形を二度測って答えが食い違うことがあります。要素はそのまま、材料もそのまま、実際に起きた変形もそのままです。変えたのは、その変形を書き留める座標系だけ。それなのにひずみエネルギーが2.5倍になります。この記事では、その2.5倍がどこから来るのかを共変・反変基底と計量テンソルで押さえ、4階の構成テンソルをどう変換すれば数字が元に戻るのかを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$ は要素の自然座標です。この3本のベクトルは互いに直交せず、長さも1ではありません。

直交しない基底の上で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$ )、それは長さではありません。可視化の最下行がその値を赤で見せています。

連続体力学が応力に反変、ひずみに共変を使う理由がここにあります。仮想仕事はスカラーでなければならず、第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 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):
    """局所デカルト座標系での等方平面応力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 の行では誤差がちょうど0です。長方形要素だけで回した検証は、このバグを絶対に捕まえられません。第二に、15°のゆがみで既に22%です。実際の曲面格子ではありふれた角度です。第三に、直交しているのに $|\mathbf{g}_2|$ を1.5倍にするだけで誤差が105%になります。原因はゆがみではなく、計量が単位行列でないことだという意味です。

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}$ をデカルト形のまま適用すれば、上の表の誤差がそのまま再現します。

三つ確認すれば片が付きます。一つ、いま手元にあるテンソル成分は物理成分か、それとも共変・反変成分か。二つ、縮約するとき上添字と下添字が組になっているか。三つ、検証ケースにゆがんだ要素が一つでもあるか。

実務上いちばん効くのは三つ目です。非直交拡散フラックス補正でもそうだったように、非直交性が生む誤差は直交格子の検証を100%通過した後に現れます。シェル要素のせん断ロッキングとMITCタイイングが要素の定式化を直す話だったとすれば、この記事はその定式化をどの座標系で読むかという話です。どちらを間違えても、平面パッチテストは通ってしまいます。