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

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

2次元四角形1次要素

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

要素剛性マトリクス

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

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

変位-ひずみマトリクス

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

応力-ひずみマトリクス

形状関数

有限要素法では未知数である変位を簡単な式(1次式か2次式)で仮定します。以下に示します。

変位を仮定する式

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

変位を仮定する式

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

ξ,η直交座標系を導入し,寸法が図2のような正方形の要素を考えます。そしてi節点の変位を(ui,vi),j節点の変位を(uj,vj),k節点の変位を(uk,vk),l節点の変位を(ul,vl)とします。

ξ-η直交座標系の正方形要素

正方形要素内の任意の点(ξ,η)の変位u(ξ,η),v(ξ,η)を次式で定義します。

変位u,vの定義

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

形状関数の条件

なぜならば

形状関数の条件

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

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

ξ-η直交座標系の正方形要素

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

正方形要素の変位の連続性

それでは,Ni(ξ,η),Nj(ξ,η),Nk(ξ,η),Nl(ξ,η)をどのような形にすればよいでしょうか。次式です。

2次元四角形1次要素の形状関数

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

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

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

ひずみと変位の関係

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

ここで,座標変換(変数変換)をします。x,y を ξ,η の関数とします。次式です。

変数変換

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

要素剛性マトリクス

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

要素剛性マトリクス[k]を求めるための積分

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

座標変換

関数f,gの条件

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

座標変換の式

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

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

u,vの微分

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

マトリクス[J]

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

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

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

∂u/∂x,∂u/∂y,∂v/∂x,∂v/∂y

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

[J]の逆行列

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

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

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

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

ガウス・ルジャンドルの数値積分公式

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

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

ガウス・ルジャンドルの数値積分公式の計算例

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

ガウス・ルジャンドルの数値積分公式の計算例

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

ガウス・ルジャンドルの数値積分公式

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

[k]を求める式

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

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

[J]の逆行列

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

∂x/∂ξ,∂y/∂ξ,∂x/∂η,∂y/∂ηの計算

ここでxi,xj,xk,xlは節点座標値であって定数なので,a1,a2,a3,a4,a5,a6,a7,a8を次式で定義します。

∂x/∂ξ,∂y/∂ξ,∂x/∂η,∂y/∂ηの計算

a2=a6,a4=a8ですね。これを使うと後で高次の項が消えますが,このままにしておきます。(13)式の偏微分は以下となります。

∂x/∂ξ,∂y/∂ξ,∂x/∂η,∂y/∂ηの計算

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

[J]の逆行列

(25)式の分母を計算しましょう。ここで,a2=a6,a4=a8を代入すると,ξηの項がなくなります。

[J]の逆行列の分母

b1,b2,b3は以下のように定義しました。

[J]の逆行列の分母

(25)式に,a2=a6,a4=a8と(26)式を代入します。分母にξとηがあることに注目しておきます。

[J]の逆行列

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

[B]の導出

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

[B]の導出

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

[B]の導出

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

[B]の導出

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

ciの定義

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

[B]の導出

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

[B]の導出

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

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

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

nと積分点

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

要素面積を求める式

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

[J]

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

det[J]

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

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

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

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

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

[B]detJ

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

[D][B]det[J]

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

[B]T

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

[B]T[D][B]det[J]

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

[B]T[D][B]det[J]11

E,ν,ciの和と積を定数diで書き換えると次式になります。

[B]det[J]ij

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

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

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

変位の計算

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

ひずみと応力の計算

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

変位-ひずみマトリクス

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

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

変位-ひずみマトリクス

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

節点ひずみと節点応力の外挿

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

四角形2次要素

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

四角形2次要素

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

2次元四角形2次要素の形状関数

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

2次元四角形2次要素の形状関数

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

ラグランジュ族要素

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

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

ラグランジュ要素族の形状関数 1次元

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

ラグランジュ型の補間多項式

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

ラグランジュ族要素の形状関数

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

ラグランジュ族要素の形状関数 2次元

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

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

三角形要素族

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

面積座標

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

1次元三角形1次要素の形状関数

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

2次元三角形2次要素の形状関数

おわりです。

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

参考文献

1)O.C.ツィエンキーヴィッツ,基礎工学におけるマトリックス有限要素法,吉識,山田 監訳,(S50)

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