前回は、有限要素法が収束性するための条件と、誤差の発生要因を解説しました。今回は、CAE での数値解法としてよく用いられる有限要素法と、他数値解法との比較として、有限差分法と境界要素法を取り上げ、それぞれの特徴を解説します。
第4回「有限要素法と他数値解法との比較」を解説していきます。
第4回もくじ
1. 有限差分法
米国の工学者S. P. ティモシェンコは、2 次元弾性問題への有限差分方程式の適用は、1908 年のドイツの数学者C. ルンゲによるねじり問題が最初だと述べています（参考：S.P.Timoshenko、J.N.Goodier：Theoryof Elasticity、Appendix、McGraw-Hill、1951 年）。差分法は、偏微分方程式の最も古い数値解法です。
・有限差分法とは
有限差分法とは、熱伝導や軸のねじれなど、対象とする物理現象を代数方程式で表し、これを解いて解を求める数値解法です。解析対象領域を格子状に分割し、各格子点での差分方程式から温度や応力関数などの未知パラメータに関する連立方程式を構成します。これに境界条件を適用して解けば、格子点での関数値が求まります。図1 に有限差分法での分割例を示します。この例では、境界上の点が格子上にのらないので、特別な処理が必要になります。
有限差分法による解析の例として、図2 に示す格子状にメッシュが切られた矩形領域での温度分布を考えます。
定常熱伝導解析の問題なので、数学的には下記のポアソン方程式で表します。
ここで、λは熱伝導率、Q は単位体積当たりの熱発生率であり、T は温度です。微分方程式を代数方程式とするため、微分商を差分商へ置き換えます。1 階の微分商は、次式で表されます。
これを用いて、2 階微分商は次式となります。
上式を、Δx=Δy=h とし、x とy 方向を考慮して、ポアソン方程式に代入すると、次式となります。
領域内部の全格子点を行列的に重ね合わせ、領域境界での格子点を境界条件である温度指定・断熱・熱流速などを考慮して連立方程式を解けば、全格子点の温度が求まります。
この例で、差分法は支配偏微分方程式の微分を、差分に置き換えただけの非常にシンプルなものと捉えることができます。正確には、微分商の差分商へ置き換えを、関数T のi 点周りのテーラー展開を使って導くことで表すことができます。
差分法は、支配微分方程式の微分を直接差分で近似する方法であり、数学的にも明快で誤差評価も容易です。しかし、複雑な形状の場合は、境界の取り扱いに難点があります。複雑な境界上では、階段状の格子になってしまいます。境界上の格子間隔を小さくすることもできますが、規則的な配置とする必要があるので、境界から遠く離れた領域にも影響が及びます。
・時間領域の差分法
CAE には、時間に伴う変化を扱う解析も多くあります。例えば、動解析における過渡応答解析や、熱解析での非定常熱伝導解析などです。これらの問題の時間領域では、主に差分法が用いられます。
過渡応答解析とは、構造物の運動方程式を時刻歴で時間積分し、時間刻みでその動的な挙動を追跡する解析です。空間離散化後は微分の変数は時間だけになるので、微分は常微分になります。運動方程式は、{} で囲ったものをベクトル、[] で囲ったものを行列とすると、有限要素法では次式で表されます。
[M]、[C]、[K] はそれぞれ質量、粘性、剛性マトリックスで、F は荷重ベクトルです。一般に、荷重ベクトルは時間の関数です。上式の時間に関する数値解法には、中央差分法を初めとして、線形加速度法、ニューマークのβ法などがあります。（参考：戸川隼人、サイエンスライブラリ情報電算機33 有限要素法による振動解析、サイエンス社、1975 年）
非定常熱伝導解析とは、熱伝導問題で温度指定や対流熱伝達、熱流速などの境界条件が時間的に変化する場合や、材料定数が温度の関数である場合の解析です。非定常熱伝導方程式を時刻歴で時間積分し、その時間的な挙動を追跡します。空間離散化後の非定常熱伝導方程式は、有限要素法では次式で表されます。
ここで[C] は熱容量マトリックス、T は温度、[K] は温度剛性マトリックス、F(t) は熱的荷重項です。上式の解法には中央差分法（クランク・ニコルソン法）がよく用いられます。
2. 境界要素法
境界要素法（Boundary Element Method：BEM）とは、境界積分方程式法（Boundary Integral Equation Method：BIEM）に有限要素法的な離散化手法を組み合わせた、種々の現象の偏微分方程式を解く数値解析技法です。（参考：C.A.Brebbia 著、神谷紀夫・田中正隆・田中喜久昭共訳、境界要素法入門、培風館、1980 年）1970 年代から盛んに研究されてきました。境界要素法では、構造物の表面境界を多数の要素に分割します。要素の次元が1 つ下がって2 次元構造では線要素に、3 次元構造では面要素になります。2 次元構造での3 種類の境界要素を図3 に示しています。一定要素では、境界要素内で変位u と表面力p の変数は要素内で一定と仮定し、要素図心で代表されていると考えます。1 次要素では、変数は境界要素内で1 次分布をしており、 2 次要素では2 次分布をしていると仮定します。2 次要素の場合には、要素の図心あたりに中間節点が必要です。
・境界要素法の原理
CAE の基礎知識の第2 回 有限要素法の原理では、仮想仕事の原理について説明しました。これは、つり合い状態にある弾性体に仮想変位（δu,δv）を与えたときに生じる仮想ひずみエネルギーは、外力のなす仮想仕事に等しいというものでした。仮想仕事の原理は、下式からも誘導できます。
応力の平衡方程式× 仮想変位の領域積分式= 力学的境界条件式×仮想変位の周辺積分式
言い換えると、仮想仕事の原理を用いて問題を解くことは、応力の平衡方程式を力学的境界条件を考慮して解くことと等価です。
ブレビアによる境界要素法（ 参考：C.A.Brebbia著、神谷紀夫・田中正隆・田中喜久昭共訳：境界要素法入門、培風館、1980 年）では、上式の仮想変位を、重み付き残差法（参考：B.A.Finlayson 著、鷲津久一郎・山本善之・川井忠彦共訳：重み付き残差法と変分原理、培風館、1974 年）での重み関数と解釈して、さらに上式右辺に下式右辺第2 項を加え、基礎関係式としています。ここで、仮想表面力とは仮想変位に対応する仮想応力の表面成分です。基礎関係式は、次式で表されます。
応力の平衡方程式× 仮想変位の領域積分式
= 力学的境界条件式× 仮想変位の周辺積分式+ 幾何学的境界条件式× 仮想表面力の周辺積分式
この基礎関係式に2 回の部分積分を行い、重み関数としての仮想変位・仮想表面力に、弾性体での基本解を用います。そして、境界要素に有限要素法と同様のアイソパラメトリック定式化（離散化）を行えば、境界要素法での解くべき方程式が得られます。基本解とは、物体内の点iに作用する単位集中荷重により、点jに生ずる変位と表面力です。2次元、3 次元の等方性弾性体については、ケルビンの解として知られています。境界要素式は次式となります。
ここで{δ}、{T} はそれぞれ節点変位、節点での分布表面力ベクトルであり、いずれも境界条件による既知の値を含んでいます。{B} は体積力がある場合の荷重項です。上式で未知の値に関する項を左辺に移項し、既知の値を右辺にまとめて荷重項とすれば、最終的に解くべき連立方程式は次式となります。｛X｝は未知数で、節点変位と節点での分布表面力から成っています。
・境界要素法の長所と短所
境界要素法には、3 つのメリットがあります。
1：要素分割は領域の表面だけを行えばよく、領域内部は必要ありません。これによって問題の次元が1つ下がります。2次元問題では線要素で、3 次元問題では面要素（三角・四角形）で表面を分割すればよいので、解析モデルの作成が格段に容易となります。
2：有限要素法や差分法などの領域型解法が苦手であった無限領域問題を扱うことができます。
3：変位と応力を同時に未知数とするので、応力は変位の微分からなるひずみに比例する変位型有限要素法より応力の精度が良くなります。
反対に、境界要素法のデメリットも、3 つあります。
1：計算時間が有限要素法と比べて、長くなることです。境界要素法では一般に表面境界だけを分割すればよいので、有限要素法などに比べ節点数はかなり少なくなります。しかし、解くべき連立方程式は密行列となり、計算時間は未知数N の3 乗に比例します。一方、有限要素法では粗行列や帯行列となることが多く、計算時間は未知数N の2 乗と行列のバンド半幅NB の積に比例します。そのため、境界要素法の連立解法に要する計算時間は、有限要素法と比べて長くなる傾向にあります。
2：境界要素法では、物体力がある場合は、その積分のための図4 のような小領域への分割が必要です。また、ポストプロセッサで応力などの濃淡図を描くときも、そのための小領域への分割と、それらの頂点で応力などの値を計算することが必要となります。
3：2、3 次元弾性問題の基本解を利用するので、対応する問題が基本解のあるポテンシャル問題や2、3 次元ソリッド問題に限定されます。有限要素法のようなはり要素に対応するものはありません。
3. 有限要素法と他解法との比較
全6回で下記の内容を解説しています。
- 第1回 CAEとは？
- 第2回 有限要素法とは
- 第3回 有限要素法の収束性と誤差
- 第4回 有限要素法と他数値解法との比較
- 第5回 CAE解析の種類
- 第6回 有限要素法のテクニック
