筆者は第54話の“CAEにモノ申す”において、有限要素法(FEM)の勉強法として、弾性学との比較検証、それに知的好奇心の涵養といったものが必要と言ってきた。そこで今回、一つの具体的問題を材料にして、その立証サンプルを提供してみたいと思う。

最初に、解析対象のモデルを紹介する。それには、FEMの泣き所の一つである無限領域の問題を採り上げる。図56-1左を見ていただきたい。これは、2次元の半無限領域の境界面に集中荷重が作用するときの応力場の問題である。いわゆる“ブシネスクの問題”と昔から呼ばれてきた問題であるが、ここでは、面外に歪が生じない平面歪問題を考える。右図はこの問題のFEMモデルである(対称性を利用した1/2モデル)。荷重値は適当に与えればいいが、以下の計算結果は、単位荷重を与えてP=1としている。

ポイント1 － 無限境界の位置

無限領域の問題では、遠方の境界をどこに設けるかで、まず悩むことになる。解析式なら距離を表す変数ｒをｒ→ ∞とすることで、事実上の無限境界の処理を施すことができるが、FEMメッシュではそうはいかない。図56-1における、dの値をどう設定するかが、知的好奇心を誘うところである－本図では、y,z軸両方向の境界位置を等距離にしているが、これは説明の便宜上のためで、もちろん違っていてもいい。念のため断わっておくと、FEMでは無限境界のために開発された無限要素という特殊要素もあるが、ここでは、標準的ユーザーを対象に話を進めているので、この要素の存在は対象外としている。

結局、経験と勘による技術的判断に任されることになるのだが、ここでは、4種の距離(本テーマでは、単位は特に意味をなさないが、一応mをイメージしてもらえればいい)を入れて解析をしてみる。その結果が表56-1である。なお、本計算時の基本メッシュサイズは、(0.1×0.1)の矩形メッシュを採用している。最下段には、理論解の値を掲載してある。

ポイント2 － 2次元弾性論の矛盾

載荷ライン上での評価点をA点にしているのは、荷重作用点直下では、応力場の特異点となるため、変位、応力ともども無限大値になり、意味をなさないからである。そのA点の沈下量を見たとき、距離dが大きくなるにつれ、ｗの値が収束の様子を示さず、増加傾向にあることが見て取れる。

これは、一見不思議な現象である。無限領域の問題を扱っているのだから、人工的に仮定する境界位置を遠くに設定すればするほど、数値の精度－精度はメッシュ密度に依存する－はともかくとして、何らかの値に収束することをわれわれは期待する。しかし、表の結果はその期待を裏切っている。

この現象を目にしたFEM初級ユーザー、または弾性学の教科書を1ページも開いたことのない人では右往左往するかもしれない。しかし、既に昔からそのことを弾性学は教えているのである。実はこの種の問題を2次元弾性論で展開していくと、矛盾を含んだ結果しか得られないのである。載荷ライン上任意点の鉛直変位を表す理論解は次の式と分かっている。

この式にある、dというのは、積分定数を決定するために、図56-1にある底辺境界を示すdと同じく、変位0の位置を仮定する距離である。dを大きく取ればとるほど、変位が大きくなることをこの式は教えている。この矛盾のため、この問題の変位場の絶対量は決められないのである。幸い、変位勾配では、dを含む定数項は消去されるので、歪あるいは応力の値は確定できることになる。

ポイント3 － mesh密度とゼロ応力

載荷ラインは構造の対称ラインでもあるので、このライン上の全ての位置で水平方向直歪とせん断歪は当然ゼロである。したがって、水平方向直応力σ_yもゼロとなる。だが、表56-1では、A点のσ_yは非零の値である。これには、2つの理由がある。一つは大きな理由で、メッシュ密度の問題である。すなわち、本問題の(0.1×0.1)のメッシュサイズでは粗すぎることを物語っている。

ところで、FEMを使用していて、理論的にはゼロになるべき応力値が、メッシュ密度を上げてもなかなか、ゼロに近づかない経験をすることがある。曲げモーメントを対象とする板やシェル構造の解析で顕著である。もっとも、この場合でも、最大値と比較すると、相対的にゼロと見なせるので、実用的には問題ないのだが。上の問題はこれと同類のことである。

二つ目の小さい理由は、節点上での応力値は通常、積分点での値を外挿計算して求められていることである。これも結局は、メッシュサイズに依存する問題なのだが、積分点での非零値を使って外挿計算するものだから、疎メッシュでは、ゼロにほど遠い値を算出してしまうことになる。

本ポイントでの話題は、経験豊かなFEMユーザーでは、改めて採り上げる内容ではないと思われるが、初級ユーザーでは、よく吟味してほしいところである。

ポイント4 － mesh密度と収束

表56-1で分かる通り、A点のσ_zとB点のσ_yの値が理論値と乖離している。特に後者の精度が悪い。B点はA点のσ_yがゼロとなるために、その代替として取った位置である。精度の悪い理由は、ポイント3で言及したメッシュサイズの問題である。そこで、メッシュサイズを細かくした計算をしてみる。

表56-1を見ると、無限境界の位置は、d=5以上のどのケースでもよさそうであるが、ここでは、一応d=10のケースで実施してみる。その結果が表56-2である。これを見ると、相当高密度のメッシュにしないと理論値に近づかないように思われるが、これは、載荷点付近の有意な応力場領域だけの問題であって、それ以外の領域では、疎メッシュのままで充分である。この辺りのことも初級ユーザーを悩ます問題であろう。

ポイント5 － 視点の変換

図56-1の問題は、本来3次元の半無限領域を平面歪モデルに簡略化して解析したため、変位の絶対量が求められないという矛盾に突き当たったわけである。ところで、平面歪モデルでの集中荷重というのは、実際は点荷重ではなく、図56-2のような、奥行きが単位幅のライン上に分布する分布荷重の総計が集中荷重Pになる問題を扱っているわけである。すなわち、本来、一直線上に等分布する線荷重の3次元問題を、平面歪問題に置き換えた際の荷重が集中荷重なのである。

ここで、再び知的好奇心が湧いてこないだろうか。それでは、3次元の半無限領域の水平境界面に点載荷する集中荷重の問題は、2次元化しては解けないなのだろうか。

そんなことはない。軸対称回転モデルとして解析すればいいのである。図56-1で使用した平面歪要素を単に軸対称要素に変換するだけでよい。ただし、荷重値はP/2→P/(2π)に変更することをお忘れなく。

メッシュサイズを(0.1×0.1)のままにして、計算した結果が表56-3である。このモデルでは、A点のσ_yが出現してくるので、B点でのそれは割愛してある。

表56-3を見ると、変位(軸対称回転モデルでは、絶対量が求まる)、応力ともd=5でほぼ値は収束していることが分かる。ただし、この値は、メッシュがあくまでも(0.1×0.1)の場合であり、解析精度云々をいう場合は、載荷点付近のメッシュ密度を上げる必要がある。

 

ついでだから、下に軸対称回転モデルの場合の、載荷ライン上の任意点の変位、応力の解析解の式を挙げておくので、FEM解析結果との比較することを、読者への宿題としておきたい。

注). 標準的なFEMでは引張応力を正値とするが、本文での数値、数式の全ては、本解析モデルが主に適用される地盤解析の慣習に合わせて圧縮応力を正値としている。

 

以上、一つのシンプルな解析モデルからでも、いろいろ知的好奇心を誘ったり、疑問が湧く話題があることを提示してきた。これらを置き去りにせず、自ら納得するまでチャレンジすることが、また、それらの積み重ねがFEM上級者への扉を開くもの、と筆者は考えている。

2009年初夏記