要素剛性マトリクスの導出
ＨＯＭＥ ＞ ＣＡＥ・有限要素法 ＞ 要素剛性マトリクスの導出
はじめに
図１，図２に示すような，三角形１次要素の要素剛性マトリクスを導出します。後述しますが三角形１次要素では仮定する変位分布は座標（ｘ，ｙ）の一次式で仮定します。この結果，応力とひずみは要素内で一定値をとなります。
説明を簡単にするために，二次元問題のひとつである平面応力状態で議論します。平面応力状態では６つの応力成分（σｘ，σｙ，σｚ，γｘｙ，γｙｚ，γｚｘ）のうちσｚ，γｙｚ，γｚｘがゼロとなります。
反時計回りにｉ，ｊ，ｋと節点に名前を付けます。ｉ節点について，座標を（ｘｉ，ｙｉ），作用するｘ方向の力をｆｉ，ｙ方向の力をｇｉ，ｘ方向変位をｕｉ，ｙ方向変位をｖｉとします。ｊ節点についてはそれぞれ（ｘｊ，ｙｊ），ｆｊ，ｇｊ，ｕｊ，ｖｊとし，ｋ節点についてはそれぞれ（ｘｋ，ｙｋ），ｆｋ，ｇｋ，ｕｋ，ｖｋとします。そして板厚をｈとします。
節点変位と節点に要する力の関係は次式となり，［ｋ］が要素剛性マトリクスです。［ｋ］はばね定数のようなものと解釈しております。
要素剛性マトリクス［ｋ］と，［ｋ］を合成して作られる全体マトリクス［Ｋ］には，以下の特徴があります。
- 正方マトリクスである。
- 対称マトリクスである。
- ある列について［Ｋ］の成分を合計するとゼロとなり，ある行についても［Ｋ］の成分を合計するとゼロとなる。
- ［Ｋ］の対角成分は正である。
- ［Ｋ］の逆マトリクスは存在しない。
これから，［ｋ］の具体的な値を導きます。
弾性体の支配方程式
弾性体の支配方程式は，つり合い方程式，ひずみ－変位関係式，応力－ひずみ関係式から構成されます。以下に記します。
つり合い方程式は次式です。
ここで，σｘ，σｙ，τｘｙは，座標（ｘ，ｙ）における応力成分，ｆは単位体積当たりに作用するｘ方向の外力，ｇはｙ方向の外力です。要素剛性マトリクスの導出では，つり合い方程式は使わず，その代わりに仮想変位の原理，ないしは，最小ポテンシャルエネルギの原理による式を使います。
ひずみ－変位関係式は次式です
ここで，εｘ，εｙ，γｘｙは座標（ｘ，ｙ）におけるひずみ成分，ｕ，ｖは座標（ｘ，ｙ）におけるｘ方向変位とｙ方向変位です。
応力－ひずみ関係式を説明する前に，二次元版一般化したフックの法則を説明します。ヤング率Ｅとポアソン比νの定義から，一般化したフックの法則は次式で表されます。
(4)，(5)式をσｘ，σｙ，τｘｙについて解くと次式となります。
(7)，(8)，(9)式をマトリクス表示にすると次式となります。
次式で定義されるマトリクス［Ｄ］を応力－ひずみマトリクスと呼ぶことにすると，応力ベクトル｛σ｝は(12)式で表記されます。
(10)式，(12)式が応力ーひずみ関係式です。
変位法
つり合い方程式，ひずみ－変位関係式，応力－ひずみ関係式を境界条件（変位拘束と外力）の下で解けば，未知数（σｘ，σｙ，τｘｙ，εｘ．εｙ，γｘｙ，ｕ，ｖ）を求まるはずですがこれがとても難しいので，少し簡単にします。ひずみ－変位関係式(3)式を応力ーひずみ関係式(10)式に代入し，さらにこれをつり合い方程式(2)式に代入すれば，未知数は変位ｕ，ｖだけになります。このように変位だけに関する方程式にして解くやり方は，一般に変位法と呼ばれています。しかし変位ｕ，ｖだけにしたとしても，方程式を解くにはまだまだ難しいです。ここで有限要素法の出番となります。有限要素法の基本となる式（(1)式）の未知数は変位｛ｕ｝ですので，ここで説明している有限要素法も変位法です。
変位の近似式（変位関数）
図３は表面の一部を固定し表面のある部分に外力が作用している連続体です。連続体の任意の点（ｘ，ｙ）の変位（ｕ，ｖ）を求める問題ですが，とても難しい問題であることは前述しました。ここで連続体を小さな要素に分割し，小さな要素の中では変位（ｕ，ｖ）を簡単なｘ，ｙの関数で表して近似してはどうでしょうか。近似ではありますが簡単な形なので，弾性体の支配方程式の近似解も簡単に求まりそうです。（ｕ，ｖ）の変化の激しいところでは要素分割を細かくすれば，よい近似ができそうです。近似式としては(13)，(14)式のような形があります。この近似式を変位関数といいます。
(13)式を使って要素剛性マトリクスを作ったものを三角形１次要素といい，三角形１次要素ではひとつの要素内の応力とひずみは一定値をとります，(14)式を使って要素剛性マトリクスを作ったものを三角形２次要素といい，三角形２次要素ではひとつの要素内の応力とひずみは１次式（線形）で変化します。
簡単な変位関数として(13)式，(14)式を例に出しましたが，他の形の変位関数も考えられます。変位関数に要求される条件を以下に記します。
１．真の解を十分よく近似できるものであること
２．物理的に正当でない関数を含まないこと。例えば，図３の連続な弾性体の内部で変位が不連続になるはずが
ないので，変位関数にも不連続関数が含まれていないこ。
３．幾何学的境界条件（変位の境界条件）を満たすことができること
４．各項が１次独立であること
５．微分や積分の計算が容易であること
(13)式，(14)式のような多項式の場合，上記２，４，５を満足しています。上記３については，要素の形状が凸多角形のように単純な形であれば問題ありません。今回は三角形なので問題ありません。１についてですが，多項式の次数をどんどん上げていくと真の解をより正確に表現できますので問題ありません。しかし，次数を上げると計算が複雑になるなどの他の弊害が出てくるので，次数は２次に抑えておき変化が急峻なところは要素分割を細かくする戦略をとることが多いです。
要素分割を細かくしていったら，有限要素法による計算結果が真の値に近づくかどうかには興味があります。真の値に収束する変位関数の条件は，文献3)によると以下のようになります。
１．節点変位が剛体変位の条件を満足するような場合に，要素内部でひずみを生じないこと
２．節点変位が一定ひずみ状態に与えるようなものであれば，変位関数は要素内において一定ひずみ状態を与える
ものであること
３．変位関数によって計算されるひずみは，要素と要素の間の境界面上で有限であること
(13)式，(14)式は上記条件を満足しています。
ここでは(13)式を使った三角形１次要素の要素剛性マトリクスを導出していきます。
ひずみ－変位関係式
(3)式によるとひずみは変位を偏微分したものですので，変位とひずみは密接な関係があります。変位－ひずみ関係式は次式で表されます。
ここで
［Ｂ］は変位－ひずみマトリクスといい，３行６列のマトリクスです。これから(13)式をベースとして［Ｂ］の具体的な値を求めていきます。
図２に示したｉ，ｊ，ｋ節点の座標と変位を(13)式に代入します。
(17)式のα０，α１，α２を求めます。(17)式をマトリクス表示にすると次式になります。
α０はクラメルの公式を使って次式で表されます。
同様にα１は次式で表されます。
同様にα２は次式で表されます。
(20)，(22)，(23)式をマトリクス表示にすると次式となります。
同様にして，β０，β１，β２を求めると次式となります。
ひずみεｘ，εｙ，γｘｙを求めます。(3)，(13)，(24)，(25)式からεｘ，εｙ，γｘｙは次式となります。
(26)，(27)，(28)式をマトリクス表示すると次式となります。
変位－ひずみマトリクス［Ｂ］は次式となります。
［Ｂ］マトリクスは，定数だけから構成されることに注目してください。このことは後で要素剛性マトリクスを求めるときの積分操作が簡単になります。変位関数を(14)式のような２次式とすると，［Ｂ］マトリクスの中に座標ｘ，ｙが入り，この積分操作で数値積分をする必要が生じます。
仮想仕事の原理
「ひとつの質点が，これに働くいくつかの力の作用の下でつり合い状態にあるとき，この質点に微小な仮想変位を与えても，質点に働いているすべての力がこの仮想変位によってなす仕事の総和はゼロである。」というのが仮想仕事の原理（principle of virtual work）です。式で表すと次式になります。式の導出は仮想仕事の原理を参照ください。
ここで，ｈは板厚，δｕ．δｖは仮想変位，Ｐｘ，Ｐｙは単位面積当たりの表面力，ｆ，ｇは単位体積当たりのｘ方向外力とｙ方向外力，δεｘ，δεｙ，δγｘｙは仮想変位の各成分δｕ．δｖに対応するひずみ成分，です。δεｘ，δεｙ，δγｘｙは次式で表されます。
ここで，Ｐｘ，Ｐｙ，ｆ，ｇ，δｕ．δｖ，δεｘ，δεｙ，δγｘｙ を以下のようにベクトル表示します。
また，｛δ＊｝Ｔは，｛δ＊｝マトリクス（ベクトル）の転置マトリクスを表します。｛ε＊｝Ｔも同様に転置マトリクスです。次式で表されます。
(31)式に(33)式，(34)式を代入すると，仮想変位の原理は次式で表すことができます。
仮想仕事の原理は(2)式のつり合い方程式と等価であるため，
１）つり合い方程式
２）ひずみ－変位関係式
３）応力－ひずみ関係式
の代りに，
１）仮想仕事の原理
２）ひずみ－変位関係式
３）応力－ひずみ関係式
としても同じ解が得られます。
これから仮想仕事の原理から要素剛性マトリクスを導出します。
仮想仕事の原理を用いた要素剛性マトリクスの導出
今，体積力を考えないことにすると，(35)式は次式となります。
要素に作用する表面力｛Ｐ｝は三角形要素の各辺に作用し，δｕ，δｖは三角形要素の各辺上の仮想変位ですが，連続体を要素分割したことに伴いこれらを離散化します。表面力｛Ｐ｝がなす仕事と離散化した節点力｛ｆ｝がなす仕事が等しくなるように節点力｛ｆ｝を定義します。離散化した節点での仮想変位と節点に作用している外力を(37)式のように表すことにすると，離散化した節点での仮想変位と節点に作用している外力がなす仕事は｛δ＊｝Ｔ｛ｆ｝となって，これが表面力｛Ｐ｝のなす仕事に等しいことになります。よって，(36)式の左辺は(38)式のようになります。
(36)式の右辺を変形していきます。
(12)式 ｛σ}＝［Ｄ］｛ε｝ を(36)式右辺に代入して
(15)式 ｛ε｝＝［Ｂ］｛ｕ｝ から，｛ε＊｝と｛ε｝は次式となって
(40)式を(39)式に代入して，(36)式（仮想仕事の原理）の右辺は次式となります。
ここで，転置マトリクスの公式 （［Ａ］［Ｂ］）Ｔ＝［Ｂ］Ｔ［Ａ］Ｔ を適用して
｛δ＊｝Ｔ，｛ｕ｝，ｈ は，座標（ｘ，ｙ）の関数でないため，積分の外に出して次式となります。
ここでは変位関数として(13)式を使用したことによって，［Ｂ］は(30)式で示したようにｘ，ｙの関数ではなく定数となるため，(46)式は次式となります。
ここで，Ｓは三角形の面積です。
以上が要素剛性マトリクスの導出過程です。なんか狐につままれたようですね。仮想変位の原理はちょっとピンとこないところがありますので，以下に最小ポテンシャルエネルギの原理を用いて要素剛性マトリクスを導出してみます。
最小ポテンシャルエネルギの原理
「ポテンシャルエネルギの総和（全ポテンシャルエネルギ）が最小になる状態は，力の静的なつり合い状態である。」というのが最小ポテンシャルエネルギの原理です。
最小ポテンシャルエネルギの原理は(2)式のつり合い方程式と等価であるため，
１）つり合い方程式
２）ひずみ－変位関係式
３）応力－ひずみ関係式
の代りに，
１）最小ポテンシャルエネルギの原理
２）ひずみ－変位関係式
３）応力－ひずみ関係式
としても同じ解が得られます。
特に，弾性体に外力がかかり変形していて，変形の過程において外力の大きさと向きは一定であるとき，関与するポテンシャルエネルギは以下のものとなります。
弾性体のひずみエネルギ Ｕ
外部のポテンシャルエネルギ Ｖ
全ポテンシャルエネルギΠは両者の和であり，これが最小になる状態が力の静的なつり合い状態です。
二次元問題で考えると，弾性体のひずみエネルギＵは，次式で表されます。
上式は，ばねに蓄えられるエネルギの式 Ｕ＝（１／２）ｋｘ２＝（１／２）ｋｘ・ｘ＝（１／２）力・変位 から類推できますね。
外部のポテンシャルエネルギＶは，外力のした仕事の分だけ減少することになりますのでマイナスがついて，変形前の状態を基準（原点）にすれば，次式で表されます。
ここで，ｈは板厚，ｕ，ｖは座標（ｘ，ｙ）におけるｘ方向変位，ｙ方向変位，Ｐｘ，Ｐｙは単位面積当たりの表面力，ｆ，ｇは単位体積当たりのｘ方向外力とｙ方向外力，です。
今，体積力を考えないことにすると，(49)式は次式となります。
全ポテンシャルエネルギΠは弾性体のひずみエネルギＵと外部のポテンシャルエネルギＶの和で次式で表されます。
Πが最小になる状態が静的平衡の状態であり，問題の解となります。
最小ポテンシャルエネルギを用いた要素剛性マトリクスの導出
(51)式の第１項 ひずみエネルギＵ を求めます。
今は図１，図２に示した三角形要素で，かつ(13)式で示した変位関数を用いているので，要素内の応力とひずみは座標（ｘ，ｙ）関数とはならず一定値です。よって，Ｓを三角形の面積とすれば(48)式のひずみエネルギは次式で表されます。
(52)式に(12)式と(15)式を代入すると次式になります。
ここで，転置マトリクスの公式 （［Ａ］［Ｂ］）Ｔ＝［Ｂ］Ｔ［Ａ］Ｔ を適用して
ここで，まだ［ｋ］が要素剛性マトリクスと決まったわけではないのですが，
とおくと，ひずみエネルギは次式で表されます。
(51)式の第２項 外部のポテンシャルエネルギＶ を求めます。
図１，図２に示した三角形１次要素の場合，外力は，ｆｉ，ｇｉ，ｆｊ，ｇｊ，ｆｋ，ｇｋ で，これらに対応する変位は，ｕｉ，ｖｉ，ｕｊ，ｖｊ，ｕｋ，ｖｋ なので，(51)式の外部のポテンシャルエネルギＶは次式で表されます。
(56)式，(57)式を(51)式に代入すると，全ポテンシャルエネルギΠは次式となります。
(58)式において，［ｋ］は節点座標と材料定数（Ｅ，ν）で決まり，｛ｆ｝は荷重条件として与えられるから，これらを固定して考えれば，Πは節点変位（ｕｉ，ｖｉ，ｕｊ，ｖｊ，ｕｋ，ｖｋ）によって決まります。Πの最小点を求めるにはΠの偏微分をゼロとおいて，以下に示した連立方程式を解けばよいことになります。
(59)式を計算していきましょう。Πを展開すると次式となります。
(59)式に従って，(60)式をｕｉ，ｖｉ，ｕｊ，ｖｊ，ｕｋ，ｖｋ で偏微分して，イコールゼロとおきます。
同様にして
(61)式～(66)式をマトリクス表示すると次式となり，(1)式と一致します。よって，(55)式で求まる［ｋ］は要素剛性マトリクスの各成分となります。(55)式と(47)式は一致しますので，仮想変位の原理からでも最小ポテンシャルエネルギの原理からでも，どちらからでも要素剛性マトリクスが求まりました。
参考文献
１）戸川，有限要素法概論，培風館，(S56)
２）三好，有限要素法入門，培風館，(S53)
３）Ｏ．Ｃ．ツィエンキーヴィッツ，基礎工学におけるマトリックス有限要素法，吉識，山田 監訳，(S50)
Information
仮想仕事の原理 を追加しました。
Copyright(C) ＣＡＥと強度計算の森 all rights reserved.