弾性マトリックス - Maximaでこうぞうりきがく

構成式のお話が出たついでといっては何ですが，ここで弾性マトリックスについて触れておきたいと思います

応力-歪関係式は一般的に81個の成分を持つ4階のテンソルDを用いて次の形で表されます
T[i,j] = D[i,j,k,l].E[k,l] （i, j, k, l はそれぞれ1〜3）
2階のテンソルT，Eの独立成分をそれぞれ一次元配列t, eとして次の形に書き換えます
t[i] = D[i,j].e[j] （i, j はそれぞれ1〜9）
この9x9の正方行列Dを"Dマトリックス"あるいは弾性マトリックス（elastic matrix）と呼びます

有限要素法ではいくつか仮定を設けて，より簡単な形にしてしまいます
elastic matrix.wxm

I : 単位テンソル
E : 微小歪テンソル
E' : 工学的せん断歪
T : 応力テンソル
%i1はトレースを計算する関数を使うための宣言です（主不変量を参照ください）
微小歪テンソルEを対称テンソルとして定義します（%o4）
Eのせん断歪成分を工学的せん断歪E'（ E'[ij] = 2*E[ij] ）を使って書き直します（%o5）
線形等方性弾性体の構成式（Hook則）を%o6式に示します（フックの法則を参照ください）
成分表示すると%o7式となります（これよりTが対称であることが解ります）

ここでEの独立成分を一次元配列に並べた歪ベクトルeを%o8式に示します
（この"歪ベクトル"は，歪テンソルに単位法線ベクトルを作用させて得られるものとは異なるので注意してください）
同様にして応力ベクトルtを%o9式に示します
ラメの定数μ,λを%o10，%o11式にそれぞれ示します
応力ベクトルtを歪ベクトルeについて括った係数マトリックスをDマトリックスとして%o12式および%o13式に示します

【平面歪問題】

いまe[3]軸方向の変形が拘束されている問題を考えます（%o14）
上式を用いて応力ベクトルtを書き直したものを%o15式に示します（拘束されている自由度は省きます）
同様に，歪ベクトルeを%o16式に示します
tをeについて括った係数マトリックスを平面歪問題のDマトリックスとして%o17式に示します