今、この稿を懺悔のつもりで書き始めている。懺悔の理由は後で告白する。今回のテーマはFEM 解析におけるある反力値のことである。

まず図1 を見ていただきたい。左右にある図は、2 種の平面要素タイプを使って解析した結果を示している。話を簡単にするため、4 節点四角形要素を使った1 要素モデルを考えている。下端2 節点を完全拘束し、上端2 節点に相互逆方向の同一値荷重をかけたときの下端での反力を示している。

読者は、図1 にあるモデルA、B で使用された要素タイプの違いがお分かりだろうか。そして、両モデルで反力値が荷重値とは異なる現象、さらにB では、左右点での反力値に差が出ている現象の解釈をどうされるだろうか。

一般に構造系全体が安定に支持された構造物では、その支持点での反力のうち、ある方向成分のトータル値が同方向の荷重合計値と釣り合うのは構造力学での常識だ。当然FEM 解析でも同様である。断っておくが、本稿で言う荷重、反力とは、節点個々でのそれらを言っているのではなく、それらを集計した、ある方向の合計値のことを言っている。図1 では、水平方向の荷重、反力の総計値のことである。以下この点にご注意願いたい。

モデルA では、上で言う常識に適っているが、ちょっと風変わりである。節点荷重、節点反力ともに向き合う節点間で異符号の同値を持つので、両者の釣り合いを考える以前に荷重と反力のそれぞれがゼロになっている。すなわちこの荷重は、荷重自身が自己平衡の状態になっている。しかしモデルB では、荷重が自己平衡になっているにもかかわらず、反力は、そうではない。

実は、モデルA は平面ひずみあるいは平面応力要素が使われており、モデルB は軸対称ソリッド要素が使われているのである。通常の弾性解析では、モデルA のように、荷重自身が自己平衡状態の場合、反力も自己平衡状態となるはずである。もし構造体が剛体であったならば、反力の値がどんな値でも矛盾は起こらない。弾性体の場合は、構造が変形するものだから、それに相応する応力場から等価節点力が一意に決まり、その集計結果として反力が求まることになる。それが、モデルA での反力値50.0 である。

それでは、モデルB での反力値は一体何を示しているのだろうか。図2 から明瞭なことだが、軸対称ソリッド要素の各節点は、それぞれ対称軸z に関する同心円の半径を背負っている。モデルB 下端の左右節点では、当然半径が違っているので、反力値に違いが出るのに不思議はない。しかし、それでは、半径方向の集計値(ここでは、この値をR 反力と呼ぶ)に非ゼロ値が残ることになり自己平衡状態でなくてってしまう。

ここで、懺悔告白の準備のために話を一転させる。軸対称シェル要素も含めた軸対称要素の剛性マトリックスについて少し言及しておく。通常、あらゆる要素タイプの剛性マトリックスは、歪エネルギー密度を要素体積範囲で体積分して求まるものである。軸対称要素もこの例外ではない。しかし軸対称要素の場合、その特質のためリング状の体積分となり(図2)、次式の通りとなる。

もちろん、dA=dzdr である。すなわち、体積分から面積分に落とせるものだから、平面応力や平面歪要素と同じ定式化手続きで済むわけである。さらに、上式右辺にある2π という係数は、平衡方程式右辺の荷重項にも同様に出現するものであるから、2π の乗算をすることは無駄な掛け算となり、実際の計算では両式から落としてしまえばいい。この場合、われわれは、1 ラジアン分の構造体を対象としていることに相当する。ただし市販のFEM コードの中には、2π ラジアンモデルのままで計算しているもある。

構造体を1 周のリングで考えるか、1 ラジアン分で考えるかは、荷重設定の際に注意を要する。特に集中荷重載荷の場合はそうであり、同じ線荷重強度q(たとえばN/m)を考えるにしても、前者では、P=2πrq であり、後者はP=rq を載荷しなければならない。

 

いよいよ懺悔物語の登場となる。昔、本エッセイの「第43 話・軸対称要素のドキッ」の最後で、とんでもない間違いの説明文を書いてしまった。モデルBの場合のように出てくる荷重値と反力の差を1 ラジアンモデルの側面での相互作用力であると言ってしまったのである。もちろん、これは大間違いである。もし1 周の2π ラジアンモデルを想定すれば分かる通り、相互作用力で説明できるものではない。バウムクーヘンモデルにつられてつい失言してしまったもので、全く汗顔の至りである。なお、この間違い記述は、既にホームページ版では削除してあるが、もし読者の中で、拙著「有限要素よもやま話」の書籍をお持ちの方がおられたら、その箇所を抹殺することを是非お願いしたい。この場をお借りしてお詫び申し上げたい。

ところで、それでは一体、R 反力の値は何なのか、という点に興味が湧くはずだ。それは、もちろん(1 ラジアン円弧長×r 方向幅)領域に分布するする反力値の集計値である。それが、モデルA のように何故ゼロにならないかのトリックは、軸対称要素では、必然的に円筒座標系参照の数値となっているからである。上で集中荷重のことを言ったが、R 反力はこれと同じ処理がなされた結果なのである。等価なソリッドモデル(1 ラジアンモデル)で3D 解析を実施してみれば、はっきりすることだが、同一半径上にある拘束節点の円筒座標系r 方向の反力値の集計された値が、軸対称要素で相当する拘束節点での分担反力値となり、さらに他の半径上にある拘束節点で求められた値と合計すればR 反力となるのである。

念押しとなるが、R 反力が荷重値と違いが出ていても、それはなんの不思議でもない。本来の1 周リングモデルを考えれば分かる通り、荷重は荷重自身で自己平衡状態であり、R 反力も同じく自己平衡となるからである。全く、R 反力は反力の中でも異端児なのである。

最後に、罪滅ぼしの意味で、この反力異端児について２つの話題を追記しておく。

 

その1：反力値がメッシュで変化する

結果がメッシュに依存するFEM 解析では、変位、応力がメッシュ数で変動するが、反力だけは不動のはずである。極端な話、1 メッシュモデルでさえも、荷重値と同値の反力が出てくるものである。それは、そもそもFEM 解析の究極の目的が構造物の平衡方程式を求めているからである。その解式の結果から反力が求まるのだから当然の成り行きである。たとえ粗いメッシュ分布のため精度の悪い応力がアウトプットされても、それから求まる等価節点力は-支持節点で集計した等価節点力がすなわち反力となる-全体的には外力すなわち荷重とは釣り合うからである。
しかるにこの異端児は、この原理に背くのである。荷重とは釣り合わなければならない、という軛から逃れるため応力と同じようにメッシュ依存性の性質を持つ。

 

その2 : 要素位置が軸から離れれば？

軸対称構造物が、もし対称軸から遠い位置にあったならば、どうなるか読者は想像できるだろうか。このときは、遠くなればなるほど、R 反力は荷重値に近づくことになる(もちろん符号は反対)。この理由は、軸対称構造体の歪テンソルのうち、円周方向成分 ε_θ = u/r を見れば明瞭であるが、分母のr が大きくなるにつれ、ε_θ がゼロに近づくからである。すなわち、歪場が平面歪の状態に近づくためである。

 

(注)以上の話は、荷重条件も拘束条件も構造体同様、軸対称となっていることを前提としている

2016年7月記