いくらメカに弱いドライバーであっても、車のエンジンについてのおおよそのことは知っていると想像する。同様に、有限要素法をブラックボックス利用するユーザも、有限要素法のエンジン部といってもいい剛性マトリックスのプロフィールぐらいは知っていてほしいというのが今回のテーマである。やや長くなるが、永久保存版のつもりで書き下してみた。

剛性マトリックスといっても、要素１つごとにできる要素剛性マトリックスとそれらを全要素について集計してできる全体剛性マトリックスの２つあるのは周知のとおりである。両者でその性質はほとんど同じといっていいが、１つだけ違う点がある。それについては後で述べる。

なお、ここでいう剛性マトリックスとは線形静弾性解析におけるそれで、一番基本となる係数マトリックスのことをいう。これの誘導については、それこそわんさとある参考書をご覧いただくとして、ここでは言及しない。以下は、剛性マトリックスが既にあるとしての話である。

まずは、剛性マトリックスの行列としての性格を列挙しておくと以下のとおりである。

• 対称行列である。
• 非対角項のほとんどがゼロの疎（スパース）行列である。
• 非ゼロの非対角項が対角項付近に集まる帯（バンド）行列である。
• 要素剛性マトリックスおよびそれを集成した全体剛性マトリックスは特異行列である。
• 最終的な全体剛性マトリックスは正定値行列となる。

これらのうち、１～３については行列の非ゼロ要素分布をみれば一目瞭然だし、一度はその模様をどこかで見た人も多いことだろう。ここでは、順番にワンポイント解説をしてみようと思う。

１．対称行列に関して

構造力学の教科書に“Maxwell-Betti の相反定理”というのがあったのを覚えておられるだろうか。「i 点に単位荷重が作用したときの j 点での変位は、j 点に単位荷重が作用したときの i 点での変位に等しい」という定理である。弾性体のこの性質が離散化モデルにしたときの剛性マトリックスの対称性に対応しているのである。

構造解析に限らず、他の工学問題においても、離散化したときにできる連立方程式の係数マトリックスが対称になることが多いという、解析者にとっては実にありがたい性質がある。そして、これがプログラミングにおいて、メモリの節約と数値解法の選択に効果をもたらすのである。ただし、対称性を前提にしたプログラミングは、そうでない場合に比べて（大次元マトリックスを相手にする場合には、おそらく１次元配列に展開しているから、それとも絡んで）複雑になることは否めない。

２．疎行列に関して

俗にスパースマトリックスといわれる剛性マトリックスの疎行列性は、任意の要素の挙動は隣接する要素と影響しあうだけで、遠く離れた要素には直接に影響しないという事実から由来している。誤解していけないのは、離散化モデルが常にこうなるというわけでなく、有限要素法のとる定式化がそうなってしまうということだ。境界要素法（BEM）では、逆にフルマトリックスといって非ゼロ要素の少ない密行列となってしまう。

この辺の事情は電磁気学分野の用語を利用すると分かりがいいかもしれない。電磁場理論では遠隔作用論、近接作用論という２つの立場がある。前者は万有引力の法則同様、逆自乗のクーロンの法則を基盤として、どこか力学に還元する古典的な立場であり、後者はファラデーに始まる“場”の理論であり、マクスウェルが数学的に完成させた。

大雑把な言い方になるが、数学的表現では前者が積分の形式であり、後者が微分の形式となる。そして、積分形式からスタートした定式化では係数マトリックスは一般にフルマトリックスとなり、微分形式（ここでいう微分形式とは外微分形式などといった数学分野の用語でないことに注意）からスタートした場合はスパースマトリックスとなる。

面白いことに、実際、市販の電磁場解析ソフトに両者の立場のものがあり、やはり、積分形式の方はフルマトリックスであり、微分形式の方はスパースマトリックスとなっている。直感的イメージとして下のように捉えていれば理解がしやすいのでは。

近接作用＝局所的な連続＝微分　→　有限要素法／スパースマトリックス
遠隔作用＝大域的＝積分　→　境界要素法／フルマトリックス

要は、有限要素法では、ある要素が手をつなぐのは隣の要素だけであり、一方、境界要素法のような積分形式では隣の要素も遠くの要素も同じく手をつなぐ。その結果が、剛性マトリックス内、非ゼロ要素の分布の違いとなる。

３．帯行列に関して

一般に、バンドマトリックスと呼ばれている剛性マトリックスの姿は本質的な姿ではない。節点番号の配置の仕方で左右されるものだから、その姿はユーザ次第のところもあるし、ソルバー内部での前処理としてバンド幅（帯の幅長）を縮小する機能が働いて変化することもある。

連立方程式の解法でガウスの消去法を基本とする直接法が採用されている場合、その処理時間はバンド幅の２乗に比例するので、この大きさの影響は実に大きい。したがって、直接法ではバンド幅の縮小化に気を遣う必要がある。ただ、直接法の中でもウェーブ・フロント法を利用する場合や大次元行列の解法でよく使用される反復法では、節点配置の違いからくるバンド幅の大小は処理時間に影響しないと思っていい。

バンド幅を最小化するテーマとなると、これは節点番号どうしの位相関係を扱うことになり、グラフ理論の登場を願うことになる。ただ、実用的にはそんな難しいことを考えなくても、節点番号を何らかの指標、例えば座標軸方向でソーティングしておけば、解決できることが多いものである。バンド幅縮小にこだわり過ぎて、方程式解法の処理時間は確かに小さくなったものの、逆に前処理に時間がかかってしまい、トータルではあまり効果がなかったという本末転倒なことにもなりかねないこともある。

４．特異行列に関して

特異行列というのは方程式が解けないこと、一意に解ベクトル X が決まらないような係数行列 A のことをいう。通常は行列 A の行列式 │A│ を求めると、それがゼロとなる行列のことをいう。数学的なことを言えば、行列式がゼロでも解を持つ行列もあるので、行列式 │A│＝０でもって、特異行列と定義すると数学の専門家には怒られるかもしれないが、本文の範囲では許されると思う。

特異行列の固有値を求めてみると、ゼロの固有値が出てくる。これに対応する固有モードが力学では剛体変形モードである。すなわち、弾性体に全く歪みが発生することなく変形するという、全体が剛体として変形するモードである。

有限要素法における要素剛性マトリックスおよび境界条件の処理が施される前の全体剛性マトリックスは特異行列となっている。

ここで、アナロジー的に想像してほしい。任意の関数が、直交関数の１つであるフーリエ級数の重ね合わせで表現できるように、要素の変位モードも固有ベクトル（固有ベクトル同士は直交している）の１次結合となる。このとき、固有ベクトルに剛体変形モードが入っていたら、実に都合悪い。変形が一意に決まらない。これはどうしてかといえば、境界条件が導入されていないからである。変形を一意に決める拘束条件がどこにも指定されていないからである。

そんなことで、境界条件を導入する以前の全体剛性マトリックスの方程式を解こうとしても行列式がゼロとなって解けないことになる。

５．正定値行列に関して

１～４の各性質に比べると、正定値という性質は少し分かりにくいかもしれない。そもそも正定値という数学用語からしてすっきりしない。筆者も数学者ではないので、ここではざっくばらんな言い方しかできない。

まずは定義から。

正定値行列とは（縦）ベクトル X とマトリックス A があったとき、ゼロベクトルではないどんな X を持ってきても

となる A のことをいう。この式は、ベクトルの成分 Xi で展開すると、Xi の２次式となる。したがって、正定値の概念は数学的には２次形式という数学で考察されるものである。

ところで、上式の X を有限要素の要素自由度（節点変位）ベクトルとみなすと、これは歪みエネルギーの式と同型となっていることが分かる。すなわち、歪みエネルギーは常に正の値を持つという、極めて常識的なことを言っているわけである。

ただし、１つ注意点がある。特異行列の項でも言ったように、要素剛性マトリックスでは歪がゼロの剛体変形モードが存在する。それで、要素剛性マトリックスの場合、上の式が成り立たず、≧０となってしまう。このようにゼロになる場合も含めた場合、半正定値といっている。正定値になるのは境界条件も考慮された全体剛性マトリックス（以下、Kと称する）の場合である。

正定値行列なんか持ち出して何のご利益があるのかと疑問の声が出てきそうだが、ここからが佳境である。K が正定値行列であるということから、次の事が言える。

５－１　K の対角項は必ず正値である

これは極めて常識的である。もし、対角項が負値であれば、ばねに引張荷重を掛けたら縮んでしまったということを意味してしまう。ゼロ値になることに関しては、ユーザの入力ミスで起こりうる。K の対角項がゼロになるケースというのは、内容的に次の２タイプある。

• 弾性係数や断面定数の入れ忘れという明確な入力ミス。
  対角項に限らず K の全要素はこれらの定数から構成されている。
• 不要自由度が K から消去されていない場合。
  一例を挙げると、板要素を使用した３次元問題では K は１節点６自由度として組み立てられる。このとき、面内回転自由度を持たない通常の板要素では隣接要素のない孤立した節点が存在すると、K の対角項にゼロが残ることになる（図36-2A）。たいがいのソルバーではこの不要自由度の処理はしているはずだが、局所座標系の設定などでユーザ自身が配慮する必要がある場合、ユーザには気付きにくい点でもある。
  筆者が見かけたもう一例は、梁要素の材端条件の１つのピン結合を使用する場合である。この場合、ピン結合の節点の回転剛性は、力学的処理をされて要素剛性マトリックスでのその対角項はゼロとなっている。当該節点が孤立節点であったり、回転剛性のない要素に接合されていると当然、K の対角項もゼロとなってしまう（図36-2B）。
  いずれにしても、梁要素や板要素といった回転自由度を持つ要素を使用した構造解析でよく見受けられる現象である。
   
  
  

  図36‒2　不要な回転自由度

なお、対角要素が正値ということでは、大きな落とし穴がある。正定値行列は必ず対角要素が正値であるが、対角項が正値である行列は正定値行列であるかといえば、全くそんなことはない。逆は真ならずである。

下の簡単な例を見ればいい。２つの対角項は正値となっているが、この行列は正定値行列ではない。これに関しては次項で続ける。

５－２　行を入れ替えなくてもピボットは常に正値である

有限要素法のブラックボックスユーザにとっては馬耳東風なことかもしれないが、プログラム開発者にとっては、実にありがたいのがこの性質である。ピボットというのは、日本語では枢軸要素という、いかめしい名前がついているが、連立方程式を消去法で解く際に各消去段階で使用される対角項のことである（任意の消去過程で、もはや最初の対角項の値でないことに注意）。消去法では大事な大事な項である。

そんなもの知らないという人も中学時代、一度は知らずに相手にしているのである。２元なり３元なりの連立方程式を消去法で解いた経験はあるはずだ。このとき、１変数の消去をシステマティクに処理すれば（これがガウスの消去法だが）ピボットを扱うことになる。

先ほどの２元の行列を使って、ここで消去法を実施してみよう。当然あるはずの右辺項（荷重ベクトル）は問題の本質に関係しないので、ここでは無視する。

２行１列にある非対角項2.0を消去するのに、一旦、１行１列の項（これがピボット）で１行の要素を除算する。

次に、１行を２倍して２行から引いて新しい２行とする。

見てのとおり、２番目のピボットが負値になってしまった。これが先ほど、話を留保していた件である。最初に行列の対角項がいかに正値であっても、消去段階でピボットすなわち対角項が正値でなくなってしまえば正定値行列とはいわない。

ピボットがゼロになった場合はどうなるのだろうか。上の計算途中にあったようにピボットは除算の分母に使用される。だから、それ以後の計算は不能となってしまう。これはその行列が特異行列であることを意味し、力学的には局所的、全体的にかかわらず、不安定構造になっていることを意味する。

数学としての行列なら、たとえ消去過程でピボットがゼロになっても行の入れ替え処理をすれば事なきを得ることも多々ある。

幸いなことに通常は、マトリックス K は正定値行列なのでピボットは正値となる。しかも大事な点は、最初の行列の行を入れ替えなくても済むという点である。もし、有限要素法で常にマトリックス K の行入れ替えを考慮する必要があるとなると、プログラミングは非常に複雑になり、大次元の行列が対象ともなれば、これは大変な計算コストがかかることになる。

ところで、実際の数値解析において、ピボットが完全にゼロとなることはまずない。数値誤差が入ってくるので極小さな値が出現することになる。そして、ユーザの皆さんが日ごろよく出くわす、下記のメッセージの登場となる。ソルバーによって表現は違うだろうが、言っている内容は同じである。

＊剛性マトリックスが特異（SINGULAR）となっている。
＊剛性マトリックスが ILL-CONDITIN である。

これらのメッセージが出るのは、消去過程のピボット値が非常に小さくなったことを言っている。ソルバーがなぜそんなお節介をやくのかは、安定した構造系（K が正定値行列）ではそんなことはあり得ないので、その警告のためである。

どういう場面でピボット値が小さくなるかといえば、拘束の指示を忘れて構造全体が不安定になっているケースか、ある節点での剛性が他の節点群の剛性に比べて極端に小さいという局所的不安定構造のケースである。

クイズの解答

遠隔作用論の立場にたてば、瞬間に引力がなくなるので、月は即座に飛んでいってしまう。近接作用論によれば、引力消滅の影響が局所領域間の連続で月まで届くということになるから、月が飛んでいくまでにその間の時間がかかる。現在の考えは重力場の理論が支配しているとおり近接作用論に軍杯があがるようだ。

2005年11月 記