有限要素法の定式化　アイソパラメトリック要素

アイソパラメトリック要素を用いた有限要素法について説明します。図1に示すような２次元四角形１次要素について説明します。表題を定式化としたのは，形状関数を変えることで２次要素，３次要素（これはほとんど使われていませんが）の有限要素法プログラムを書くこともできるし，３次元問題へと拡張することもできるからです。セレンディピティ族要素，適合要素・非適合要素，低減積分・完全積分についても触れております。

有限要素法詳細を復習します。要素剛性マトリクスは次式でした。

hは板厚，[B]は変位－ひずみマトリクス，[D]は応力－ひずみマトリクスです。

連立方程式を解いて節点変位が求まれば，ひずみは次式で求まります。

ひずみが求まれば，平面応力問題の応力は次式で求まります。[D]は応力－ひずみマトリクスです。

形状関数

有限要素法では未知数である変位を簡単な式（１次式か２次式）で仮定します。以下に示します。

例えば(4)式のα₀は次式となります。

かなり複雑ですね。(5)式や(6)式の場合もっと複雑になりα₀，β₀などを直接求めることは困難です。あと(1)式の積分も面倒なものになります。別の発想がないでしょうか。ここでは図1に示した四角形１次要素で考えましょう。

ξ，η直交座標系を導入し，寸法が図2のような正方形の要素を考えます。そしてi節点の変位を(u_i,v_i)，j節点の変位を(u_j,v_j)，k節点の変位を(u_k,v_k)，l節点の変位を(u_l,v_l)とします。

正方形要素内の任意の点（ξ，η）の変位ｕ（ξ，η），ｖ（ξ，η）を次式で定義します。

N_i(ξ,η)，N_j(ξ,η)，N_k(ξ,η)，N_l(ξ,η)を形状関数といいます。それでは形状関数の形を決めるわけですが，下記条件が成立すれば都合がいいですね。

なぜならば

となって，節点位置での変位を再現できるからです。

図3に示す中点m（ξ=-1,η=0）の場合を考えましょう。中点mのx方向変位が (u_i + u_l)／2，中点mのy方向変位が (v_i + v_l)／2となればなおさら都合がいいですね。なぜならば，mは中点なのだから節点mの変位は直感的に節点iと節点lの値の平均値であるべきだからです。

もうひとつ注目すべき点があります。節点mの変位が節点iと節点lだけで決まり，節点jと節点kの影響を受けないことで，隣の要素とで変位の連続性が保証されることです。図4を使って説明しましょう。位置cの変位は節点aと節点bの変位だけで決まります。一方，隣の要素の位置dの変位も節点aと節点bの変位だけで決まります。ということは位置cの変位と位置dの変位が等しくなります。つまり要素間での変位の連続性が保証されます。このような要素を適合要素といいます。

それでは，N_i(ξ,η)，N_j(ξ,η)，N_k(ξ,η)，N_l(ξ,η)をどのような形にすればよいでしょうか。次式です。

これらの式の導出過程はありません。ぱっとひらめいたものです。これらの式を形状関数として使った要素は，形状関数がぱっとひらめいて導出されたことからセレンディピティ族要素と呼ばれています。「セレンディピティ」という言葉とは，「セレンディップの3人の王子」という童話で，セレンデップ（セイロン，今のスリランカですね）の3人の王子が，その洞察力からいろいろなことを発見したことからきています。

変位－ひずみマトリクス[B]の導出

それでは２次元四角形１次要素のBマトリクスを求めましょう。Bマトリクスとは節点変位から要素内のひずみを求めるマトリクスでしたね。(2)式であらわされます。そしてひずみは次式であらわされます。

(10)式では，∂u/∂x，∂v/∂y，∂u/∂y，∂v/∂x を求める必要があります。しかし，∂u/∂ξなどは(8)式，(9)式から求められますが，xやyで偏微分したものはまだ計算できていません。

ここで，座標変換（変数変換）をします。x，y　を　ξ，η　の関数とします。次式です。

f(ξ,η)，g(ξ,η)の形を決めましょう。要素剛性マトリクス[k]の計算式を持ち出します。次式です。

この式は２重積分です。もし愚直に積分するとしたら，図5のように領域を分割して計算することになります。何とかできそうですが，例えば３次元問題では３重積分となって領域分割が複雑になりますね。後述する２次要素では要素の辺が曲線となってますます計算が複雑になりますね。

関数f(ξ,η)，g(ξ,η)を積分領域が図6のように正方形で-1≦ξ，η≦1となるようにすれば，積分が楽になり都合がいいですね。当然ですが，次式が成立しなければなりません。

(12)式の形どこかで見ましたね。(8-3)式です。x，yの座標変換を次式で行うことにします。

変位分布を仮定する式と座標変換の式を同じ形にしました。このような要素をアイソパラメトリック要素といいます。アイソ（iso）とは「同じ」という意味です。

変位　u，v　が　x，y　の関数で，x，y　が　ξ，η　の関数であるとき，u，v　の微分は次式であらわされます。

マトリクス[J]は次式で定義しました。

変数変換したときの２重積分では，dxdyがdet[J]dξdηとなることはお約束でした。det[J]はヤコビアンですね。要素剛性マトリクス[k]の積分は次式に変換されます。

積分が簡単になりました。変位－ひずみマトリクス[B]の導出に戻ります。(14)式において左側から[J]マトリクスの逆マトリクス[J]^-1を掛け算すると次式が得られます。

2行2列の行列の逆行列は簡単に計算できて，次式となります。

xとyは(13)式で示したようにξとηの関数なので，∂x/∂ξ，∂x/∂η，∂y/∂ξ，∂y/∂ηが求まります。よって，(2)式と(10)式を使って変位－ひずみマトリクス[B]が求まります。そして(16)式から要素剛性マトリクスを計算することができます。できるといってもこの積分は容易ではありません。積分は解析的に行うのではなく数値積分で行います。

ガウス・ルジャンドル (Gauss-Legendre) の数値積分公式

数値積分の方法を1次元の場合で説明します。関数f(x)の定積分は次式で近似的に計算することができます。w_iは重み，x_iは積分点で，表1のような値となります。この方法をガウス・ルジャンドルの数値積分公式といいます。

n=3の場合では，図7に示したように3点のf(x)の値を求めるだけで定積分ができます。

(19)式と表1の値を使えば，2n-1次の多項式の積分を厳密に求めることができます。2つ例を示しましょう。

表2は２次の多項式の積分です。n=2，3，4で厳密解が求まっています。n=2の場合，2n-1=2*2-1=3 なので，n=2で十分です。

表3は分子が２次の多項式，分母が１次の多項式の積分です。n=4でも厳密解が求まっていせん。分数の形では厳密解が求まらないことを覚えておいてください。

２次元問題での積分公式は次式となります。n=3を選択した場合，図8に示すように9点のf(ξ，η)を求めることで積分値が求まります。図示した点を積分点といいます。

ガウス・ルジャンドルの数値積分を使うと要素剛性マトリクス[k]の近似値は次式となります。

変位－ひずみマトリクス[B]の導出　数値積分版

数値積分する場合の変位－ひずみマトリクス[B]を求めます。ヤコビアン[J]の逆マトリクスは次式でした。この各成分を求めましょう。

(13)式をξとηで偏微分します。

ここでx_i，x_j，x_k，x_lは節点座標値であって定数なので，a₁，a₂，a₃，a₄，a₅，a₆，a₇，a₈を次式で定義します。

a₂=a₆，a₄=a₈ですね。これを使うと後で高次の項が消えますが，このままにしておきます。(13)式の偏微分は以下となります。

上式に適切なξとηの数値を代入すれば，[J]^-1が求まります。後の考察のために[J]^-1を数値化せずにξ，ηの関数として表して，[B]の成分もξ，ηの関数として表してみましょう。(18)式に(24)式を代入します。

(25)式の分母を計算しましょう。ここで，a₂=a₆，a₄=a₈を代入すると，ξηの項がなくなります。

b₁，b₂，b₃は以下のように定義しました。

(25)式に，a₂=a₆，a₄=a₈と(26)式を代入します。分母にξとηがあることに注目しておきます。

(17)式を計算しますが，その成分∂u/∂ξ，∂u/∂η，∂v/∂ξ，∂v/∂ηを計算しておきましょう。

(17)式に(28)式～(32)式を代入します。

∂u/∂x，∂u/∂y，∂v/∂x，∂v/∂yを計算していきましょう。

ξηの高次項が消えましたね。この式に(23)式を代入します。

以下に記す定数c_iを導入します。

(35)式は以下となります。同様の手順で計算すると∂u/∂y，∂v/∂x，∂v/∂yは以下となります。

(37)式を(10)式に代入すると[B]マトリクスができあがります。以下に示します。

[B]マトリクスができました。実際の計算では(39)式を計算するのではなく，ξとηを決めて数値化できる段階で早めに計算しておきます。

要素剛性マトリクス[k]の計算

[B]が求まりましたのでガウス・ルジャンドルの数値積分公式を使って[k]を計算します。でもここで疑問が生じます。積分公式では2n-1次の多項式は厳密に積分できるのですが，今回の場合nをいくらにすればよいのでしょうか。n=1，n=2，n=3の場合の積分点を図9に示します。

文献¹⁾では，「変位法による弾性問題の有限要素解析では，積分の精度が要素体積を厳密に計算するうえで十分であれば，解の収束が得られる。」とあります。しかし今は２次元問題を議論しているので要素面積となります。要素面積Aは次式で表されます。

[J]を求めましょう。(15)式に(24)式を代入します。

次はdet[J]です。(26)式と同じ形ですので次式となります。

det[J]はξ，ηの１次式なので，n=1　つまり積分点数は1点でよいことになります。後述する２次要素では det[J]はξ，ηの２次式になるので，n=2　つまり積分点数は2×2=4点となります。結論から言いますと，このようにして作られた要素は低減積分要素といって広く使われています。確かに低減積分要素では要素分割を細かくしていったら厳密解に近づきます。でも，[k]の積分はこれで十分なのでしょうか。

文献²⁾では四角形１次要素の場合n=2，四角形２次要素の場合n=3とされていることが多いとのことです。そして，このようにして計算した要素を一般に完全積分要素といいます。積分の次数がひとつ増えました。これで厳密な積分が行われるのでしょうか。少し考察しましょう。

考察の前に書いておきますが，実務では完全積分要素を使うか低減積分要素を使うかはユーザが選択できるようになっていて，完全積分要素，低減積分要素共に要素分割を細かくしていくと厳密解に近づきます。１次要素に限定すれば，完全積分要素では曲げ変形する場合ロッキングと言って剛性が高めの計算結果となり，これを回避するために低減積分要素が使われるのですが，低減積分要素ではアワーグラスモードという現象が現れ都合が悪くなることがあります。アワーグラスモードは要素分割を細かくすると回避できます。２次要素では，低減積分要素がディフォルトとなっていることが多く，ロッキングもアワーグラスモードも発生しません。

そして商用の有限要素法ソフトでは公開されていない創意工夫がなされていて，ロッキングとアワーグラスモードが現れにくい要素が実装されています。これらには図4で述べた変位の連続性を犠牲にしたものががあります。非適合要素と呼ばれていますが適合要素より計算精度がよくなっています。

それでは，(21)式の[B]^T[D][B]det[J]を計算していきましょう。繰り返しになりますが有限要素法プログラムでは[B]^T[D][B]det[J]の計算は，[B]単体，[D]単体，[J]単体で数値化しますので，以下の式の展開は不要です。nを決める考察のためだけに式を展開していきます。[B]det[J]は(39)式と(42)式から次式となります。

[D][B]det[J]は次式となります。

[B]^Tは次式となります。

[B]^T[D][B]det[J]は次式となります。

[B]^T[D][B]det[J]の1行1列の成分だけ計算しましょう。

E，ν，c_iの和と積を定数d_iで書き換えると次式になります。

以上のようになりました。[k]の各成分は，分母はξ，ηの１次式，分子はξ，ηの２次式となります。分子だけに注目するとn=2で十分なようですが，分子が１次式なので表3で述べたようにnをいくらにしても厳密な積分はできません。一般的に完全積分要素と呼んでいるのですが，「完全」な積分ではないようです。

繰り返しになりますがプログラミングでは(48)式を計算するのではなく，ξ＝±√(1/3)，η＝±√(1/3)などを[B]，[J]に代入して数値化しマトリクス計算をします。

２次要素での[k]の各成分は，分母はξ，ηの２次式，分子はξ，ηの４次式となります。分子だけに注目するとn=3（2n-1=5）で十分です。

変位の計算

やっと要素剛性マトリクス[k]が求まりました。これからは有限要素法詳細で述べたように，[k]を足し合わせて全体マトリクス[K]を構築し，境界条件に従って[K]を分解し，[K]の逆マトリクスを求め，有限要素法詳細のページの(22)式で変位を求めます。

ひずみと応力の計算

ひずみは[B]と各要素の変位ベクトルを掛け算することで求められます。次式です。

[B]はξとηの関数なのでξ＝±1，η＝±1として[B]の各成分を求め，(2)式によって各節点のひずみが求まりまるはずなのですが，一般にそのように方法はとられていません。積分点の位置つまりξ＝±√(1/3)，η＝±√(1/3)での位置でひずみが計算されています。剛性マトリクスは積分点で求めたのだからひずみも積分点で求めるのがもっとも計算精度がよいという発想です。このことを数学的に証明する力量は私にはありませんので，このような表現でとどめておきます。

積分点での応力は平面応力問題の場合，次式で求めます。

節点ひずみと節点応力は図10に示すように外挿して求めます。

有限要素法ソフトが表示する節点応力については，ＦＥＭソフトによる節点応力のページで述べます。

四角形２次要素

２次元問題の四角形２次要素は，図11に示すように辺の中点に節点があり，計8個の節点があります。

形状関数は以下の条件だと都合がいいですね。

セレンディピィティ族要素では，セイロンの3人の王子のように創意工夫でぱっとひらめかないといけないのです。次式です。

この後の手順は上述しました。

ラグランジュ族要素

セレンディピティ族要素を説明してきましたが他にもあります。ラグランジュ族要素と三角形要素族です。ラグランジュ族要素を説明します。

1次元の場合の形状関数を考えてみましょう。形状関数L_iⁿ(ξ)は図12のようにi番目の節点で1，それ以外の節点で0の値をとる必要があります。

これは次式で表されます。ラグランジュ型の補間多項式と呼ばれています。

上式はi番目の節点（ξ=ξ_i）で1となり，i番目以外の節点（ξ≠ξ_i）で0となりますね。２次元問題での形状関数は(50)式を２つ用意し掛け算することで得られます。次式となります。

２次元四角形２次要素の場合，節点のレイアウトは図13のようになり，節点iの形状関数は次式となります。(50)式，(51)式で使っているi，j，m，nと以下の図と式で使っているi，j，m，nとは，意味が全然違っていることに注意してください。

セレンディピィティ族３次要素，４次要素では形状関数を導くのは簡単ではなくセイロンの王子の出番でしたが，ラグランジュ族要素では(50)式を使って形状関数を機械的に作ることができます。

(49)式と(52)式を比較すると，ラグランジュ族要素はξ，ηの高次項ばっかりで低次項がないことがわかります。変位の分布は座標の０次，１次，２次に比例すると考えるのが自然で，セレンディピィティ族要素ではこれを表現できますが，ラグランジュ族要素では表現できませんね。また図から節点がひとつ増えています。このような理由からかラグランジュ族要素は普及しませんでした。

三角形要素族

任意の３次元形状に対して６面体要素だけで要素分割する方法はまだ見つかっていません。解析業務では四面体要素が日常的に使われています。ここでは四面体要素の２次元版の三角形要素の形状関数を紹介します。三角形要素では図14と(53)式で示すような面積座標（L_i，L_j，L_k）を導入します。３次元問題では体積座標になります。

２次元問題では独立変数はxとyの二つですが，面積座標では変数が3つになってしまいました。これを2つに減らすために(54)式を導入します。１次要素の形状関数は直感的に以下のようになります。

２次要素の形状関数は漸化式を使って求めることができ¹⁾，セイロンの王子の出番はありません。次式となります。

おわりです。

長くなりましたね。以上のように要素剛性マトリクス導出の手順は確立されています。商用有限要素法ソフトでは四面体２次要素の低減積分版が使われることが多く，六面体２次要素は要素分割な可能な形（六面体とか多角形を掃引した形）に対して使われています。そして，低減積分でかつ非適合要素がディフォルトになっていることを付け加えておきます。

参考文献

１）Ｏ．Ｃ．ツィエンキーヴィッツ，基礎工学におけるマトリックス有限要素法，吉識，山田　監訳，(S50)

２）泉，酒井，理論と実践がつながる　有限要素法シミュレーション，森北出版，(2010)