角度0で鉛直の枠は0になり、壁の4自由度を残すと16×16は特異のままです。

橋はなぜあんなにたわむのか — トラス有限要素の直接剛性法

1956年、ボーイングの M. J. Turner は後退翼の応力を手で追って行き詰まりました。部材が数百本あると、つり合いを一本ずつ解けません。同僚と出した答えが直接剛性法(direct stiffness method)です。要素剛性を作り、自由度が重なる位置で足し、 $Ku = f$ を一度解きます。

角度から見ます。棒が水平に寝ると、鉛直の枠が死にます。角度スライダを0まで下げると、4×4剛性の鉛直自由度の行と列がまとめて0になります。

Orientation θ: 30°Auto-rotate

Teal cells are positive, pink negative. At θ = 0° only the horizontal DOFs carry stiffness; rotate the bar and the entries redistribute as c², s², and cs — every 2×2 block is rank one.

$\theta$ が0のとき空になる枠を見てください。棒は自分の軸方向にしか硬くありません。角度を上げると、その0が $x$ と $y$ に分かれます。

剛性は一本のばねから始まる

棒(bar)要素は、軸方向の力しか受けないばねです。長さ $L$ 、断面積 $A$ 、弾性係数 $E$ なら、軸剛性は $k = EA/L$ です。伸びた長さあたりに必要な力です。両端の軸変位 $u_1, u_2$ に対する節点力は次のとおりです。

\begin{pmatrix} f_1 \\ f_2 \end{pmatrix} = \frac{EA}{L} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}

$f_1, f_2$ は両端の節点力、括弧の中が局所剛性行列 $k_e$ です。行の和は0です。二節点が一緒に動く剛体移動には力が要りません。境界条件を外すと、そのために行列が特異になります。

局所から全体へ — 回転変換

部材は四方に傾いています。局所の軸変位を全体の $x, y$ に直さないと組み立てられません。部材が $x$ 軸となす角を $\theta$ 、 $c = \cos\theta$ 、 $s = \sin\theta$ とします。軸変位は全体変位の射影です。 $u_{\text{axial}} = c\,u_x + s\,u_y$ 。この変換を両側から掛けると、4×4の全体要素剛性になります。

k_e = \frac{EA}{L} \begin{pmatrix} c^2 & cs & -c^2 & -cs \\ cs & s^2 & -cs & -s^2 \\ -c^2 & -cs & c^2 & cs \\ -cs & -s^2 & cs & s^2 \end{pmatrix}

各項は節点1の $x, y$ と節点2の $x, y$ をつなぎます。 $\theta = 0$ なら $s = 0$ 、 $c = 1$ なので、上の行列は $EA/L$ を掛けたまま一行目が 1, 0, −1, 0、二行目は全部0です。三行目は一行目の符号を反転したもの、四行目も0です。2×2ブロックがランク1なのは、この0のせいです。座標が傾くと剛性テンソルも一緒に回ります。その変換はメトリックで構成テンソルを回す場合に分けて書いてあります。

重なる自由度を足す — 全体組み立て

組み立てとは、重なる自由度で剛性を足すことです。節点 $i$ の自由度は全体インデックス $2i, 2i+1$ に置きます。 $x$ の次が $y$ です。要素の局所自由度4個をこのインデックスへ送り、その位置に要素剛性を積みます。

K = \sum_{e} \mathbf{L}_e^\top \, k_e \, \mathbf{L}_e

$\mathbf{L}_e$ は要素自由度を全体へ送る選択行列、 $K$ は全体剛性行列です。コードは選択行列を掛けるのではなく、インデックスへ直接足します。一つの節点に部材が何本も付けば、その対角ブロックに寄与が積み上がります。

境界条件：支点を消す

組み立てただけの $K$ は特異行列です。構造がまだ浮いていて、剛体移動が残っているからです。支点で変位を0に縛らないと解けません。自由度を自由集合 $f$ と固定集合 $c$ に分けます。

\begin{pmatrix} K_{ff} & K_{fc} \\ K_{cf} & K_{cc} \end{pmatrix} \begin{pmatrix} u_f \\ u_c \end{pmatrix} = \begin{pmatrix} F_f \\ R_c \end{pmatrix}

$u_c = 0$ なら上の行だけが残り、 $K_{ff}\,u_f = F_f$ です。この縮小系は非特異です。反力はあとから $R_c = K_{cf}\,u_f$ で回収します。 $R_c$ は支点が構造へ返す力です。

荷重スライダとEAを一緒に上げると、先端のたわみがどう縮むかが見えます。

Load P: 24 kNAngle: 0°EA: 8 MNPulse load

Red members are in tension, blue in compression; thickness scales with axial force. Raise EA and the same load bends the truss far less — stiffness is literally the matrix that maps load to displacement.

EAを大きくすると、同じ荷重でもたわみは減ります。剛性行列が、荷重から変位への写像そのものです。荷重を固定してEAだけ二倍にし、先端たわみが半分へ向かうかを見てください。

Python — 12部材のトラスを解く#

壁に固定した3パネルの片持ちです。節点8個、部材12本。座標と接続、荷重を入れて、変位と軸力を出します。

import numpy as np
 
nodes = np.array([[0,0],[0,1],[1,0],[1,1],[2,0],[2,1],[3,0],[3,1]], float)
members = [(0,2),(2,4),(4,6),(1,3),(3,5),(5,7),
           (2,3),(4,5),(6,7),(1,2),(3,4),(5,6)]
EA = 8.0e6           # 軸剛性 EA [N]
fixed = [0, 1]       # 壁に固定された節点
 
def bar_stiffness(p1, p2, EA):
    d = p2 - p1
    L = np.hypot(*d)
    c, s = d / L
    k = EA / L * np.array([[ c*c,  c*s, -c*c, -c*s],
                           [ c*s,  s*s, -c*s, -s*s],
                           [-c*c, -c*s,  c*c,  c*s],
                           [-c*s, -s*s,  c*s,  s*s]])
    return k, L, (c, s)
 
ndof = nodes.shape[0] * 2
K = np.zeros((ndof, ndof))
geom = []
for a, b in members:                       # 全体剛性を組み立てる
    k, L, cs = bar_stiffness(nodes[a], nodes[b], EA)
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    K[np.ix_(dof, dof)] += k
    geom.append((L, cs))
 
F = np.zeros(ndof)                          # 自由端（節点6,7）に24 kNの荷重
for n in (6, 7):
    F[2*n+1] -= 12.0e3
 
fixed_dof = [d for n in fixed for d in (2*n, 2*n+1)]
free_dof = [d for d in range(ndof) if d not in fixed_dof]
 
u = np.zeros(ndof)                          # 縮約系 K_ff u_f = F_f
u[free_dof] = np.linalg.solve(K[np.ix_(free_dof, free_dof)],
                              F[free_dof])
 
for m, (a, b) in enumerate(members):        # 部材軸力 N = (EA/L)[-c,-s,c,s]·u
    L, (c, s) = geom[m]
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    N = EA / L * np.array([-c, -s, c, s]) @ u[dof]
    print(f"member {a}-{b}: N = {N/1e3:+7.2f} kN "
          f"({'tension' if N > 0 else 'compression'})")
 
print(f"free-end drop = {u[2*6+1]*1e3:.3f} mm")

np.ix_ が要素剛性を全体の枠へ足します。このループが商用FEAの組み立てと同じ形です。

節点間隔は1 mです。水平・鉛直の部材9本は $L = 1$ で、コードのEAが $8.0\times 10^6$ Nなら $k = EA/L = 8.0\times 10^6$ N/mです。先端節点一つの鉛直荷重は $12\times 10^3$ Nなので、この棒一本だけを押すと変位は $1.5\times 10^{-3}$ m、つまり1.5 mmです。トラス先端はこの割り算ではありません。部材が荷重を分け合い、実際に解く行列は縮小した $K_{ff}$ です。

斜材は三本です。(0, 1)と(1, 0)をつなぐ部材の長さは $\sqrt{2}$ なので、剛性は水平部材の $1/\sqrt{2}$ 倍です。EAが同じでも斜材の方が柔らかいです。その部材の角度は45°です。 $c^2$ と $s^2$ は $1/2$ で、 $cs$ の絶対値も $1/2$ です。4×4で0でない枠はすべて $EA/L$ の半分、長さが $\sqrt{2}$ なので枠の大きさは $EA/(2\sqrt{2})$ です。水平部材は鉛直の枠が0のまま空き、斜材は四方向がその大きさで埋まります。同じEAが角度だけでこう分かれます。節点8個の自由度は16個です。壁の節点0と1が4個を拘束するので、解く系は12×12です。二次元の最小拘束は3個ですが、この例は4個を拘束します。壁節点一つの $y$ を外すと数は最小にかかり、両方外すと $K_{ff}$ は再び特異になります。

剛性行列が特異になるとき

現場でいちばん見る失敗は特異行列です。分かれ道は三つです。

境界条件が欠けると剛体移動が残り、 $K_{ff}$ が特異です。二次元は自由度3個、三次元は6個を少なくとも拘束します。

三角形になっていない四角形パネルは、部材不足のメカニズムです。トラスは三角形で埋めないと揺れます。

長さが0、または節点が重なると $L = 0$ で除算が破裂します。座標が同じ点は、結合の前に除きます。

壁の節点を一つ外すと、そこで組み立ては止まる

水平棒の1.5 mmは部材一本の話です。12部材片持ちの先端たわみは $K_{ff}\,u_f = F_f$ を解いて初めて出ます。壁の4自由度を外さないうちは、その式が成り立ちません。支持が足りず剛体モードが残ると、変位は数値ではなく零空間です。

棒を面要素に替えた瞬間、同じ組み立てが別の理由で固まります。薄いシェルでその固まりがどう出るかは、せん断ロッキングとMITCが剛性を戻す過程にあります。トラスでは支持4個を先に確認し、厚さを落としたモデルのときだけそちらへ移ります。

橋はなぜあんなにたわむのか — トラス有限要素の直接剛性法

橋はなぜあんなにたわむのか — トラス有限要素の直接剛性法

剛性は一本のばねから始まる

局所から全体へ — 回転変換

重なる自由度を足す — 全体組み立て

境界条件：支点を消す

Python — 12部材のトラスを解く#

剛性行列が特異になるとき

壁の節点を一つ外すと、そこで組み立ては止まる

関連記事