科学

有限要素法の手計算による求解方法は?基本例題で学ぶ計算手順!(剛性マトリックス:節点変位:要素応力:連立方程式など)

当サイトでは記事内に広告を含みます
いつも記事を読んでいただきありがとうございます!!! これからもお役に立てる各情報を発信していきますので、今後ともよろしくお願いします(^^)/

有限要素法を手計算で解くことは、FEMの根本的なアルゴリズムを深く理解するための最も効果的な学習方法です。

ソフトウェアを使った解析は便利ですが、裏側で何が行われているかを理解しないと結果の正しい解釈ができません。

この記事では、最も基本的なトラス・バー要素を用いた1次元問題の手計算例題を通じて、要素剛性マトリックスの組み立て・全体剛性マトリックスの構築・境界条件の適用・連立方程式の求解・応力の計算という一連の手計算プロセスを詳しく解説していきます。

有限要素法の手計算プロセスの基本:バー要素を使った例題

それではまず、最も基本的な1次元バー要素(トラス要素)を使った単純な問題で手計算の流れを解説していきます。

バー要素は軸力のみを伝える最もシンプルな有限要素で、FEMの基本アルゴリズム(要素剛性マトリックスの組み立て・アセンブリ・境界条件適用・求解・応力計算)を全て含むため手計算学習に最適な要素です。

例題の設定

以下の単純な二要素トラス問題を手計算で解いていきます。

例題:長さL₁=0.5m・断面積A₁=1×10⁻³m²・ヤング率E₁=200GPaのバー要素1(節点1〜節点2)と、長さL₂=0.5m・断面積A₂=0.5×10⁻³m²・ヤング率E₂=200GPaのバー要素2(節点2〜節点3)が直列につながった一次元問題。境界条件:節点1は完全固定(u₁=0)。節点3に軸方向荷重F₃=10kNを適用。節点2・3の変位u₂・u₃を求めよ。

Step1:要素剛性マトリックスの計算

各バー要素の要素剛性マトリックス[k]は以下の式で計算されます。

バー要素の剛性マトリックス:[k] = (EA/L) × [ 1 -1 / -1 1 ]。要素1:k₁ = E₁A₁/L₁ = (200×10⁹ × 1×10⁻³) / 0.5 = 4×10⁸ N/m = 400 MN/m。要素2:k₂ = E₂A₂/L₂ = (200×10⁹ × 0.5×10⁻³) / 0.5 = 2×10⁸ N/m = 200 MN/m。

要素1の剛性マトリックス(節点1・2に対応、単位:MN/m):[k₁] = 400×

要素2の剛性マトリックス(節点2・3に対応):[k₂] = 200×

Step2:全体剛性マトリックスのアセンブリ

全体剛性マトリックス[K](3×3)は各要素剛性マトリックスを対応する節点位置に組み込んで構築します。

全体剛性マトリックス[K](単位:MN/m):K₁₁=400、K₁₂=-400、K₁₃=0 / K₂₁=-400、K₂₂=400+200=600、K₂₃=-200 / K₃₁=0、K₃₂=-200、K₃₃=200。全体方程式:[400 -400 0 / -400 600 -200 / 0 -200 200]×{u₁,u₂,u₃}ᵀ = {F₁,F₂,F₃}ᵀ

アセンブリでは共有節点(節点2)の剛性が要素1と要素2の両方から加算(400+200=600)されることがポイントです。

全体剛性マトリックスは帯状スパース行列(非ゼロ成分が対角帯に集中)となり、この構造が大規模問題での効率的な求解を可能にします。

境界条件の適用と連立方程式の求解

続いては、境界条件を全体方程式に適用して連立方程式を解く手順を確認していきます。

Step3:境界条件の適用

境界条件として節点1が完全固定(u₁=0)、節点3に荷重F₃=10kN=0.01MNを適用します。

節点2は自由節点で外力なし(F₂=0)、節点1の反力はF₁(未知の反力)です。

固定境界条件の処理として、u₁=0に対応する第1行・第1列を消去(またはペナルティ法で処理)して2×2の縮約方程式を得ます。

縮約後の方程式(u₁=0を代入):[600 -200 / -200 200]×{u₂,u₃}ᵀ = {0, 0.01}ᵀ(単位:MN/m × m = MN)

Step4:連立方程式の求解

2×2の連立方程式を解きます。

600u₂ – 200u₃ = 0 ···(1) / -200u₂ + 200u₃ = 0.01 ···(2)。(1)式より:u₂ = 200u₃/600 = u₃/3 ···(3)。(3)を(2)に代入:-200(u₃/3) + 200u₃ = 0.01 → u₃(-200/3 + 200) = 0.01 → u₃ × (400/3) = 0.01 → u₃ = 0.01×3/400 = 7.5×10⁻⁵ m = 0.075mm。u₂ = u₃/3 = 2.5×10⁻⁵ m = 0.025mm。

節点2の変位u₂=0.025mm・節点3の変位u₃=0.075mmが得られました。

Step5:要素応力・ひずみの計算

節点変位から各要素の応力・ひずみを計算します。

要素1の応力σ₁:ε₁=(u₂-u₁)/L₁=(0.025×10⁻³ – 0)/0.5 = 5×10⁻⁵。σ₁ = E₁×ε₁ = 200×10⁹ × 5×10⁻⁵ = 10×10⁶ Pa = 10 MPa(引張)。要素2の応力σ₂:ε₂=(u₃-u₂)/L₂=(0.075-0.025)×10⁻³/0.5 = 10⁻⁴。σ₂ = E₂×ε₂ = 200×10⁹ × 10⁻⁴ = 20 MPa(引張)。

軸力に換算すると:要素1の軸力N₁ = σ₁×A₁ = 10×10⁶ × 1×10⁻³ = 10 kN、要素2の軸力N₂ = σ₂×A₂ = 20×10⁶ × 0.5×10⁻³ = 10 kN。

両要素の軸力が10kNで一致しており、外力との釣り合いが満足されていることが確認でき、計算の正しさが検証されます。

手計算で学ぶFEMの重要な学習ポイント

続いては、手計算を通じて理解が深まるFEMの重要な概念と学習ポイントを確認していきます。

剛性マトリックスの物理的意味

剛性マトリックスの各成分Kᵢⱼは「節点jに単位変位を与えたときに節点iに生じる力」という物理的意味を持ちます。

対角成分(Kᵢᵢ)は常に正で節点iの自己剛性、非対角成分(Kᵢⱼ:i≠j)は節点間の結合剛性を表します。

全体剛性マトリックスは対称(Kᵢⱼ=Kⱼᵢ:マクスウェルの相反定理)・正定値(全節点を固定しない限り)という数学的特性を持ちます。

境界条件処理の方法の比較

固定変位境界条件の処理方法には消去法・ペナルティ法・ラグランジュ乗数法があります。

消去法は対応する行・列を削除して次元を縮小する方法で、本手計算例で採用した最も直接的な方法です。

ペナルティ法は剛性に大きな値(ペナルティ係数)を対角成分に加算して境界条件を強制する方法で、プログラム実装が容易です。

ラグランジュ乗数法は境界条件を厳密に満足させる最も数学的に厳密な方法ですが、方程式の次元が増加します。

手計算からFEMソフトウェアへの橋渡し

手計算でFEMの本質(剛性マトリックス・アセンブリ・境界条件・求解・後処理)を体感することで、FEMソフトウェアが内部で何をしているかを深く理解できます。

この理解があることで、ソフトウェアのブラックボックス的な操作から脱却し、解析結果の妥当性を自ら判断できる力が身につきます。

PythonやMATLABを使って本手計算例をプログラム化することも、FEMアルゴリズムの理解を深める効果的な学習方法として推奨されます。

手計算によるFEM学習は解析エンジニアとして最も重要な「結果を正しく判断する力」の育成において不可欠なプロセスといえるでしょう。

まとめ

この記事では、有限要素法の手計算による求解方法として二要素バー問題を例題に、要素剛性マトリックスの計算・全体剛性マトリックスのアセンブリ・境界条件の適用・連立方程式の求解・応力計算という一連の手計算プロセスを解説しました。

手計算を通じて剛性マトリックスの物理的意味・アセンブリの仕組み・境界条件処理の方法という有限要素法の本質的な概念を体感的に理解することができます。

FEMソフトウェアの結果を正しく解釈し信頼性を評価するためには、手計算レベルでのアルゴリズム理解が解析エンジニアとしての基礎力として欠かせません。

本例題を出発点として、梁要素・2次元平面要素・熱要素へと学習を拡張していくことで、実用的な有限要素解析能力を体系的に習得できるでしょう。