有限要素法(FEM)が大衆化してしまった現在、「えっ」と思わず声を漏らしてしまいそうなユーザーの声を聞くことがある。下の質問は、FEMベースのある市販解析ソフトを使用しているユーザーからの質問を拾ったものである。

• Q1．変位が節点で出力されるのに、どうして応力は要素で出力されるのか？
• Q2．隣同士の要素の応力を節点位置でみると、同じ値でないのはどうして？

この声を聞いて、筆者は、ついにこんな質問が出てくる時代になったのだと、ため息が出てしまった。この件は、現在の主流である変位型FEMの基本中の基本である。FEM理論を学習し始めた人の質問ならば別に驚きもしないが、既に市販ソフトを使用するユーザーが発する質問であってはいけないもの、と私は感想を持つ。FEMソフトをブラックボックス利用する人たちが増えた証拠である。

これには、ポストプロセッサにも責任の一端があると言えるかもしれない。応力評価の便宜目的に、ポストプロセッサはFEMソルバーの結果を加工することが多く、FEMのオリジナル特質を隠してしまうからである。

また、この質問は、ひょっとしたら、構造解析の経験者でも骨組解析オンリーでやってきた人の中からも出てくるかもしれない。要素剛性マトリックスを生成してから後のFEMの手続きは、骨組要素でも連続体要素でも変わりはないが、要素剛性マトリックスを生成するまでの定式化では、両要素で一線を画するものがある。

質問への回答を一言で言ってしまえば、－隣接要素間で、応力の連続性は犠牲にするが、変位の連続性は考慮した、節点変位(接続要素で共通)をパラメータとする変位関数を各要素内領域で仮定するため、応力は各要素で独立に求まるから、ということになる。

 

これでは、あまりに素っ気ない説明なので、理論的背景にFEM教養も挿入しながら、以下、筆者なりの解説を試みたいと思う。ちょっと長くなるが、関心ある読者は続きをお読みください。

なお、話を簡潔にするため、しばらくは2次元弾性体を材料にして進める。初期歪、初期応力、体積力といったものは、ここでは考慮しない。すなわち、この弾性体は境界辺で規定応力と規定変位を受けて釣り合いを保っているものとする。ここの話では、この条件でも、テーマの本質を全く失うものではない。話しの展開途中で3次元弾性体に移る予定である。

数式も入れていくが、数式があるからといって、腰を引かないでいただきたい。採用する数式は話の説明上、入れた方が、理解が得やすいと思う式だけを採用し、テーマの本質にあまり影響しない式の展開は省いている。

 

平衡状態にある2次元弾性体では、下の3種類の支配方程式(x-yデカルト座標系で考慮)と境界条件が満足されていることを弾性学は教えている。

上式は弾性学の教科書の冒頭部で必ず記載されている式であり、記号の使い方も慣習的なものを使用しているので、ここでは改めて説明はしない。

これらの式を一度も見たことない、という読者がおられたら、それは、考えものですよ。きついことを言うようですが、その状態のままでは、FEMソフトを使用する資格はないと思うべきである。もし、FEMソフトを利用中であるならば、一旦休止して、序の部分でも結構だから(上に掲げた式はそれに相当)、弾性学の本を開くことを、是非勧めます。

 

さて、前置きはそれぐらいにして、本論に入りたい。式(1)～式(3)の8式を眺めると、未知変数が、変位成分2成分（ u, v ）、歪3成分（ε_x, ε_y, γ_xy ）、応力3成分（ σ_x, σ_y. τ_xy ） の計8個ある。8変数の8方程式だから、この偏微分方程式は境界条件を考慮すれば、理論的には解けるはずのものである。しかし、実在する複雑な構造物の応力解析に対して解く術をわれわれは全く知らない。偏微分方程式の解が求まるのは、シンプルな形状(矩形や円形構造)とシンプルな境界条件の場合だけである。多くの弾性学の参考書やハンドブックに掲載されている解析例はそういったものに限られているものである。この状況は、何も弾性体の応力解析に限らず、一般に物理問題で現れる偏微分方程式の境界値問題でも同様の問題である。解析的に解けるのは、シンプルな条件の場合に限られているものである。

それでは、近似的にでもいいから何とか解を求めるには、どうするかの解決法の一つが数値的解法のFEMなのである。ここでは、市販ソフトのほとんどが採用している、変位法をベースにしたFEMで話を進める。以後、式(1)～式(5)を引きずっていくのは煩雑になるので、まず、これらをマトリックス(ベクトル)表現してみる。

これらのマトリックス表現を使えば、式(1)～式(5)は次のように簡潔化された式となる。

釣合条件と変形条件を同時に満足する厳密解を求めることが困難ということで、偉い先人たちは、どちらか一方の条件を緩和するという手段で近似解を求めることを着想した。このとき、釣合条件(6)を緩めたのが変位法の出発点である。

すなわち、弾性体内の変形条件(7)と変位境界条件(9)が満足されていることを前提にする一方、釣合条件(6)の成立の厳密性を犠牲にするわけである。もちろん、釣合条件を全く無視しては何の意味もないので、なるべく成立するような手段を考えるのである。それには、次の式を考えてみる。

上式に出てくるδ記号は、変分記号を表す。“変分”て何だという初めて出会う読者は、本エッセイ中では、微分と同じような演算子と思ってくださればいいかと思う。変数一つひとつに作用するのが微分なら、変分は変数からなる関数に作用するものである。式(11)では、変位ベクトルに作用している。このδuは物理的に言えば、仮想変位と呼ばれるものである。仮想変位とは、載荷によって実際に現れる変位ではなく、荷重を受けて変形している弾性体に、さらに勝手に仮想する変位のことである。但し、与えられた拘束条件を崩すものであってはいけない(上図参照)。

式(11)の意味するところは、応力の釣合い式 ∂σ＝0 がなるべく弾性体内で成立するように、重みδuを掛けて全領域で体積分していることである。ちょっと平均的な考えが入っているわけである。これは、個々に誤差が生じるというものの、全体的に見渡して誤差が最小になる分布関数を求めるという“最小二乗法”の考えに通ずるものがある。

余談となるが、こういう平均的な考えというのは、下図に示すように、一般に対象とする分布関数に得意、不得意がある。なめらかな関数に対しては良好に振る舞う一方、局部的に変化が大きい関数は苦手とするものである。われわれが応力集中問題でメッシュ切りに苦労するゆえんでもある。

閑話休題。ところで、断わりなしに体積分という言葉を使ってきたが、先に述べたように、式(11)からは2次元から3次元に格上げする。ここからは、引用する式で2次元の簡潔さのメリットは何もないので、より一般的な3次元問題に切り替えて話を進めていくことにする。先に導入したベクトル表現の応力成分なども3次元用に切り替えれば済むことである。

さて、式(11)の意味合いをなんとか分かってもらえたとしても－いくぶん天下り的な理解でも仕方ないと思う－このままの状態では一歩も進めない。そこで、式の変形を考えていきたい。ここで、登場願うのが、有名な“ガウスの発散定理”である。体積分を面積分に置き換えるという、たいへん重宝する定理であり、物理分野における“場”の問題では広範囲に利用されているものである。この定理については、次回でも話の材料とする予定でいるので、ここでは、その公式の姿を記載するに留める。

発散定理を利用するのは、平衡状態にある弾性体を考えるとき、外力(体積力除くとして)や拘束条件といったものは、通常弾性体表面で指定されるものだから、どうしても表面に関する表示がほしいためである。

本話は教科書ではなく、あくまでもエッセイなので、式(11)の具体的展開を記述するつもりは毛頭なく、最終目標とする式の変形には、発散定理だけでなく、

という条件も使用することだけは言っておく。すると、最終的に、次の式が出てくることになる。

式(13)こそが、“仮想仕事の原理”と言われている式である。左辺は発生応力がなす仕事を、右辺は外力による仕事量を表しています。実に、この式が変位型FEMの定式化のスタート式なのである。左辺項が剛性マトリックスを生みだし、右辺項が荷重ベクトルを生み出すものである。実際のFEMの定式化においては、式(8)を式(13)に導入した次の式が使用される。

既に読者はお気付きだろうが、式(13)までは、特にFEMの話は出てきていない。そう、これまでの話はFEM理論が拠り所とする弾性体の変分原理を述べたもので、弾性体全体に当てはまる原理である。式(13)あるいは、式(14)をメッシュ分割された要素一つひとつに当てはめて進めていくのがFEMの定式化である。したがって、ここから先がFEM独特の定式化なのでる。

ところで、FEMのように人工的に弾性体を要素分割してしまうと、元の弾性体には無かった隣接要素間の境界面が出てくることになる。式(14)の由来では、そんな境界面の存在は考慮してなかったはず。この点は、次のように考える。変位型FEMでは、隣接要素間の境界面でも、やはり応力の連続性は無視して、変位の連続性のみを考慮している。しかし、この隣接要素間で共通の変位が寄与する面積分(辺積分)項が発生したとしても、われわれにとってラッキーなことがある。右図を見ていただきたい。これは、2次元の場合の例図ですが、要素間での共通辺の辺積分は、被積分項の変位が共通である上、隣接要素で積分方向がお互い反対になり(全要素について、積分方向は統一されている必要があるため)、この積分値は打ち消し合って消えてしまうのである。

さて、以上の準備を終えて、具体的なFEMの定式化に入りことにする。これからは、単なるマトリックス演算があるだけだ。式(14)では、変数に相当するのが、変位だけである(歪も変位の関数)。そこで、要素内の変位分布をその要素を構成する節点での変位をパラメータ(未知変数)として、下のように仮定する。

ｑは要素の構成節点の変位を並べた縦ベクトルであり、Ｎはその節点変位値から要素内部の変位を補間する関数で、座標変数からなる関数となる。式(15)を式(7)に代入すれば、次の変位－歪を表現するBマトリックスが求まる。

このBマトリックスこそ、FEM開発者の念願の代物であり、FEMの発展段階において、これを求めるために四苦八苦した歴史があったのである。ここでは、標準的なFEM話をしているので、ＮマトリックスからBマトリックスへの手続きが一見簡単なように見えるため-実際、教科書的要素ならそれで済むが、実務場面で使用される要素開発となると、なかなかスタンダードな手続きだけでは済まないことが多いのである。

それはさておき、式(15)、式(16)を式(14)に代入して、右辺項を左辺にシフトすれば、次式となる。

式(17)において、δｑは任意の仮想変位なので、式が常に成立するためには括弧内の式がゼロである必要がある。それで、次のようにFEM要素一つの基本方程式ができあがることになる。

式(18)を全要素について集成すれば、全体剛性マトリックスと言われる連立方程式の係数マトリックスができあがることになる。式(18)を見れば一目瞭然、未知変数は節点変位ｑである。

結構、長い道のりを経てきたが、これで、節点変位が未知数の連立方程式を解きさえすれば、節点変位が直接求まることになる。一方、応力の方は、求まった節点変位から個々の要素で式(16)を使って求めるため、各要素間で違った値になるというわけである。ついでに言っておくと、直接求まる変位の精度に比べると、1階微分の作用を経る応力の方は、一段と精度が落ちる宿命にある。

ここで、読者は疑問が湧かないだろうか。構造物の応力解析に携わる者にとって、どちらかと言えば、変位よりも応力の方が重要と考える設計者も多いことを想像する。そんな人たちからは、変位分布を仮定するよりも、応力分布こそ仮定しては、という声が聞こえてきそうである。

この考えは、確かに昔からあった。応力型FEMと言って、変位型FEMとは対照的に、釣合条件(6)と応力境界条件(10)が満足されていることを前提に、変形条件(7)の成立の厳密性を犠牲にする手段である。そして、なるべく変形条件が成立するように、下の式を考えるものである。

式(19)が応力型FEMの出発原理である。本文のところどころで、わざわざ“変位型”という用語を使用してきたのは、応力を未知数とする応力法の存在があったからである。

応力法は、骨組解析が構造解析の主流だった時代には、変位法に比べて、全体未知数が少なくて済み、結構もてていた解析法だったのだが、連続体要素を主対象とするFEMでは、その発展段階において、淘汰された感がある。それは、不静定力という明確な存在を未知数にすれば済む骨組解析に対して、連続体要素の解析では、なかなか困難な考察を必要とし、要素のタイプ毎に違ったアプローチが要求され、定式化も複雑となり、開発者からは敬遠されたからである。変位法のように全ての要素で統一化された手続きを取れないため、特にいろんな要素を扱い、汎用性を謳う流通パッケージソフトではほとんど採用されていないと思われる－ここのところ、多少、筆者の独断が入っている。

応力型FEMは汎用ソフトよりも個別目的ソフトで採用されていたと思われる。また、板要素やシェル要素のように、隣接要素間の変位連続性を満たすことが困難な要素開発では、変位法と混在した形の“ハイブリッド要素”と呼ばれる要素タイプで応力法の考えが採用されていたものである。

なお、今回の話を終えるにあたって、一つ注釈しておくことがある。変位にしても応力にしても、今まで要素境界面(辺)での連続性の条件を言ってきたが、実際問題においては、FEMの発展途上において、むしろその連続性条件を緩和した方が良好な性能を示すこともあることが発見され、現在では、この条件が金科玉条のものではなくなっていることも付け加えておく。

2009年11月記