有限要素法解析プログラムの中では、最終的な剛性マトリックスを得る過程で一種のフラクタル現象を見ることがある。ただし、教育目的の学習プログラムや単目的の専用プログラムではおそらく見ることはなく、各種の実務に耐える目的で開発された汎用プログラム内での話だが。

今回は、このフラクタル現象を紹介してみたいと思う。それには、下図にあるように一つのソリッド要素と3 個のデカルト座標系を用意する – 最後には、さらにもう一つの座標系も追加されるが。

まず3 座標系の説明をしておくと、次の通りである。

• X_G − Y_G − Z_Z (L_G)：構造系全体で参照される全体座標系
• X_E − Y_E − Z_E (L_E)：ソリッド要素を定義する際に使われる要素座標系
• X_M − Y_M − Z_M (L_M)：直交異方性材料を定義する材料主軸

ここで考える構造物の材料特性は等方性材料ではなく、直交異方性材料である。しかも、材料主軸が都合よく全体座標系と一致するというわけでないことを前提にする。すなわち、かなり一般的な剛性マトリックスを作成する手順を考える。まずは要素一つひとつの要素剛性マトリックスを考えてみると、どの教科書を見ても下の式が記述されているはずである。

ところが、式(1)は、等方性材料の場合はそのまま通用するのだが、異方性材料の場合では、BマトリックスもDマトリックスも同一の座標系で定義されていることが前提になる。図1 のように、材料主軸(L_M系が要素座標系L_Eと違っている異方性材料では、L_M系で定義されているDマトリックスをL_E系でのものに変換するする必要がある¹。その変換式とは次の式である。

上式のD_M、D_EはそれぞれL_M系、L_E系で定義されたDマトリックスである。T_MはDマトリックス用の座標変換マトリックスというべきもので、これについては最後でもういちど登場してもらうゆえ、中身をここでは言及しない。

さて、式(1)の要素剛性マトリックスK_Eは、構造全体に集計するに当たって、全体座標系参照の剛性マトリックスに変換する必要がある。その変換式が周知の次の式である。

ここで出てくる変換マトリックスT_Eは、これまたよく知られるL_E系とL_G系の座標軸間の方向余弦ベクトルで構成される座標変換マトリックスである。具体的内容は、式(4)中の右辺マトリックスの中身がそうであるが、記号の意味は力学や応用数学での慣習に則った記号なので、説明しなくても理解できるかと思うので省かしてもらう。

さて、式(1)から式(3)までを一つの式にまとめてみると、

となる。式(5)を眺めると、本エッセイ第51話で採り上げたマトリックス3 重積が3 度も連続して出現していることが見える。この姿を見ると、筆者はついフラクタル模様を連想してしまうのだが、読者の皆さんはどうだろうか。

まだ、続く。もし、当該要素が持つ節点が支持点にあたり、しかも全体座標系とは一致しない斜め拘束支持であったケースでは、式(5)のK_Gは極端なとき、さらに次のごとく変換されることになる。

もっとも、式(6)の表現には、ちょっと説明が必要で、この式では、要素の全自由度が斜め支持されているように取られかねないが、実際問題、そんなことはほとんどなく、一部の自由度が該当するだけだろうから、現実には、K_Gの中の小マトリックスが式(6)のような変換を受けるだけである。

 

以上、マトリックス3 重積の姿が次々湧きだしてくる様子をフラクタル模様だと感じる筆者の感想を紹介した次第だが、せっかくだから座標変換マトリックスの追記話もしておく。

式(5)に出てくる座標変換マトリックスT_Eと式(6)に出ている座標変換マトリックスT_Sは本質的に同質のものだが、式(2)にある座標変換マトリックスT_Mはこれらと違うという話である。

T_Eマトリックス、T_Sマトリックスのコンテンツは、式(4)に見る通り、方向余弦変数の1 次関数であることはよく知られていることである。しかし、T_Mマトリックスのそれは方向余弦変数の2 次関数となる。2 次関数の具体的表現は、2座標系間の歪成分の座標変換を微分のチェーンルールを使って誘導することができるが、その詳細は適当な弾性学の教科書を閲覧していただきたい。煩雑な式となる3 次元の場合を避けて、ここでは、2 次元の場合の座標変換マトリックスを式(7)に記載しておく-式(7)は工学歪を採用した変換式。

読者の皆さん、T_EマトリックスあるいはT_SマトリックスとT_Mマトリックスの本質的差異がお分かりだろうか。それは、両者で座標変換対象になるものを並べて比較すると理解の助けになるかと思う。

式(8)にある剛性マトリックスK対象の座標変換マトリックスがT_Eであり、式(9)にあるDマトリックス対象の座標変換マトリックスがT_Mであることはもちろんだが、今、K、D をそれぞれ右側にある物理量の作用素と考えてみていただきたい。そうすると、Kは節点変位ベクトルに作用した結果、やはり節点力という別のベクトルを生み出すとみなせる。また、Dは歪という2 階テンソルに作用した結果、やはり応力という別の2 階テンソルを生み出すとみなせる。すなわち、ベクトル(1 階テンソル)作用素の座標変換マトリックスは方向余弦変数の1 次関数から成り立ち、2 階テンソル作用素の座標変換マトリックスは方向余弦変数の2 次関数から成り立つと言えることになる。

実は、式(9)をもっとも一般的に表現すれば、添え字表現のテンソル式では、式(10)にようになる。

いわゆる、Dは4 階テンソルなのである。弾性体の構成則を組み立てるDマトリックスがこんな高階のテンソルなんてことには驚きだろうが、普段の解析業務で、そんなことを意識する機会はほとんどないだろう。

2015年4月記

1. 実際の定式化においては、ソリッド要素のような連続体要素の剛性マトリックスは、現在アイソパラメトリック要素が主流なので、要素座標系を経ずして、直接全体座標系に変換してしまう。ここでは、テーマの都合上、あえて一旦要素座標系に変換している。ただし、プレート/シェル要素の場合は、要素座標系への変換が必要となる。 [↩]