目的
有限要素法を学ぶにあたって、
森北出版株式会社 竹内則雄他著 「計算力学第２版 有限要素法の基礎」
を読み進めながら、自分なりの理解ポイントをまとめていく。
参考書籍
「計算力学第２版 有限要素法の基礎」森北出版株式会社
日本計算工学会編
竹内則雄、樫山和男、寺田賢二郎共著
第1章 計算力学概論
1.1 計算力学とモデル化
固体力学や流体力学電磁気学などの分野で、現象理解のため、計算力学という、対象を離散化してコンピュータでシミレーションする分野がある。
- 差分法 FDM
- 有限要素法 FEM
- 有限体積法 FVM
- 境界要素法 BEM
CAE=Computer-Aided Engineering
対象物を支配する方程式が存在して、それぞれ相互連関する。
- 固体：平衡方程式、構成方程式
- 流体：連続の式、運動方程式
- 電磁場：Maxwell方程式
- 物質移動：輸送方程式
現象をシミュレートするのに、複数の方程式を組み合わせて、条件づけて解くことになる。
連続の式：
$$ \frac{D \rho}{D t} + \rho \nabla \cdot v = 0 $$
運動方程式
$$ \rho \left( \frac{\partial v}{\partial t}
+ v \cdot \nabla v \right) = \nabla \cdot \sigma + \rho f$$
構成方程式
$$ \sigma = - \nabla p - \mu \nabla^2 v
- (\lambda + \mu) \nabla (\nabla \cdot v) $$
1.2 計算力学の手法
対称方程式を離散化してコンピュータで解く。
空間と時間の両方を離散化するが、空間の離散化には、代表的に差分法、有限要素法、有限体積法、境界要素法、またメッシュフリー法などがある。時間の離散化では、陽解法や陰解法がある。
差分法は、微分を差分によって近似して、離散化した時空間の格子上の値の差分で近似して解く。高階微分は、差分の定義を繰り返して定義する。時空間を正方格子に分割することが多い。時間と空間の刻みを変えることも可能。
\frac{du}{dx} \sim \frac{u_{i+1} - u_i}{h} \\
\frac{d^2 u}{dx^2} \sim \frac{u_{i+2} - 2 u_{i+1} + u_i}{h^2}
有限要素法は、領域を多数の小領域に分割して領域内を単純な関数の重ね合わせとして表現する。有限要素法は、微分方程式を直接解く差分法と異なり、積分形式に変換して解く。領域を要素、節点をノードと呼ぶ。三角形分割することが多いが、四角形でも可能。
有限体積法は、任意形状のセルに分割してセル節点を中心とするコントロールボリューム内で積分して得られる積分形で扱う。有限要素法等は、ドロネー分割とボロノイ図のような双対関係な感じか。
境界要素法は、境界上を分割してその上の未知量を求める。
第2章 数学的基礎
2.1 行列の代数
省略
2.2 場とその演算
省略
ガウスの定理は必要
\int_\Omega \nabla f dV = \int_\Gamma f d {\bf{S}}
$V$はその次元の体積、$\bf{S}$は境界での法線。
第3章
3.1 支配方程式の種類
支配方程式の多くは２階の微分方程式で書かれ、係数の関係による分類される。
A \frac{\partial^2 \phi}{\partial x^2}
+ B \frac{\partial^2 \phi}{\partial x \partial t}
+ C \frac{\partial^2 \phi}{\partial t^2}
+ D \frac{\partial \phi}{\partial x}
+ E \frac{\partial \phi}{\partial t}
+ F \phi + G = 0
- $B^2 - 4 A C < 0 $: 楕円型
- $B^2 - 4 A C = 0 $: 放物型
- $B^2 - 4 A C > 0 $: 双曲型
型の違いにより適切な離散化，解法がある。
定常・非定常かどうか。
3.2 代表的な初期値・境界値問題
楕円型問題
楕円型の代表は、Poisson方程式。
\frac{\partial^2 \phi}{\partial x^2}
+ \frac{\partial^2 \phi}{\partial y^2}
= f(x,y)
定常熱伝導問題などが代表的。
２階微分を含むので、ある点の値は周囲の点の値の平均値となる。
ポテンシャル流れとは、完全流体の渦無し流れ。
以下のポテンシャル$\phi$が存在。
$$\frac{\partial \phi}{\partial x} = -u,
\frac{\partial \phi}{\partial y} = -v $$
$u, v$は速度、$\phi$は速度ポテンシャルと呼ぶ。
質量保存の法則は密度を用いて、オイラーの連続方程式を得る。
$$ \frac{\partial \rho}{\partial t}
+ \frac{\partial (\rho u)}{\partial x}
+ \frac{\partial (\rho v)}{\partial y} = 0 $$
境界条件
- Dirichlet境界条件：境界上で物理量の値を指定
- Neumann境界条件：物理量の微分係数を指定
弾性問題
微小領域の応力釣り合いを考えると、釣り合い方程式・平衡方程式を得る。
\frac{\partial \sigma_x}{\partial x}
+ \frac{\partial \tau_{xy}}{\partial y} + b_x = 0 \\
\frac{\partial \tau_{xy}}{\partial x}
+ \frac{\partial \sigma_y}{\partial y} + b_y = 0
応力とひずみ
ひずみは、変位 $(u, v)$ に対して、
\epsilon_x = \frac{\partial u}{\partial x}, \epsilon_y = \frac{\partial v}{\partial y}
であらわされ、せん断ひずみは、
\gamma_{xy} = \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y}
となる。
平面応力状態：薄い板の面に荷重がかかり面の直交方向に伸長される状態
平面ひずみ状態：長い円筒が軸に垂直な力がかかって断面がゆがむような状態
放物型問題
熱伝導方程式が１例：
\frac{\partial \phi}{\partial t} = \frac{\partial^2 \phi}{\partial x^2}
$\phi$ は温度や濃度など、時間拡散するようなものを扱う。
方程式は、時間変化量は拡散による収支であると言っている。
定常解は楕円型になるので、急変動している部分をおさえる、平均化するように働く。
時間方向は、指数関数的な解となるので、一瞬で無限遠まで影響が届く。
解くときは、方向性がないので、空間方向はガラーキン法、時間方向は陰的方法がよいらしい。
熱伝導問題
熱流束 $q$ は、温度勾配に比例して、
$$ q = - k_x \frac{\partial T}{\partial x} $$
であらわす。 $k_x$ は熱伝導率。
微小領域での熱流束と、単位時間あたりの温度上昇にあたる熱量を等価とすれば、
\rho c \frac{\partial T}{\partial t}
= \frac{\partial }{\partial x} \left(k_x \frac{\partial T}{\partial x} \right)
+ \frac{\partial }{\partial y} \left(k_y \frac{\partial T}{\partial y} \right) + Q
$c$ は比熱、 $\rho$ は密度とする。
双曲型問題
典型例は、波動方程式
\frac{\partial^2 u}{\partial t^2}
= c^2 \frac{\partial^2 u}{\partial x^2}
$c$ は波速。
\left( \frac{\partial}{\partial t}
- c \frac{\partial }{\partial x} \right)
\left( \frac{\partial}{\partial t}
+ c \frac{\partial}{\partial x} \right) u = 0
と変形すれば、
\left( \frac{\partial}{\partial t}
- c \frac{\partial }{\partial x} \right)u = 0 \\
\left( \frac{\partial}{\partial t}
+ c \frac{\partial}{\partial x} \right) u = 0
の解のどちらも、元の方程式を満たす。一般解は、
$$u = u(x \pm ct)$$
の形とできるので、初期形状を保ったまま、位置だけが正または負の方向に移動していく解となる。
方向性があるので解くときは、風上化を施した有限要素法を使い、時間方向は陽的方法がよいらしい。
双曲型と放物型を組み合わせた、移流・拡散方程式がある。
$$
\frac{\partial u}{\partial t}
+ c \frac{\partial u}{\partial x}
- k \frac{\partial^2 u}{\partial x^2}
$$
第４章 マトリックス構造解析
有限要素法の基礎
それぞれの自由度をベクトル、行列化して関係式を求めることで、行列の演算規則を駆使して、システマティックに解を求めることができる。
4.1 マトリックス変位法
変位法による剛性行列、剛性方程式の誘導
各接点ごとの自由度を考えて、①変形条件、②釣り合い条件、③適合条件を構成して、行列の満たす式を作る。変形条件は各ばねの変位と力の関係を、釣り合い条件は外力と力、内部の力のつりあいの条件を、適合条件は、構造的に一体になって動くようなときの変位の整合性の条件を与える。
一方で、部材ごとの視点で条件を考えても、同様の方程式を得ることができる。
4.1 ２次元トラス構造解析への応用
全体座標系と部材座標系
各部材ごとの座標で計算されたものを座標変換をして、全体座標系に統一して計算を行う。
第５章 重み付き残差法と有限要素法
ここが本題。
近似解と誤差
近似解を $U(x)$ として、一時独立な $N$ 個の関数の線形結合を考える。
$$
U(x) = \sum^N_{k=1} \Psi_k (x) c_k
$$
仮に、$\Psi_k = x^{k-1}$ と取って、ディリクレ条件から、係数 $c_k$ の関係性が求められる。係数を整理すると、
U(x) = g(x) + \sum_{k=1}^M \Phi_k(x) \alpha_k \\
g(x) = \left( 1 - \frac{x}{l} \right) u_0 + \frac{x}{l} u_l \\
\Phi_k(x) = x (x^k - l^k)
と、問いのディリクレ境界条件に起因する項と、新たな基底 $Psi_k(x)$ が得られる。一方、 $Phi_k(x)$ は、
$$
\Phi_k(0) = \Phi_k(l) = 0
$$
のように、境界でゼロとなるように整理されている。
近似解に対して、境界条件部分を一部に押し付けて、後半は形状の生成を担っていると考えることができそう。この性質は、重み付き残差法などで基底関数が満たすべき条件も満たしているとのこと。
次に、以下のような関数を考える。
$$ r(U(x)) = \frac{d^2 U(x)}{d x^2} - m^2 U(x) $$
厳密解ならば、方程式の両辺がイコールとなるが、今、近似解を考えるので、方程式の両辺の差は一般にはゼロにはならない。これから、この差分、残差が極力ゼロになるように、係数 $\alpha_k$ を決めていくことになる。
最小二乗法
全ての $x$ においてゼロとできればよいが、近似解のため、条件を緩めて、領域全体の残差の合計を小さくすることを考える。ここでは、2乗したものの領域全体に対する和を考える。
$$
R(U) = \frac{1}{2} \int^l_0 (r(U))^2 dx
$$
まだ係数 $\alpha_k$が決まっていない。この未知パラメータに対して、2乗誤差を最小となるものを考えようとすると、$ \alpha_k$ の停留値を考えればよく、
$$
\frac{\partial R(U)}{\partial \alpha_k} = \int^l_0 \frac{\partial r(U)}{\partial \alpha_k} r(U) dx = 0
$$
となる。たとえば、$k=1, \cdots, n$ とすると、$\alpha_k$ に関する式が $n$ 個得られるので、$\alpha_k$ について解けば、任意の個数の近似項が求められる。
重み付き残差法
2乗誤差による計算は、局所的な重要度を無視して、全体の和が最小となるような値を求めるため、平均的な解となってしまう。求めたい厳密解が緩やかで、値の変動、値の大小差が小さいようなものであればよいが、もし局所的に値が大きくなる、領域全体で大きな値の変化があるようなときだと、大きな値の影響度が大きくなる。データ平均値を求めるときや、線形回帰の線を求めるときに、外れ値の影響が大きく利いてしまうようなイメージ。全員平等に、値を小さくしましょう！な感じ。
そこで、一般的な関数$\xi_i(x)$を考えて、
$$
\int^l_0 \chi_i r(U) dx =
\int^l_0 \chi_i \left[ \sum_{k=1}^M \left(
\frac{d^2 \Phi_k}{dx^2} - m^2 \Phi_k \right) \alpha_k - m^2 g \right] dx = 0
$$
を満たすように、$\alpha_k$ を決めることを考える。