橋梁やロケットの機体、自動車のフレーム――これらの構造物にかかる力と変形を正確に知りたいとき、私たちはどうすればよいのでしょうか。単純な棒1本であれば、応力とひずみの関係とフックの法則を使って手計算で解くことができます。しかし現実の構造物は、形状が複雑で、荷重や拘束条件も一様ではありません。このような問題に対して、解析解（紙と鉛筆で得られる厳密解）を求めることは、ほとんどの場合不可能です。
この「解析解が得られない問題を、コンピュータの力で系統的に解く」ための最も強力な手法が有限要素法（Finite Element Method, FEM）です。有限要素法は構造解析の世界を根本的に変えた手法であり、以下のような幅広い分野で日常的に使われています。
- 構造解析: 航空機の翼やロケットの胴体が荷重下でどのように変形し、どこに応力集中が生じるかを予測します。設計段階で試作を繰り返す代わりに、コンピュータ上で安全性を評価できます
- 自動車の衝突シミュレーション: 車体が衝突したとき、エネルギーがどのように吸収され、乗員にどの程度の力が伝わるかを解析します
- 熱解析: 電子機器の放熱設計や、溶接時の温度分布の予測にもFEMが使われています
- 振動・固有値解析: 建築物の地震応答や、機械部品の共振周波数を求めるときにも有限要素法は欠かせません
本記事では、有限要素法がなぜ機能するのかという原理を変分法と最小ポテンシャルエネルギー原理から丁寧に導出し、最も基本的な1次元バー要素から2次元トラス要素までを段階的に解説します。さらに、全ての理論をPythonでスクラッチ実装することで、商用FEMソフトのブラックボックスの中身を自分の手で理解できるようになることを目指します。
本記事の内容
- FEMが必要な理由と基本的な考え方（離散化・近似・組立）
- 変分法と最小ポテンシャルエネルギー原理による定式化
- 1次元バー要素の剛性行列の導出
- 要素剛性行列から全体剛性行列の組立
- 境界条件の適用と連立方程式の求解
- 2次元トラス要素への拡張
- Pythonによるスクラッチ実装とメッシュ収束性の検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 応力とひずみ — 応力 $\sigma$ とひずみ $\varepsilon$ の定義と物理的意味
- フックの法則 — $\sigma = E \varepsilon$ の線形弾性の基礎
- 軸力と伸び — 棒の引張・圧縮における力と変形の関係
- 行列の基礎 — 行列の積・逆行列の基本操作
- 連立方程式 — 連立1次方程式の解法（FEMの最終ステップ）
- 偏微分の基礎 — 多変数関数の微分（変分法の理解に必要）
FEMが必要な理由
解析解の限界
材料力学の教科書で扱う問題は、均一な断面の棒、等分布荷重を受けるはりなど、比較的単純な形状と荷重条件を持つものがほとんどです。これらの問題では微分方程式を解いて厳密解を得ることができます。しかし、現実の構造物はそう簡単ではありません。
たとえば、穴の空いた板に引張荷重がかかる場合を考えてみましょう。穴の周囲では応力集中が生じ、応力分布は非常に複雑になります。このような問題の支配方程式（偏微分方程式）を解析的に解くことは、特殊な形状を除いてほぼ不可能です。
ではどうするか。問題を小さな「要素」に分割し、各要素では単純な近似関数を使い、それらを組み合わせて全体の近似解を得る――これが有限要素法の基本的なアイデアです。
FEMの3つの基本ステップ
有限要素法の手順は、大きく分けて以下の3つのステップで構成されます。
ステップ1: 離散化（メッシュ分割）
連続体（構造物全体）を有限個の小さな要素に分割します。1次元問題では棒を短い区間に分け、2次元問題では三角形や四角形の要素に分割します。分割の結果得られる節点（node）と要素（element）の集合をメッシュと呼びます。
ステップ2: 要素レベルの近似（剛性行列の構築）
各要素の中で変位（変形）を簡単な関数（形状関数）で近似します。この近似を使って、各要素について「力と変位の関係」を行列で表現します。これが要素剛性行列です。バネの「$F = kx$」を行列に拡張したものと思えば直感的です。
ステップ3: 全体系の組立と求解
全ての要素の剛性行列を「重ね合わせ」て、構造全体の剛性方程式を組み立てます。境界条件（固定点や既知変位）を適用した後、得られた連立方程式を解けば、各節点の変位が求まります。変位がわかれば、ひずみ、応力も順に計算できます。
このように、FEMは「分割して、近似して、組み立てる」という明快な手順で動作します。では、この手順がなぜ正しい答えに近づくのか、その数学的な根拠を見ていきましょう。
変分法と最小ポテンシャルエネルギー原理
エネルギー最小の直感
FEMの数学的基礎を理解するために、まず変分法と最小ポテンシャルエネルギー原理について押さえておきましょう。
日常的な直感として、「物体は最もエネルギーの低い状態に落ち着く」という経験を私たちは持っています。ボールを器に置けば底に転がり落ちますし、ゴムバンドを引っ張って手を離せば自然な長さに戻ります。この「エネルギーが最小になる位置に落ち着く」という原理を構造力学に適用したものが、最小ポテンシャルエネルギー原理です。
構造物が外力を受けて変形するとき、内部にはひずみエネルギーが蓄えられます。同時に、外力は仕事をします。構造物の全ポテンシャルエネルギー $\Pi$ は次のように定義されます。
$$ \begin{equation} \Pi = U – W \end{equation} $$
ここで $U$ はひずみエネルギー（内部に蓄えられるエネルギー）、$W$ は外力による仕事です。最小ポテンシャルエネルギー原理は、全ての可能な変形のうち、実際に起こる変形は $\Pi$ を最小にするものであると述べています。
1次元弾性体のひずみエネルギー
具体的に、1次元の弾性棒（断面積 $A$、ヤング率 $E$、長さ $L$）を考えてみましょう。棒の変位を $u(x)$ とすると、ひずみは次のようになります。
$$ \begin{equation} \varepsilon(x) = \frac{du}{dx} \end{equation} $$
フックの法則より応力は $\sigma = E \varepsilon$ ですから、単位体積あたりのひずみエネルギー密度は次式で与えられます。
$$ \begin{equation} u_0 = \frac{1}{2} \sigma \varepsilon = \frac{1}{2} E \varepsilon^2 \end{equation} $$
これを棒全体にわたって積分すると、棒全体のひずみエネルギーが得られます。
$$ \begin{equation} U = \int_0^L \frac{1}{2} E A \left( \frac{du}{dx} \right)^2 dx \end{equation} $$
外力 $f(x)$（単位長さあたりの分布荷重）による仕事は次のようになります。
$$ \begin{equation} W = \int_0^L f(x) \, u(x) \, dx \end{equation} $$
したがって、全ポテンシャルエネルギーは次式です。
$$ \begin{equation} \Pi[u] = \int_0^L \frac{1}{2} E A \left( \frac{du}{dx} \right)^2 dx – \int_0^L f(x) \, u(x) \, dx \end{equation} $$
変分法による平衡方程式の導出
最小ポテンシャルエネルギー原理によると、実際の変位 $u(x)$ は $\Pi[u]$ を最小にします。関数空間上での「最小化」を行うのが変分法です。偏微分の概念を関数空間に拡張したものと捉えると理解しやすくなります。
$\Pi$ を最小にする条件は、$u(x)$ に微小な変分 $\delta u(x)$ を加えたとき、$\Pi$ の変化がゼロになることです。すなわち、$\Pi$ の第一変分が消えるという条件です。
$$ \begin{equation} \delta \Pi = 0 \end{equation} $$
具体的に $\Pi$ の変分を計算してみましょう。$u$ を $u + \delta u$ に置き換えて1次の項を取り出します。
$$ \begin{equation} \delta \Pi = \int_0^L E A \frac{du}{dx} \frac{d(\delta u)}{dx} dx – \int_0^L f(x) \, \delta u(x) \, dx = 0 \end{equation} $$
第1項に部分積分を適用します。$\frac{du}{dx}$ を微分する側、$\delta u$ を積分しない側として部分積分すると、
$$ \begin{equation} \left[ E A \frac{du}{dx} \delta u \right]_0^L – \int_0^L \frac{d}{dx}\left(EA \frac{du}{dx}\right) \delta u \, dx – \int_0^L f(x) \, \delta u \, dx = 0 \end{equation} $$
境界で $\delta u = 0$（固定端の場合）とすると、境界項は消えます。$\delta u$ は任意の関数なので、被積分関数がゼロでなければなりません。こうして次の微分方程式が得られます。
$$ \begin{equation} \frac{d}{dx}\left(EA \frac{du}{dx}\right) + f(x) = 0 \end{equation} $$
これは1次元弾性棒の平衡方程式そのものです。つまり、ポテンシャルエネルギーを最小化するという条件から、力のつり合い式が自然に導かれるのです。カスティリアーノの定理もエネルギー原理の一種ですが、FEMではこの変分法的なアプローチを体系的に使います。
ここまでで、FEMの理論的な基盤となるエネルギー原理を理解しました。次に、この原理を使って具体的にFEMの中核である「剛性行列」を導出していきましょう。
1次元バー要素の剛性行列の導出
バー要素のモデル
FEMの出発点として、最も単純な要素である1次元のバー要素（棒要素）を考えます。バー要素とは、2つの節点を持ち、軸方向の引張・圧縮のみを伝える要素です。日常的には、トラス橋の個々の部材をイメージするとよいでしょう。
1本のバー要素を考えます。左端の節点を節点1（座標 $x_1 = 0$）、右端の節点を節点2（座標 $x_2 = L_e$）とし、要素の長さを $L_e$、断面積を $A_e$、ヤング率を $E_e$ とします。各節点での変位を $u_1$, $u_2$ とします。
形状関数の導入
FEMでは、要素内部の変位を節点変位の「内挿」として近似します。バー要素は2節点なので、変位を1次関数（線形）で近似するのが自然です。
$$ \begin{equation} u(x) = N_1(x) \, u_1 + N_2(x) \, u_2 \end{equation} $$
ここで $N_1(x)$, $N_2(x)$ は形状関数（shape function）と呼ばれ、次の条件を満たします。
- $N_1(0) = 1, \; N_1(L_e) = 0$（節点1では $N_1$ だけが効く）
- $N_2(0) = 0, \; N_2(L_e) = 1$（節点2では $N_2$ だけが効く）
この条件を満たす最もシンプルな1次関数は次の通りです。
$$ \begin{equation} N_1(x) = 1 – \frac{x}{L_e}, \quad N_2(x) = \frac{x}{L_e} \end{equation} $$
行列形式で書くと、次のようにまとめられます。
$$ \begin{equation} u(x) = \begin{bmatrix} N_1 & N_2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \bm{N} \bm{u}_e \end{equation} $$
形状関数のイメージとしては、「節点1に変位1を与えて節点2を固定すると、要素内の変位は $N_1$ の形で直線的に減少する」というものです。$N_2$ も同様で、これら2つの基本パターンの重ね合わせで任意の変位状態を表現します。
ひずみ-変位関係（Bマトリクス）
ひずみは変位の微分ですから、形状関数を微分することで得られます。
$$ \begin{equation} \varepsilon(x) = \frac{du}{dx} = \frac{dN_1}{dx} u_1 + \frac{dN_2}{dx} u_2 \end{equation} $$
形状関数の微分を計算すると、次のようになります。
$$ \begin{equation} \frac{dN_1}{dx} = -\frac{1}{L_e}, \quad \frac{dN_2}{dx} = \frac{1}{L_e} \end{equation} $$
これを行列形式で表現します。
$$ \begin{equation} \varepsilon = \begin{bmatrix} -\frac{1}{L_e} & \frac{1}{L_e} \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \bm{B} \bm{u}_e \end{equation} $$
この $\bm{B}$ をひずみ-変位マトリクス（B-matrix）と呼びます。1次の形状関数を使っているため、$\bm{B}$ は定数（$x$ に依存しない）になっていることに注目してください。これは要素内でひずみが一様であることを意味しています。
要素剛性行列の導出
いよいよ、ポテンシャルエネルギーの最小化から要素剛性行列を導出します。1つの要素のひずみエネルギーは、先ほどの一般式に形状関数の近似を代入すると次のようになります。
$$ \begin{equation} U_e = \int_0^{L_e} \frac{1}{2} E_e A_e \varepsilon^2 \, dx = \frac{1}{2} \bm{u}_e^T \left( \int_0^{L_e} E_e A_e \bm{B}^T \bm{B} \, dx \right) \bm{u}_e \end{equation} $$
$\bm{B}$ が定数なので、積分は簡単に実行できます。
$$ \begin{equation} \int_0^{L_e} E_e A_e \bm{B}^T \bm{B} \, dx = E_e A_e L_e \bm{B}^T \bm{B} \end{equation} $$
具体的に $\bm{B}^T \bm{B}$ を計算しましょう。
$$ \begin{equation} \bm{B}^T \bm{B} = \begin{bmatrix} -1/L_e \\ 1/L_e \end{bmatrix} \begin{bmatrix} -1/L_e & 1/L_e \end{bmatrix} = \frac{1}{L_e^2} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \end{equation} $$
$L_e$ を掛けると、次のようになります。
$$ \begin{equation} \bm{k}_e = \frac{E_e A_e}{L_e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \end{equation} $$
これがバー要素の要素剛性行列です。この行列の物理的な意味を考えてみましょう。剛性行列の各成分は「節点 $j$ に単位変位を与えたとき、節点 $i$ に生じる力」を表しています。
- $k_{11} = EA/L$: 節点1を右に単位変位させると、節点1には右向き（正方向）の反力が $EA/L$ だけ生じる
- $k_{12} = -EA/L$: 節点2を右に単位変位させると、節点1には左向き（負方向）の力が $EA/L$ だけ生じる
軸力と伸びで学んだ $\delta = FL/(EA)$、すなわち $F = (EA/L)\delta$ というばねの関係式が、行列の形で一般化されていることがわかります。
要素の剛性行列を得ることができたので、次のステップとして、複数の要素を「組み立てる」方法を見ていきましょう。
全体剛性行列の組立
組立の考え方
現実の構造物は1つの要素だけではなく、多数の要素で構成されています。全体剛性行列の組立とは、個々の要素の剛性行列を、節点の接続関係に基づいて「足し合わせる」操作です。
この操作の物理的な意味は明快です。2つの要素が1つの節点を共有している場合、その節点にかかる力は両方の要素からの寄与の合計です。数学的には、全ポテンシャルエネルギーが各要素のポテンシャルエネルギーの和になることから導かれます。
$$ \begin{equation} \Pi = \sum_{e=1}^{n_e} \Pi_e = \sum_{e=1}^{n_e} \left( \frac{1}{2} \bm{u}_e^T \bm{k}_e \bm{u}_e – \bm{u}_e^T \bm{f}_e \right) \end{equation} $$
具体例: 3節点2要素の組立
具体的な例で理解しましょう。3つの節点（節点1, 2, 3）と2つのバー要素（要素A: 節点1-2、要素B: 節点2-3）からなる構造を考えます。
要素Aの剛性行列は節点1と節点2に関係し、$2 \times 2$ の行列です。
$$ \begin{equation} \bm{k}_A = \frac{E_A A_A}{L_A} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \end{equation} $$
要素Bの剛性行列は節点2と節点3に関係し、同様に $2 \times 2$ の行列です。
$$ \begin{equation} \bm{k}_B = \frac{E_B A_B}{L_B} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \end{equation} $$
全体系は3つの節点（自由度3）を持つので、全体剛性行列は $3 \times 3$ です。各要素の剛性行列を、全体の節点番号に対応する位置に「配置」して足し合わせます。
要素Aは節点1, 2に対応するので、全体行列の(1,1), (1,2), (2,1), (2,2) 成分に寄与します。要素Bは節点2, 3に対応するので、(2,2), (2,3), (3,2), (3,3) 成分に寄与します。
$$ \begin{equation} \bm{K} = \begin{bmatrix} k_A & -k_A & 0 \\ -k_A & k_A + k_B & -k_B \\ 0 & -k_B & k_B \end{bmatrix} \end{equation} $$
ここで $k_A = E_A A_A / L_A$, $k_B = E_B A_B / L_B$ と略記しました。注目すべき点として、節点2に対応する対角成分 $(2,2)$ が $k_A + k_B$ になっていることが挙げられます。これは、節点2が要素Aと要素Bの両方に接続されているため、両要素の剛性が「並列」に寄与していることを反映しています。
全体剛性方程式
全体剛性行列を使うと、構造全体の力と変位の関係が次の連立方程式で表されます。
$$ \begin{equation} \bm{K} \bm{U} = \bm{F} \end{equation} $$
ここで $\bm{K}$ は全体剛性行列、$\bm{U}$ は全節点の変位ベクトル、$\bm{F}$ は全節点の外力ベクトルです。この式はバネの $F = kx$ を多自由度に一般化したものに他なりません。
ただし、このままでは $\bm{K}$ は特異行列（行列式がゼロ）であり、連立方程式を解くことができません。なぜなら、構造全体が自由に剛体移動できる状態だからです。実際に解を得るには、境界条件を適用して剛体運動を拘束する必要があります。
境界条件の適用
なぜ境界条件が必要か
全体剛性行列 $\bm{K}$ の特異性（逆行列が存在しないこと）は、物理的には構造物が空間中に浮いている状態に対応しています。力を加えると構造全体が加速度的に移動してしまうため、平衡状態が定まりません。現実の構造物は必ずどこかで支持（固定）されているので、その拘束を数学的に表現する必要があります。
適用方法: 行と列の削除法
最もシンプルな方法は、既知変位（通常はゼロ）に対応する行と列を剛性方程式から除去する方法です。
たとえば先ほどの3節点系で、節点1が固定（$u_1 = 0$）されている場合、全体方程式の第1行と第1列を削除して、縮小された方程式を解きます。
$$ \begin{equation} \begin{bmatrix} k_A + k_B & -k_B \\ -k_B & k_B \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} F_2 \\ F_3 \end{bmatrix} \end{equation} $$
この方程式の剛性行列は正則（逆行列が存在する）になっているので、$u_2$ と $u_3$ を求めることができます。
適用方法: ペナルティ法
もう一つの方法として、ペナルティ法があります。これは、固定したい節点の対角成分に非常に大きな値（ペナルティ係数 $\alpha$）を加える方法です。たとえば節点1を固定する場合、$K_{11}$ に $\alpha$（例: $10^{20}$）を加え、同時に $F_1$ にも $\alpha \times 0 = 0$ を加えます（ゼロ変位の場合）。
$$ \begin{equation} \begin{bmatrix} k_A + \alpha & -k_A & 0 \\ -k_A & k_A + k_B & -k_B \\ 0 & -k_B & k_B \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} 0 \\ F_2 \\ F_3 \end{bmatrix} \end{equation} $$
$\alpha$ が十分大きければ、$u_1 \approx 0$ となり、実質的に固定と同じ効果が得られます。ペナルティ法の利点は、行列のサイズを変えずに境界条件を適用できる点にあります。
非ゼロの既知変位 $u_1 = \bar{u}$ を与えたい場合は、$F_1 = \alpha \bar{u}$ とすればよいのです。
ここまでで、1次元バー要素のFEMの全手順を理解しました。Pythonで実際に実装して、理論が正しく機能することを確かめてみましょう。
Python実装1: 1次元バー要素のFEM解析
問題設定
まず、最も基本的な問題として、直列に接続された2本のバー要素（3節点）の問題を解きます。節点1を固定し、節点3に引張力 $P$ を加えます。各要素の断面積やヤング率は異なるものとします。
以下のコードでは、要素剛性行列の生成、全体剛性行列の組立、境界条件の適用、連立方程式の求解、そして応力の計算までをスクラッチで実装します。
import numpy as np
# === 問題パラメータ ===
# 要素1: 節点1-2, 要素2: 節点2-3
E = np.array([200e9, 100e9]) # ヤング率 [Pa] (鋼, アルミ)
A = np.array([1e-4, 2e-4]) # 断面積 [m^2]
L = np.array([1.0, 1.5]) # 要素長さ [m]
n_nodes = 3 # 節点数
n_elements = 2 # 要素数
# 要素の接続情報（各要素の節点番号, 0始まり）
connectivity = np.array([[0, 1],
[1, 2]])
# === ステップ1: 要素剛性行列の生成 ===
def bar_element_stiffness(E_e, A_e, L_e):
"""1次元バー要素の剛性行列を返す"""
k = (E_e * A_e / L_e) * np.array([[ 1, -1],
[-1, 1]])
return k
# === ステップ2: 全体剛性行列の組立 ===
K_global = np.zeros((n_nodes, n_nodes))
for e in range(n_elements):
k_e = bar_element_stiffness(E[e], A[e], L[e])
nodes = connectivity[e] # この要素の節点番号
for i in range(2):
for j in range(2):
K_global[nodes[i], nodes[j]] += k_e[i, j]
print("全体剛性行列 K [N/m]:")
print(K_global)
# === ステップ3: 境界条件と荷重 ===
F = np.zeros(n_nodes)
F[2] = 10000.0 # 節点3に10 kNの引張力
# 境界条件: 節点1を固定 (u_1 = 0)
fixed_dofs = [0]
free_dofs = [i for i in range(n_nodes) if i not in fixed_dofs]
# === ステップ4: 縮小系を解く ===
K_reduced = K_global[np.ix_(free_dofs, free_dofs)]
F_reduced = F[free_dofs]
U_free = np.linalg.solve(K_reduced, F_reduced)
# 全変位ベクトルの復元
U = np.zeros(n_nodes)
U[free_dofs] = U_free
print(f"\n節点変位 [m]: {U}")
# === ステップ5: 各要素のひずみと応力 ===
for e in range(n_elements):
n1, n2 = connectivity[e]
strain_e = (U[n2] - U[n1]) / L[e]
stress_e = E[e] * strain_e
print(f"要素{e+1}: ひずみ = {strain_e:.6e}, 応力 = {stress_e/1e6:.2f} MPa")
# === 検算: 理論解との比較 ===
# 直列バネの伸び: delta = P*L1/(E1*A1) + P*L2/(E2*A2)
delta_1 = F[2] * L[0] / (E[0] * A[0])
delta_2 = F[2] * L[1] / (E[1] * A[1])
print(f"\n--- 理論解との比較 ---")
print(f"理論: u2 = {delta_1:.6e} m, u3 = {delta_1 + delta_2:.6e} m")
print(f"FEM: u2 = {U[1]:.6e} m, u3 = {U[2]:.6e} m")
このコードの実行結果を確認しましょう。全体剛性行列 $\bm{K}$ の $(2,2)$ 成分は、要素1と要素2の剛性 $k_1 = E_1 A_1 / L_1 = 2.0 \times 10^7$ N/m と $k_2 = E_2 A_2 / L_2 \approx 1.333 \times 10^7$ N/m の和になっているはずです。これは節点2が両要素に接続されているという物理的事実を反映しています。
FEM解と理論解が完全に一致することが確認できます。これはバー要素の解が厳密解に一致するという重要な性質です。1次のバー要素は、要素内で断面積やヤング率が一定であり、分布荷重がなければ、実際の変位分布（線形）と一致するため、要素数に関わらず厳密解を与えます。
では次に、もう少し規模の大きな問題で全体剛性行列の組立手順を詳しく見ていきましょう。
Python実装2: 多要素バーの組立と収束性
分布荷重がある場合の問題
先ほどは集中荷重の問題でしたが、今度は分布荷重がかかる1次元問題を考えます。分布荷重がある場合、線形のバー要素だけでは厳密解を捉えきれないため、要素を増やす（メッシュを細かくする）ことで解の精度が向上します。この収束性こそがFEMの本質です。
問題として、長さ $L = 1$ m、断面積 $A = 1 \times 10^{-4}$ m$^2$、ヤング率 $E = 200$ GPaの棒が左端で固定され、自重による一様分布荷重 $q = 1000$ N/mを受ける場合を考えます。
この問題の厳密解は次のとおりです。
$$ \begin{equation} u(x) = \frac{q}{EA}\left(Lx – \frac{x^2}{2}\right) \end{equation} $$
FEMで分布荷重を扱うには、分布荷重を等価節点力に変換する必要があります。要素長 $L_e$ の一様分布荷重 $q$ に対する等価節点力は次のようになります。
$$ \begin{equation} \bm{f}_e = \frac{q L_e}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \end{equation} $$
これは分布荷重を両端の節点に均等に分配するという直感に一致しています。以下のコードでは、要素数を変えて解の収束を確認します。
import numpy as np
import matplotlib.pyplot as plt
def solve_bar_fem(n_elements, L_total, E_val, A_val, q_load):
"""1次元バーFEM解析（一様分布荷重）"""
n_nodes = n_elements + 1
L_e = L_total / n_elements # 各要素の長さ
# 全体剛性行列と荷重ベクトルの組立
K = np.zeros((n_nodes, n_nodes))
F = np.zeros(n_nodes)
for e in range(n_elements):
# 要素剛性行列
ke = (E_val * A_val / L_e) * np.array([[1, -1], [-1, 1]])
# 等価節点力（一様分布荷重）
fe = (q_load * L_e / 2) * np.array([1, 1])
n1, n2 = e, e + 1
K[n1, n1] += ke[0, 0]
K[n1, n2] += ke[0, 1]
K[n2, n1] += ke[1, 0]
K[n2, n2] += ke[1, 1]
F[n1] += fe[0]
F[n2] += fe[1]
# 境界条件: 節点0を固定
free = list(range(1, n_nodes))
K_red = K[np.ix_(free, free)]
F_red = F[free]
U_free = np.linalg.solve(K_red, F_red)
U = np.zeros(n_nodes)
U[free] = U_free
x_nodes = np.linspace(0, L_total, n_nodes)
return x_nodes, U
# パラメータ
L_total = 1.0 # 棒の全長 [m]
E_val = 200e9 # ヤング率 [Pa]
A_val = 1e-4 # 断面積 [m^2]
q_load = 1000.0 # 分布荷重 [N/m]
# 厳密解
x_exact = np.linspace(0, L_total, 200)
u_exact = (q_load / (E_val * A_val)) * (L_total * x_exact - x_exact**2 / 2)
# 要素数を変えてFEM解を計算
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
element_counts = [2, 4, 8, 16]
colors = ['#e74c3c', '#3498db', '#2ecc71', '#9b59b6']
for n_elem, color in zip(element_counts, colors):
x_fem, u_fem = solve_bar_fem(n_elem, L_total, E_val, A_val, q_load)
axes[0].plot(x_fem, u_fem * 1e6, 'o-', color=color,
label=f'{n_elem} elements', markersize=5)
axes[0].plot(x_exact, u_exact * 1e6, 'k--', linewidth=2, label='Exact')
axes[0].set_xlabel('Position x [m]')
axes[0].set_ylabel('Displacement u [μm]')
axes[0].set_title('FEM Solutions vs Exact Solution')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# メッシュ収束性の検証
n_list = [2, 4, 8, 16, 32, 64, 128]
errors = []
for n_elem in n_list:
x_fem, u_fem = solve_bar_fem(n_elem, L_total, E_val, A_val, q_load)
# 先端（x=L）での変位の誤差
u_tip_exact = (q_load / (E_val * A_val)) * (L_total**2 / 2)
error = abs(u_fem[-1] - u_tip_exact) / u_tip_exact
errors.append(error)
axes[1].loglog(n_list, errors, 'bo-', markersize=8, linewidth=2)
axes[1].set_xlabel('Number of elements')
axes[1].set_ylabel('Relative error at tip')
axes[1].set_title('Mesh Convergence')
axes[1].grid(True, alpha=0.3, which='both')
# 収束次数の表示
if len(errors) > 1:
slope = np.log(errors[-1]/errors[-2]) / np.log(n_list[-1]/n_list[-2])
axes[1].text(0.05, 0.95, f'Convergence rate ≈ {abs(slope):.1f}',
transform=axes[1].transAxes, fontsize=12,
verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.tight_layout()
plt.savefig('bar_fem_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
左のグラフから、要素数が増えるにつれてFEM解が厳密解（放物線）に近づいていく様子が見て取れます。要素数2のときはまだ明らかな差がありますが、8要素もあればほぼ厳密解と区別がつかなくなります。
右のグラフは、先端変位の相対誤差を両対数プロットしたものです。直線的に減少していることから、誤差が要素数（あるいはメッシュサイズ）のべき乗で減少することがわかります。1次バー要素の場合、先端変位の収束次数は2次（要素数を2倍にすると誤差が $1/4$ になる）であることが理論的に知られており、グラフからもそれが確認できます。この収束特性は、FEMの信頼性を保証する重要な性質です。
1次元問題でFEMの基本が理解できたところで、いよいよ2次元のトラス問題に拡張していきましょう。
2次元トラス要素への拡張
トラス構造とは
トラス構造とは、棒材（部材）がピン接合で結ばれた構造物のことです。橋梁、送電塔、宇宙ステーションの構造体などがトラスの例です。各部材は軸力（引張または圧縮）のみを伝え、曲げモーメントは伝えません。これは先ほどの1次元バー要素と同じ性質です。
では、1次元バー要素と何が違うのでしょうか。2次元トラスでは、部材が様々な角度で配置されるため、局所座標系（部材方向に沿った座標系）と全体座標系（構造全体に共通の $x$-$y$ 座標系）の間の座標変換が必要になります。
座標変換
長さ $L_e$ のトラス部材が、全体座標系の $x$ 軸と角度 $\theta$ をなしているとします。この部材の節点1の座標を $(x_1, y_1)$、節点2の座標を $(x_2, y_2)$ とすると、方向余弦は次のようになります。
$$ \begin{equation} c = \cos\theta = \frac{x_2 – x_1}{L_e}, \quad s = \sin\theta = \frac{y_2 – y_1}{L_e} \end{equation} $$
ここで要素長さは次式で計算されます。
$$ \begin{equation} L_e = \sqrt{(x_2 – x_1)^2 + (y_2 – y_1)^2} \end{equation} $$
全体座標系での要素剛性行列
局所座標系（部材軸方向）でのバー要素の剛性行列は先ほど導出した $2 \times 2$ 行列ですが、2次元では各節点に $x$ 方向と $y$ 方向の2つの自由度があるため、要素の剛性行列は $4 \times 4$ になります。
座標変換行列 $\bm{T}$ を定義します。局所変位 $\bm{u}^{\text{local}}$ と全体変位 $\bm{u}^{\text{global}}$ の関係は次式です。
$$ \begin{equation} \bm{u}^{\text{local}} = \bm{T} \bm{u}^{\text{global}} \end{equation} $$
変換行列は次の通りです。
$$ \begin{equation} \bm{T} = \begin{bmatrix} c & s & 0 & 0 \\ 0 & 0 & c & s \end{bmatrix} \end{equation} $$
全体座標系での要素剛性行列は、$\bm{k}^{\text{local}}$ を座標変換して得られます。
$$ \begin{equation} \bm{k}_e^{\text{global}} = \bm{T}^T \bm{k}_e^{\text{local}} \bm{T} \end{equation} $$
この計算を実行すると、次の $4 \times 4$ 行列が得られます。
$$ \begin{equation} \bm{k}_e = \frac{E_e A_e}{L_e} \begin{bmatrix} 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{bmatrix} \end{equation} $$
この行列の構造を観察すると、いくつかの重要な性質がわかります。
- 対称性: $\bm{k}_e = \bm{k}_e^T$ — 作用反作用の法則を反映しています
- 特異性: 行列式はゼロ — 単独の要素は剛体移動が可能です
- 角度依存性: $\theta = 0$（水平部材）のとき $s = 0$ となり、$y$ 方向の自由度に関する成分がすべてゼロになります。これは水平バー要素が $y$ 方向の力を伝えないという物理を正しく表現しています
全体系の組立
2次元トラスの全体剛性行列の組立は、1次元のときと同じ「足し合わせ」の手順で行います。ただし、各節点が $x$, $y$ の2つの自由度を持つため、節点番号 $i$ のグローバル自由度番号は $2i$ と $2i+1$ に対応します。
組立の手順をまとめると次の通りです。
- 全体剛性行列 $\bm{K}$（サイズ $2n \times 2n$, $n$ は節点数）をゼロで初期化
- 各要素 $e$ について: – 方向余弦 $c$, $s$ と要素長 $L_e$ を計算 – $4 \times 4$ の全体座標系要素剛性行列 $\bm{k}_e$ を計算 – 要素の節点番号からグローバル自由度番号を取得 – $\bm{k}_e$ の各成分を $\bm{K}$ の対応する位置に加算
- 境界条件を適用して縮小方程式を解く
理論が揃ったので、2次元トラスの完全なFEMソルバーをPythonで実装しましょう。
Python実装3: 2Dトラスの変位・応力解析
問題設定
古典的な2次元トラス問題として、以下の構造を解析します。4つの節点と5本の部材からなる平面トラスで、上部節点に斜め荷重が作用します。
節点2 -------- 節点3
/| |
/ | |
/ | |
節点0 ---------- 節点1
(固定) (ローラー: y方向固定)
import numpy as np
import matplotlib.pyplot as plt
# === 2Dトラス FEMソルバー ===
def truss2d_element_stiffness(E, A, node1, node2):
"""2次元トラス要素の全体座標系での剛性行列（4x4）を返す"""
dx = node2[0] - node1[0]
dy = node2[1] - node1[1]
L = np.sqrt(dx**2 + dy**2)
c = dx / L # cos(theta)
s = dy / L # sin(theta)
k_local = (E * A / L)
k_e = k_local * 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_e, L, c, s
def solve_truss2d(nodes, elements, E_list, A_list, bc_dict, forces):
"""
2Dトラス構造をFEMで解析する
Parameters
----------
nodes : ndarray, shape (n_nodes, 2) — 節点座標
elements : list of (i, j) — 要素の接続情報
E_list : list — 各要素のヤング率
A_list : list — 各要素の断面積
bc_dict : dict — {自由度番号: 既知変位値}
forces : ndarray, shape (2*n_nodes,) — 荷重ベクトル
Returns
-------
U : ndarray — 全節点変位
stresses : list — 各要素の応力
"""
n_nodes = len(nodes)
n_dofs = 2 * n_nodes
K = np.zeros((n_dofs, n_dofs))
# 全体剛性行列の組立
for e, (n1, n2) in enumerate(elements):
k_e, L_e, c_e, s_e = truss2d_element_stiffness(
E_list[e], A_list[e], nodes[n1], nodes[n2])
# グローバル自由度番号
dofs = [2*n1, 2*n1+1, 2*n2, 2*n2+1]
for i in range(4):
for j in range(4):
K[dofs[i], dofs[j]] += k_e[i, j]
# 境界条件の適用（行列縮小法）
all_dofs = set(range(n_dofs))
fixed_dofs = sorted(bc_dict.keys())
free_dofs = sorted(all_dofs - set(fixed_dofs))
K_red = K[np.ix_(free_dofs, free_dofs)]
F_red = forces[free_dofs]
# 既知変位の寄与を右辺に移項
for fd in fixed_dofs:
if bc_dict[fd] != 0.0:
for i, fi in enumerate(free_dofs):
F_red[i] -= K[fi, fd] * bc_dict[fd]
# 連立方程式の求解
U_free = np.linalg.solve(K_red, F_red)
U = np.zeros(n_dofs)
for fd in fixed_dofs:
U[fd] = bc_dict[fd]
for i, fi in enumerate(free_dofs):
U[fi] = U_free[i]
# 各要素の応力を計算
stresses = []
for e, (n1, n2) in enumerate(elements):
dx = nodes[n2][0] - nodes[n1][0]
dy = nodes[n2][1] - nodes[n1][1]
L_e = np.sqrt(dx**2 + dy**2)
c_e = dx / L_e
s_e = dy / L_e
# 局所変位（軸方向成分のみ抽出）
u_local_1 = c_e * U[2*n1] + s_e * U[2*n1+1]
u_local_2 = c_e * U[2*n2] + s_e * U[2*n2+1]
strain_e = (u_local_2 - u_local_1) / L_e
stress_e = E_list[e] * strain_e
stresses.append(stress_e)
return U, stresses, K
# === 問題の定義 ===
# 節点座標 [m]
nodes = np.array([
[0.0, 0.0], # 節点0: 左下（ピン支持）
[4.0, 0.0], # 節点1: 右下（ローラー支持）
[0.0, 3.0], # 節点2: 左上
[4.0, 3.0], # 節点3: 右上
])
# 要素の接続情報 (節点番号のペア)
elements = [
(0, 1), # 要素0: 下弦材（水平）
(0, 2), # 要素1: 左柱（鉛直）
(1, 3), # 要素2: 右柱（鉛直）
(2, 3), # 要素3: 上弦材（水平）
(0, 3), # 要素4: 斜材（対角）
]
# 材料特性（全要素同一）
E_val = 200e9 # ヤング率 200 GPa（鋼）
A_val = 1e-4 # 断面積 1 cm^2
E_list = [E_val] * len(elements)
A_list = [A_val] * len(elements)
# 境界条件: 節点0はピン支持（x,y固定）, 節点1はローラー（y固定）
bc = {0: 0.0, 1: 0.0, # 節点0: u_x=0, u_y=0
3: 0.0} # 節点1: u_y=0
# 荷重: 節点3に下向き20kN, 右向き10kN
forces = np.zeros(2 * len(nodes))
forces[2*3] = 10000.0 # 節点3: F_x = 10 kN
forces[2*3 + 1] = -20000.0 # 節点3: F_y = -20 kN
# === 解析実行 ===
U, stresses, K_global = solve_truss2d(nodes, elements, E_list, A_list, bc, forces)
# 結果表示
print("=== 節点変位 ===")
for i in range(len(nodes)):
print(f" 節点{i}: u_x = {U[2*i]*1e3:.4f} mm, u_y = {U[2*i+1]*1e3:.4f} mm")
print("\n=== 要素応力 ===")
elem_names = ['下弦材(0-1)', '左柱(0-2)', '右柱(1-3)',
'上弦材(2-3)', '斜材(0-3)']
for e, name in enumerate(elem_names):
state = "引張" if stresses[e] > 0 else "圧縮"
print(f" {name}: σ = {stresses[e]/1e6:.2f} MPa ({state})")
この結果から、トラスの変形状態と応力分布の全貌がわかります。各節点の変位（mm単位）と各部材の応力（MPa単位）が出力されます。応力の正負から、各部材が引張状態か圧縮状態かも判別できます。
斜材（対角部材）が最も大きな応力を受けることが多いのは、斜材が荷重を斜めに伝える役割を持ち、軸力成分が集中するためです。また、不静定問題のように冗長な部材がある場合、荷重は剛性に応じて分配されるため、FEMでなければ簡単に解けない問題も体系的に処理できます。
次に、変形形状を可視化して、結果が物理的に妥当かどうかを直感的に確認しましょう。
Python実装4: トラスの変形可視化
変位の数値だけではなく、変形前後の形状を図示することで、解が物理的に正しいかどうかを視覚的に判断できます。FEMの結果をポストプロセスして可視化する能力は、実務でも非常に重要です。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
def plot_truss_deformation(nodes, elements, U, stresses, scale=50,
elem_names=None):
"""トラスの変形前後の形状と応力分布を可視化する"""
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# --- 左図: 変形前後の形状 ---
ax = axes[0]
n_nodes = len(nodes)
# 変形後の節点座標
nodes_def = nodes.copy()
for i in range(n_nodes):
nodes_def[i, 0] += U[2*i] * scale
nodes_def[i, 1] += U[2*i+1] * scale
# 変形前（灰色の破線）
for n1, n2 in elements:
ax.plot([nodes[n1, 0], nodes[n2, 0]],
[nodes[n1, 1], nodes[n2, 1]],
'k--', linewidth=1, alpha=0.4)
# 変形後（応力に応じた色）
stress_arr = np.array(stresses)
max_stress = np.max(np.abs(stress_arr))
for e, (n1, n2) in enumerate(elements):
# 引張=青, 圧縮=赤
if stresses[e] > 0:
color = plt.cm.Blues(0.3 + 0.7 * abs(stresses[e]) / max_stress)
else:
color = plt.cm.Reds(0.3 + 0.7 * abs(stresses[e]) / max_stress)
ax.plot([nodes_def[n1, 0], nodes_def[n2, 0]],
[nodes_def[n1, 1], nodes_def[n2, 1]],
color=color, linewidth=3, solid_capstyle='round')
# 節点番号
for i in range(n_nodes):
ax.plot(nodes[i, 0], nodes[i, 1], 'ko', markersize=6)
ax.annotate(f' {i}', nodes[i], fontsize=11, fontweight='bold')
# 支持条件の記号
ax.annotate('△ Pin', xy=(0, -0.3), fontsize=9, color='gray',
ha='center')
ax.annotate('▽ Roller', xy=(4, -0.3), fontsize=9, color='gray',
ha='center')
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_title(f'Truss Deformation (scale × {scale})')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
# --- 右図: 応力の棒グラフ ---
ax2 = axes[1]
if elem_names is None:
elem_names = [f'Elem {e}' for e in range(len(elements))]
colors_bar = ['#3498db' if s > 0 else '#e74c3c' for s in stresses]
bars = ax2.barh(range(len(stresses)),
[s/1e6 for s in stresses],
color=colors_bar, edgecolor='white', linewidth=0.5)
ax2.set_yticks(range(len(stresses)))
ax2.set_yticklabels(elem_names)
ax2.set_xlabel('Stress [MPa]')
ax2.set_title('Element Stresses (Blue: Tension, Red: Compression)')
ax2.axvline(x=0, color='black', linewidth=0.5)
ax2.grid(True, alpha=0.3, axis='x')
plt.tight_layout()
plt.savefig('truss_deformation.png', dpi=150, bbox_inches='tight')
plt.show()
# 前のセクションで計算した結果を使って可視化
plot_truss_deformation(
nodes, elements, U, stresses, scale=50,
elem_names=['下弦材(0-1)', '左柱(0-2)', '右柱(1-3)',
'上弦材(2-3)', '斜材(0-3)']
)
左図では、変形前の形状（灰色の破線）と変形後の形状（色付き実線）が重ねて表示されます。変形は実際には非常に小さい（mm単位）ため、50倍にスケールアップして描画しています。変形の方向が荷重の作用方向と整合していること、固定点（節点0）では変位がゼロであること、ローラー支持（節点1）では $y$ 方向の変位がゼロで $x$ 方向にはスライドしていることが確認できます。
右図の棒グラフは、各部材の応力を可視化しています。青色が引張、赤色が圧縮を示します。荷重が下向きにかかっているため、鉛直部材は圧縮を受け、斜材は引張力で荷重を支えていることがわかります。これはトラス構造の典型的な力の流れです。
ここまでで、2次元トラスの解析と可視化が完成しました。最後に、メッシュの細かさと解の精度の関係を、連続体の問題で定量的に検証してみましょう。
Python実装5: メッシュ収束性の検証
収束性の重要性
FEMでは要素を細かくすれば解の精度が上がることを直感的に述べましたが、実用上は「どの程度細かくすれば十分か」を定量的に評価する必要があります。これがメッシュ収束性（mesh convergence）の検証です。
メッシュ収束性の検証は、FEMを使う全ての技術者が行うべき基本的な作業です。要素数を増やしていったとき、着目する物理量（変位、応力など）が安定した値に収束しているかどうかを確認します。
収束のメカニズム
1次のバー要素を使った場合、先端変位の誤差は要素サイズ $h$ に対して $O(h^2)$ で減少します。これは要素数 $n$ に対して $O(1/n^2)$ の収束です。より高次の形状関数（2次関数など）を使えば、$O(h^3)$ やそれ以上の高速な収束が得られます。
以下のコードでは、分布荷重を受ける棒の問題で、1次要素と2次要素の両方を実装して収束速度を比較します。2次要素は要素内部に1つの中間節点を追加し、変位を2次関数で近似する要素です。
import numpy as np
import matplotlib.pyplot as plt
def solve_bar_fem_linear(n_elements, L_total, E_val, A_val, q_load):
"""1次バー要素（線形要素）でのFEM解析"""
n_nodes = n_elements + 1
h = L_total / n_elements
K = np.zeros((n_nodes, n_nodes))
F = np.zeros(n_nodes)
for e in range(n_elements):
# 要素剛性行列（1次）
ke = (E_val * A_val / h) * np.array([[1, -1], [-1, 1]])
# 等価節点力
fe = (q_load * h / 2) * np.array([1, 1])
n1, n2 = e, e + 1
K[n1, n1] += ke[0, 0]; K[n1, n2] += ke[0, 1]
K[n2, n1] += ke[1, 0]; K[n2, n2] += ke[1, 1]
F[n1] += fe[0]; F[n2] += fe[1]
# 境界条件
free = list(range(1, n_nodes))
U_free = np.linalg.solve(K[np.ix_(free, free)], F[free])
U = np.zeros(n_nodes)
U[free] = U_free
return U[-1] # 先端変位
def solve_bar_fem_quadratic(n_elements, L_total, E_val, A_val, q_load):
"""2次バー要素（3節点要素）でのFEM解析"""
# 2次要素: 各要素に3節点（端点2つ + 中点1つ）
n_nodes = 2 * n_elements + 1
h = L_total / n_elements # 要素長
K = np.zeros((n_nodes, n_nodes))
F = np.zeros(n_nodes)
for e in range(n_elements):
n1 = 2 * e # 左端節点
n2 = 2 * e + 1 # 中点節点
n3 = 2 * e + 2 # 右端節点
# 2次要素の剛性行列（解析積分）
ke = (E_val * A_val / (3 * h)) * np.array([
[ 7, -8, 1],
[-8, 16, -8],
[ 1, -8, 7]
])
# 等価節点力（2次の形状関数に基づく）
fe = (q_load * h / 6) * np.array([1, 4, 1])
local_nodes = [n1, n2, n3]
for i in range(3):
for j in range(3):
K[local_nodes[i], local_nodes[j]] += ke[i, j]
F[local_nodes[i]] += fe[i]
# 境界条件
free = list(range(1, n_nodes))
U_free = np.linalg.solve(K[np.ix_(free, free)], F[free])
U = np.zeros(n_nodes)
U[free] = U_free
return U[-1] # 先端変位
# パラメータ
L = 1.0
E = 200e9
A = 1e-4
q = 1000.0
# 厳密解（先端変位）
u_exact = q * L**2 / (2 * E * A)
# 収束性の検証
n_list = [1, 2, 4, 8, 16, 32, 64]
errors_linear = []
errors_quadratic = []
for n in n_list:
u_lin = solve_bar_fem_linear(n, L, E, A, q)
errors_linear.append(abs(u_lin - u_exact) / u_exact)
u_quad = solve_bar_fem_quadratic(n, L, E, A, q)
err_q = abs(u_quad - u_exact) / u_exact
errors_quadratic.append(err_q if err_q > 1e-16 else 1e-16)
# 可視化
fig, ax = plt.subplots(figsize=(9, 6))
ax.loglog(n_list, errors_linear, 'bo-', markersize=8, linewidth=2,
label='Linear element (1st order)')
ax.loglog(n_list, errors_quadratic, 'rs-', markersize=8, linewidth=2,
label='Quadratic element (2nd order)')
# 参考線（理論的な収束次数）
n_ref = np.array([2, 64])
ax.loglog(n_ref, 0.5 * (n_ref[0]/n_ref)**2,
'b--', alpha=0.4, label=r'$O(n^{-2})$ reference')
ax.loglog(n_ref, 0.1 * (n_ref[0]/n_ref)**4,
'r--', alpha=0.4, label=r'$O(n^{-4})$ reference')
ax.set_xlabel('Number of elements', fontsize=12)
ax.set_ylabel('Relative error in tip displacement', fontsize=12)
ax.set_title('Mesh Convergence: Linear vs Quadratic Elements', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
ax.set_ylim(1e-16, 1)
plt.tight_layout()
plt.savefig('mesh_convergence_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
# 収束次数の数値確認
print("=== 収束次数の確認 ===")
for i in range(1, len(n_list)):
if errors_linear[i] > 0 and errors_linear[i-1] > 0:
rate_lin = np.log(errors_linear[i-1] / errors_linear[i]) / np.log(n_list[i] / n_list[i-1])
else:
rate_lin = float('inf')
if errors_quadratic[i] > 1e-16 and errors_quadratic[i-1] > 1e-16:
rate_quad = np.log(errors_quadratic[i-1] / errors_quadratic[i]) / np.log(n_list[i] / n_list[i-1])
else:
rate_quad = float('inf')
print(f" n={n_list[i-1]}→{n_list[i]}: "
f"Linear rate={rate_lin:.2f}, Quadratic rate={rate_quad:.2f}")
このグラフから、FEMの収束特性に関する2つの重要な事実が読み取れます。
-
1次要素（線形要素）は $O(n^{-2})$ で収束する: 要素数を2倍にすると、誤差はおよそ $1/4$ に減少します。これは理論予測と完全に一致しています。数値積分における台形則の精度と同じ $O(h^2)$ の収束です
-
2次要素は $O(n^{-4})$ で収束する: 要素数を2倍にすると、誤差は $1/16$ に減少します。たった数個の2次要素で、数十個の1次要素と同等以上の精度が得られることがわかります。この問題では2次要素1つでも厳密解（機械精度）に到達しています
この結果は実務上も重要な意味を持ちます。精度が不足している場合、単純にメッシュを細かくする（h-refinement）か、要素の次数を上げる（p-refinement）かの2つの戦略があり、問題の性質に応じて適切な方を選ぶ必要があります。
FEMの理論的基盤: 重み付き残差法からの見方
ガラーキン法
ここまではエネルギー原理（変分法）からFEMを導出しましたが、もう一つの重要な視点として重み付き残差法、特にガラーキン法からFEMを理解することができます。この視点は、エネルギー原理が直接適用できない問題（例: 熱伝導、流体問題）にもFEMを拡張する際に不可欠です。
考え方はこうです。支配方程式（平衡方程式など）を近似的に満たす解を求めたいとき、方程式の残差（residual）をゼロにすることはできませんが、残差に重み関数を掛けて積分したものをゼロにすることならできます。
棒の平衡方程式を例にとりましょう。
$$ \begin{equation} \frac{d}{dx}\left(EA\frac{du}{dx}\right) + f(x) = 0 \end{equation} $$
近似解 $u_h(x)$ をこの方程式に代入すると、一般には残差 $R(x)$ が生じます。
$$ \begin{equation} R(x) = \frac{d}{dx}\left(EA\frac{du_h}{dx}\right) + f(x) \neq 0 \end{equation} $$
重み付き残差法では、適切な重み関数 $w(x)$ を使って次の積分条件を課します。
$$ \begin{equation} \int_0^L w(x) \, R(x) \, dx = 0 \end{equation} $$
ガラーキン法は、重み関数を近似関数と同じ形状関数から選ぶ方法です。つまり $w(x) = N_i(x)$ とします。これを代入して部分積分を行うと、先ほどの変分法から得られたのと全く同じ剛性行列が得られます。
具体的に、重み関数 $w = N_i$ を使って残差の重み付き積分を計算し、部分積分を適用すると次のようになります。
$$ \begin{equation} \int_0^{L_e} EA \frac{dN_i}{dx} \frac{dN_j}{dx} dx \cdot u_j = \int_0^{L_e} N_i f(x) \, dx \end{equation} $$
左辺の積分が剛性行列の $(i,j)$ 成分に対応し、右辺が等価節点力に対応しています。この結果が変分法から得た式と一致するのは偶然ではなく、自己随伴な問題ではガラーキン法と変分法は数学的に等価であることが証明されています。
2つの視点の使い分け
変分法とガラーキン法は、同じFEMの定式化に至る2つの道筋です。
- 変分法（エネルギー原理）: 物理的な直感が得やすく、解が「エネルギー最小」を満たすという明確な最適化問題として理解できます。構造力学の問題に自然に適用できます
- ガラーキン法（重み付き残差法）: より一般的で、エネルギー原理が存在しない問題（対流・拡散方程式など）にも適用できます。数学的な汎用性が高いのが利点です
どちらの視点から出発しても、最終的に得られる連立方程式 $\bm{K}\bm{U} = \bm{F}$ は同じです。問題の性質と目的に応じて、理解しやすい方を選べばよいでしょう。
この重み付き残差法の視点は、FEMを熱伝導や電磁場の問題に適用する際にも本質的な役割を果たします。
FEMの精度と要素選択の指針
誤差の種類
FEMの解には、大きく分けて3種類の誤差が含まれます。
離散化誤差: 連続体を有限個の要素で近似することによる誤差です。メッシュを細かくすれば減少します。これは先ほどの収束性検証で確認した誤差です。
モデリング誤差: 物理モデルの仮定（線形弾性、微小変形など）が現実と乖離することによる誤差です。これはFEM固有の誤差ではなく、物理モデルの選択の問題です。
数値誤差: 浮動小数点演算の丸め誤差、連立方程式の解法に起因する誤差です。通常は他の誤差に比べて十分小さいですが、条件数の悪い問題では注意が必要です。
アプリオリ誤差評価
FEMの離散化誤差に対しては、次のようなアプリオリ（事前）誤差評価が理論的に得られています。
$$ \begin{equation} \|u – u_h\| \leq C h^p \end{equation} $$
ここで $u$ は厳密解、$u_h$ はFEM解、$h$ は要素サイズ、$p$ は形状関数の次数に依存する定数、$C$ は問題に依存する定数です。この式は、要素サイズ $h$ を半分にすると誤差が $1/2^p$ になることを意味しています。
実務的な注意点
メッシュ設計において注意すべき点をいくつか挙げておきます。
-
応力集中部の局所的な細分化: 応力集中が生じる穴や切欠きの周囲では、メッシュを局所的に細かくする必要があります。全体を一律に細かくするのは計算コストの無駄です
-
要素のアスペクト比: 極端に細長い要素や歪んだ要素は精度を悪化させます。理想的には正方形や正三角形に近い形状の要素を使うべきです
-
収束検証の必要性: どんな問題でも、解が収束しているかを確認することが不可欠です。要素数を増やしても結果が変わらなくなる点を見つけることが、FEM解析の信頼性を担保します
Python実装6: 1次元問題の総合解析
複合棒材の解析
最後に、より実践的な問題として、異なる断面積と材料を持つ複合棒材の問題を解きます。この問題は、段付き軸や異種材料の接合部を模擬しています。棒は一端を固定し、他端に集中荷重を受けます。さらに途中に集中荷重も加えます。
import numpy as np
import matplotlib.pyplot as plt
def solve_composite_bar(sections, loads, fixed_node=0):
"""
複合棒材（段付き軸）のFEM解析
Parameters
----------
sections : list of dict — 各セクションの {E, A, L, n_elem}
loads : dict — {節点番号: 荷重値}
fixed_node : int — 固定節点番号
"""
# 節点とメッシュの構築
total_nodes = sum(s['n_elem'] for s in sections) + 1
K = np.zeros((total_nodes, total_nodes))
F = np.zeros(total_nodes)
# 節点座標と要素情報を構築
x_coords = [0.0]
elem_info = [] # (E, A, L_e, node1, node2)
node_idx = 0
for sec in sections:
L_e = sec['L'] / sec['n_elem']
for _ in range(sec['n_elem']):
x_coords.append(x_coords[-1] + L_e)
elem_info.append({
'E': sec['E'], 'A': sec['A'], 'L_e': L_e,
'n1': node_idx, 'n2': node_idx + 1
})
node_idx += 1
x_coords = np.array(x_coords)
# 全体剛性行列の組立
for info in elem_info:
ke_val = info['E'] * info['A'] / info['L_e']
ke = ke_val * np.array([[1, -1], [-1, 1]])
n1, n2 = info['n1'], info['n2']
K[n1, n1] += ke[0, 0]; K[n1, n2] += ke[0, 1]
K[n2, n1] += ke[1, 0]; K[n2, n2] += ke[1, 1]
# 荷重の適用
for node, load in loads.items():
F[node] += load
# 境界条件
free = [i for i in range(total_nodes) if i != fixed_node]
U_free = np.linalg.solve(K[np.ix_(free, free)], F[free])
U = np.zeros(total_nodes)
U[free] = U_free
# 要素ごとの応力
stresses = []
x_stress = []
for info in elem_info:
n1, n2 = info['n1'], info['n2']
strain = (U[n2] - U[n1]) / info['L_e']
stress = info['E'] * strain
stresses.append(stress)
x_stress.append((x_coords[n1] + x_coords[n2]) / 2)
return x_coords, U, np.array(x_stress), np.array(stresses), elem_info
# === 問題定義: 3セクションの段付き軸 ===
sections = [
{'E': 200e9, 'A': 2e-4, 'L': 0.5, 'n_elem': 4}, # 鋼, 太い
{'E': 200e9, 'A': 1e-4, 'L': 0.3, 'n_elem': 3}, # 鋼, 細い
{'E': 70e9, 'A': 1.5e-4, 'L': 0.4, 'n_elem': 4}, # アルミ
]
# 節点4（セクション1と2の境界）に5kN, 先端に10kN
n_sec1 = sections[0]['n_elem']
n_total = sum(s['n_elem'] for s in sections)
loads = {
n_sec1: 5000.0, # 中間荷重 5 kN
n_total: 10000.0 # 先端荷重 10 kN
}
x, U, x_s, stresses, elem_info = solve_composite_bar(sections, loads)
# === 可視化 ===
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True)
# 断面積の分布
A_plot = [info['A'] * 1e4 for info in elem_info] # cm^2
x_elem = [(x[info['n1']] + x[info['n2']]) / 2 for info in elem_info]
axes[0].bar(x_elem, A_plot, width=0.03, color='steelblue', alpha=0.7,
edgecolor='white')
axes[0].set_ylabel('Cross-section area [cm²]')
axes[0].set_title('Composite Bar Analysis: Cross Section, Displacement, Stress')
axes[0].grid(True, alpha=0.3)
# セクション境界を表示
boundaries = [0]
for sec in sections:
boundaries.append(boundaries[-1] + sec['L'])
for b in boundaries[1:-1]:
for ax in axes:
ax.axvline(x=b, color='gray', linestyle=':', alpha=0.5)
# 変位分布
axes[1].plot(x, U * 1e6, 'b-o', markersize=4, linewidth=2)
axes[1].set_ylabel('Displacement [μm]')
axes[1].grid(True, alpha=0.3)
# 荷重の矢印
for node, load in loads.items():
axes[1].annotate(f'{load/1e3:.0f} kN',
xy=(x[node], U[node]*1e6),
xytext=(x[node], U[node]*1e6 + 20),
arrowprops=dict(arrowstyle='->', color='red'),
fontsize=9, color='red', ha='center')
# 応力分布
colors_stress = ['#2ecc71' if s > 0 else '#e74c3c' for s in stresses]
axes[2].bar(x_s, stresses / 1e6, width=0.025,
color=colors_stress, edgecolor='white', alpha=0.8)
axes[2].set_xlabel('Position x [m]')
axes[2].set_ylabel('Stress [MPa]')
axes[2].axhline(y=0, color='black', linewidth=0.5)
axes[2].grid(True, alpha=0.3)
# セクション注釈
sec_labels = ['Steel\n(A=2cm²)', 'Steel\n(A=1cm²)', 'Aluminum\n(A=1.5cm²)']
for i, sec in enumerate(sections):
x_mid = boundaries[i] + sec['L'] / 2
axes[0].text(x_mid, max(A_plot) * 0.9, sec_labels[i],
ha='center', fontsize=9, style='italic')
plt.tight_layout()
plt.savefig('composite_bar_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
# 結果の数値出力
print("=== 複合棒材解析結果 ===")
print(f"先端変位: {U[-1]*1e6:.2f} μm")
print(f"最大応力: {max(abs(stresses))/1e6:.2f} MPa "
f"(位置 x = {x_s[np.argmax(abs(stresses))]:.3f} m)")
print(f"最小応力: {min(abs(stresses))/1e6:.2f} MPa")
# セクション別の応力表示
print("\n--- セクション別応力 ---")
idx = 0
for i, sec in enumerate(sections):
sec_stresses = stresses[idx:idx+sec['n_elem']]
print(f" セクション{i+1} ({sec_labels[i].replace(chr(10), ', ')}): "
f"σ = {sec_stresses[0]/1e6:.2f} MPa")
idx += sec['n_elem']
この結果から、段付き軸の力学的特性が明確にわかります。
上段のグラフは各セクションの断面積を示しています。鋼の太い部分（セクション1: $A = 2$ cm$^2$）から鋼の細い部分（セクション2: $A = 1$ cm$^2$）を経て、アルミニウム部分（セクション3: $A = 1.5$ cm$^2$）へと材質・断面が変化しています。
中段の変位分布は、セクションごとに傾きが異なる区分的線形関数になっています。各セクションでの傾き（ひずみ）は $\varepsilon = \sigma / E = F / (EA)$ であり、$EA$ の小さいセクションほど急な傾きを持ちます。中間荷重が加わる節点では傾きが不連続に変化していることも確認できます。
下段の応力分布から、断面積が変わるところで応力が不連続にジャンプしていることが読み取れます。これは、軸力が同じでも断面積が異なれば $\sigma = F/A$ が変わるためです。セクション2（鋼、$A = 1$ cm$^2$）で応力が最大になっているのは、この部分が最も細く、かつ全荷重（先端10 kN + 中間5 kN = 15 kN）を受けているからです。設計上、この部分が強度のボトルネックになります。
FEMの拡張と発展
ここで扱った要素の位置づけ
本記事で扱った1次元バー要素と2次元トラス要素は、FEMの最も基礎的な要素です。実際のFEMソフトウェアでは、これらを発展させた多様な要素が用意されています。
はり要素: バー要素に曲げの自由度を追加した要素です。各節点に変位に加えて回転角の自由度を持ちます。はりのたわみの問題を解くことができます。
平面応力・平面ひずみ要素: 2次元の連続体を三角形や四角形の要素で離散化する要素です。板に穴が空いた構造の応力集中なども解析できます。
シェル要素: 薄い殻構造（自動車のボディ、航空機の胴体など）を解析するための要素で、面内力と曲げを同時に扱います。
ソリッド要素: 3次元の任意形状を四面体や六面体で離散化する要素です。最も汎用的ですが、計算コストも大きくなります。
非線形問題への拡張
本記事では線形弾性問題（フックの法則が成り立つ範囲）のみを扱いましたが、実際の問題では非線形性が重要になることがあります。
材料非線形: 塑性変形（永久変形）、クリープ、超弾性などの非線形応力-ひずみ関係を扱う問題です。反復計算（ニュートン-ラフソン法など）が必要になります。
幾何学的非線形: 大変形により変形後の形状と初期形状が大きく異なる場合の問題です。座屈解析もこのカテゴリに含まれます。
接触問題: 複数の物体が接触・離脱を繰り返す問題です。自動車の衝突シミュレーションなどが代表例です。
動的問題
時間に依存する問題（振動、波動伝播、衝撃応答など）もFEMで扱うことができます。空間の離散化にFEMを使い、時間方向には差分法（ニューマーク法、中心差分法など）を組み合わせます。
まとめ
本記事では、有限要素法（FEM）の基礎を、原理の導出からPython実装まで一貫して解説しました。
- FEMの基本思想: 構造物を有限個の要素に分割し、各要素で変位を形状関数で近似し、全体を組み立てて連立方程式を解くという3ステップの手法です
- エネルギー原理: 最小ポテンシャルエネルギー原理から、実際の変形がエネルギーを最小にすることを利用して、剛性行列を体系的に導出できます
- 1次元バー要素: $\bm{k}_e = (EA/L)\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$ という $2 \times 2$ の要素剛性行列が全ての出発点です
- 全体剛性行列の組立: 節点の接続関係に基づいて要素剛性行列を「足し合わせる」ことで全体系を構築します
- 境界条件: 行列縮小法やペナルティ法で拘束を適用し、特異行列を正則にして解を得ます
- 2次元トラス: 座標変換行列 $\bm{T}$ を導入することで、任意の角度の部材を統一的に扱えます
- メッシュ収束性: 要素を細かくするほど解は厳密解に収束し、その収束速度は形状関数の次数で決まります
FEMは「分割して、近似して、組み立てる」という単純な考え方に基づいていますが、その適用範囲は構造力学にとどまらず、熱伝導、電磁場、流体力学など、偏微分方程式で記述されるあらゆる物理現象に拡張できます。本記事で身につけた基礎は、これら全ての応用の出発点となります。
次のステップとして、以下の記事も参考にしてください。
- はりのたわみ — はり要素の導出と曲げ問題への拡張
- 応力集中 — FEMが威力を発揮する典型的な問題
- カスティリアーノの定理 — エネルギー法のもう一つの重要な定理
- 不静定問題 — FEMの強みが活きる冗長構造の解析
- 数値積分 — 高次要素で必要となるガウス求積法の基礎