薄くするほど硬くなる — シェル要素のせん断ロッキングとMITCタイイング

ロッキングは要素が弱いからではなく、せん断ひずみを微分で得ているために起きます。MITCはその値をタイイング点で読み直します。

流体構造連成解析に使うシェルコードを書き、片持ち梁で検証しました。板厚 0.2 m では理論値と小数点以下3桁まで一致します。板厚を 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まで上げると、目に見える範囲ではまったく曲がりません。青い梁はせん断ひずみを要素中央の1点だけで読んだもので、細長比をいくら上げても形が変わりません。右下のエネルギー棒の赤い領域がせん断エネルギーです。

この現象をせん断ロッキング(shear locking)と呼びます。要素が弱いために起きる問題ではありません。逆に、あってはならない剛性が生じて構造を押さえつけているのです。格子を細かくすれば緩和しますが、消えはしません。要素を4個から32個に増やしても、 $L/t = 500$ では依然として正解の1.3%です。

2つのエネルギーの比は細長比の2乗で開く#

Reissner–Mindlin シェル理論は、板厚方向のファイバー(ディレクター)が変形後も直線であるが中立面に垂直である必要はない、と置きます。その結果、ひずみエネルギーに項が2つ残ります。

\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$ はせん断補正係数です。曲げ剛性は板厚の3乗、せん断剛性は板厚の1乗に比例します。両者の比は次のとおりです。

\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$ は25分の1になり、ペナルティは875000まで上がります。 $L/t = 2000$ なら $t = 5\times 10^{-4}$ で、ペナルティは $1.4\times 10^{7}$ まで行きます。連続体では $\gamma$ が板厚に合わせて減り、 $\gamma^2$ がこの $1/t^2$ を打ち消します。 $\gamma_h = -\kappa r$ がゼロでなければ打ち消しは切れます。偽のせん断に、その係数がそのまま掛かります。

線形要素はせん断ひずみゼロを作れない

1次元に落とすと原因が一行で見えます。たわみ $w$ と断面回転 $\phi$ を同じ線形形状関数で補間する2節点要素を考えます。横せん断ひずみの定義は次のとおりです。

\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

回転は1次関数なので正確に再現されますが、たわみは2次なので再現されません。その不一致がすべて $\gamma_h$ に流れ込みます。真の $\gamma$ はゼロなのに要素は $-\kappa r$ を持っており、その値に $3.5/t^2$ のペナルティが掛かります。トラス有限要素では軸ひずみ1つだけなので、このような組み合わせは生じません。ロッキングは、ひずみが複数あってそのうち1つが消えなければならないときに現れます。

要素の中で目が留まる点が1つあります。 $\gamma_h = -\kappa r$ は $r = 0$ でちょうどゼロです。誤った値が要素全体に一様に広がっているのではなく、1点だけは正しいのです。ガウス点 $\pm 1/\sqrt{3}$ では偽のせん断が $\mp \kappa/\sqrt{3}$ として残り、コードの GAUSS がまさにその2点です。完全積分はそこで剛性を作ります。タイイングが $r = 0$ だけを拾う理由は、この式にあります。

シェルには座標系が3つある — どこで手を入れるか#

1次元の梁なら「要素中央で読む」で終わりです。曲面シェルでは、その文をどの座標系で書くかから決めなければなりません。シェル要素には座標系が3つあります。自然座標 $(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)$ の2点で、 $\gamma_s$ を $B(1,0)$ と $D(-1,0)$ の2点で読みます。

\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$ について1次でしたが、この補間にはそのような項を入れる場所がありません。

低減積分(reduced integration)と混同しやすいのですが、結果が一致するのは1次元の特殊な場合だけです。低減積分は積分点を減らして剛性行列を柔らかくする操作なので、ゼロエネルギーモード(アワーグラス)を招きやすくなります。MITC は積分をそのままにして、ひずみの補間空間そのものを変えます。剛性は依然として正確に積分され、ランク落ちも生じません。

MITC3+ と MITC4 が点を置く位置#

四角形要素の MITC4 は、先の2組で終わりです。三角形は事情が悪くなります。3辺を対称に扱いながら等方性を保つタイイング配置が自明ではなく、初期の MITC3 は歪んだ格子で収束が悪かったのです。

MITC3+ は回転場に要素中央のバブル(bubble)自由度を1つ加え、タイイング点を3辺の中点ではなく内側に引いた3点に置きます。自由度が1つ増える代わりに、歪んだ三角形格子でも板厚に依存しない(uniformly optimal)収束が得られます。任意形状の曲面を三角形で覆う必要のある実務では、この差は大きいです。

Python で数えたロッキング係数#

2つの定式化を同じコードに入れ、片持ち梁の自由端たわみを理論値と比べました。違いはせん断則の一行だけです。

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):
    """2節点線形要素のB行列。tied=True ならせん断を xi=0 の1点だけで読む。"""
    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:                                                         # 2点ガウス:正確積分
        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 は約22倍、100 から 500 は約25倍に減り、ペナルティ $(L/t)^2$ もその区間で25倍になります。5 から 20 だけは約6倍で止まります。理論たわみ $P L^3/(3 E I) + P L/G A$ のせん断項がまだ残っていて、ロッキングだけの2乗則がむき出しにならないからです。

同じ行のせん断エネルギー比は 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 まで上げても広がりません。要素4個の離散化誤差です。

32個に増やすとタイイング列は 0.99976 です。そのときも完全積分は $L/t = 100$ で 0.24213、 $L/t = 500$ で 0.01262、 $L/t = 2000$ で 0.00080 です。格子を8倍細かくしても、細長比 500 では正解の1.3%です。ロッキングは格子を細かくして消せる誤差ではありません。

流体構造連成コードにシェルを載せる前に

CFD 側でシェルを扱う場面は、たいてい連成解析です。薄い板や膜が流れに乗って揺れ、その変位が格子や仮想境界マーカーに渡ります。このときロッキングは静かに誤った答えを返します。構造が必要以上に硬くなると固有振動数が上にずれ、フラッター開始速度や付加質量効果もすべて一緒に狂います。残差はきれいに落ち、反復も収束します。間違っているのは剛性行列です。

ですからシェル要素を組み込む前に確認することは3つです。細長比を10倍にしたとき正規化たわみが保たれるか、要素を意図的に歪ませてもその値が持ちこたえるか、そして大変形が必要なら初期応力(幾何剛性)項まで同じ座標変換規則に従っているか。前の2つは片持ち梁1本で30分あれば終わります。それを飛ばすと、後で流体側から原因を探すことになります。