第37話 ある疑惑
過日、ある FEM ユーザより有限要素法解析で求めた平板の曲げモーメントの値が公式集の表にある数値と違うという指摘を受けた。
その問題点というのは、3辺固定1辺自由の矩形板に等分布荷重をかけたときの固定端中央部の曲げモーメント――図37-1にある O 点の曲げモーメント―― My の値である。
実際、筆者も手元にある FEM ソルバーで計算してみたところ、たしかに、公式集にある値とは明らかに差がある。スパン長等の具体的な諸元数値はこの際、本質的ではないので記述を省くとして、そのときの結果は次表のとおりである。
O 点の My | ||||
b/a | 公式集 | ソルバー A | ソルバー B | ソルバー C |
---|---|---|---|---|
0.6 | 6515 | 6354 | 5879 | 5822 |
1.5 | 4915 | 6681 | 6185 | 6117 |
これを見て、読者はどう思われるだろうか。ソルバー間に大きな違いがないことから、プログラム内部の問題ではなく、メッシュ数に疑問が向くのが自然かもしれない。しかし、この計算時のメッシュは b/a=0.6のモデルで56×32、b/a=1.5のモデルで56×80なので充分なメッシュを張ってある。事実、これらを細分割した4倍のメッシュモデルでも結果に変化がないことは確かめている。
数値自身の差異はさておくとして、1つ、大きなミステリーがある。b/a が大きくなるにつれて FEM 結果は My が大きくなっているのに対して、公式集では逆に小さくなっている。混迷が増すところだ。
今度は、板要素から離れてソリッド要素を利用して、その応力値から曲げモーメントを計算するアプローチをとってみた。ただ、応力レベルでは当該位置は固定端で応力特異点となってしまい、はなはだ都合悪いのだが、あるメッシュで固定した b/a 対 My の傾向をみることは許されるだろう。この計算の結果は次の通りである。数値にやや違いがあるが、これは後で記すミンドリン板の結果と一致し(せん断変形が考慮されるので一致は当然)、しかも、b/a に対する傾向は平板要素の結果と同じである。
ソリッド要素の応力から求めた My | |
b/a | My |
---|---|
0.6 | 5828 |
1.5 | 6129 |
さて、これだけの材料が揃うと、疑惑の目が公式集の方へ向かってしまう。はたして、公式集の数値は正しいのだろうか。印刷ミスはないのだろうか。市販書物で掲載されている各種の境界条件での平板の計算値は大概、チモシェンコの名著“Theory of Plate and Shells(1959)”から引用されていることが多い。
この本は邦訳本もあり、昔、大学での教科書で採用されているケースも多かったのではないだろうか。そして、この本の該当ページを手繰ってみると、公式表の値とぴったり一致している。
恐れ多くも、チモシェンコの名著に疑いの目を持っていいのだろうか。しかし、この本をよくよく見ると、掲載した数表はどうも他人が計算したものを引用した節のあることが脚注にある。
そうなると、計算値というのは戦前になされたものになる。有限要素法が誕生する以前、平板理論の具体的数値解法というのは、2つしかなかった。差分式を使って微分方程式を直接解く方法か、二重フーリエ級数を利用する方法である。チモシェンコの本では後者の方法がよく採用されている。フーリエ級数の方法でも、結局は、展開項数分の係数を連立方程式で解く必要がある。
そこで、疑惑がわく。コンピュータの無かった昔、連立方程式の解式は至難の技術である。はたして、差分式だと、充分の格子を張ったのか、フーリエ級数だと充分な項数を取ったのかという疑問点である。この疑問視を裏付けるデータがある。次表を見ていただきたい。これは今回の平板モデルに対して、b/a=0.6のモデルで12×20、b/a=1.5のモデルで18×12という疎メッシュで計算した結果である。
疎メッシュでの My | |||
b/a | ソルバー A(*) | ソルバー B | ソルバー C |
---|---|---|---|
0.6 | 6332 | 5148 | 5148 |
1.5 | 6606 | 4672 | 4646 |
1つのソルバーを除いて(*)、他は b/a 対 My の傾向が公式集のものと一致した。粗いメッシュの FEM 結果と一致する公式集というのは、実は、差分式であろうと、フーリエ級数であろうと、充分な密度および項数を取っていなかったのではと思いたくなってしまうのである。
もちろん、浅学菲才な筆者には、有限要素法の結果に絶対の自信を持っているわけではない。読者の皆さんは今回の現象をどう判断されるだろうか。手元に平板を扱える解析プログラムをお持ちの読者は一度、計算を試されてはいかがだろうか。そして、何か有用な情報が得られたなら、是非、筆者までお知らせ願えないでしょうか。
ここからは、補足説明。
平板の曲げ問題を扱った公式集はもちろん、教科書、参考書の類はすべてキルヒホッフ理論をベースとした薄板を対象としたものと言っていい。したがって、有限要素法の結果をそれらと比較する場合は、ちょっと注意がいる。現在、市販 FEM コードは板要素に面外せん断変形を考慮できる厚板、すなわちミンドリン板要素を採用していることが多い。それゆえ、比較する場合は、せん断変形無視のオプションフラグを立てて実行する必要がある。
ところで、上で*と注記したソルバーでは薄板専用の板要素があり、それを実行した結果である。
ついでに言っておくと、ミンドリン板要素を使用した計算でも b/a 対 My の傾向は薄板の場合と同じであった。
2005年12月 記
読者からの寄せられたコメント
平板の最大曲げ応力の公式で、分子の部分に長手方向の長さの二乗で表示されているものは、分かりますが、短手方向の長さの二乗で表示されている理由を説明したものがありません。前者は梁の罫線ン延長線上だから分かります。後者は短手方向の引っ張りがあると思われます。ご説明下されば幸いです。
連載を終了してから長らく本サイトを覗いていなかったもので、今回久しぶりに覗いてみて初めていただいたコメントに気づきました。返信遅れて申し訳ありません。平板力学は、文字通り平面的広がりが対象の応力場の力学ですから、ご指摘の問題、なんの不思議もないと思うのですが?円板での曲げモーメントを考えれば、その特殊形としてイメージすれば納得できませんか。厳密に納得したければ、やはり数式的考察が必要になるかと思いますが、その際、「平板の基礎理論」(半谷祐彦著、彰国社)の閲覧を推奨します。本書は、もはや絶版になっているようですが、どここかの図書館にあるかと思います。比較的読みやすい本です。