FEMINGWAY 〜有限要素法解析など構造設計にまつわる数理エッセイ〜

第55話 公式集には間違いもある

本エッセイ第37話で“ある疑惑”というタイトルがあったことを読者は覚えておられるだろうか。そこでは、等分布荷重を受ける3辺固定1辺自由の矩形板の曲げモーメントの問題を話題にしていた。

図55-1 3辺固定1辺自由の矩形板

図55-1 3辺固定1辺自由の矩形板

ここで改めて、問題点を再記しておく。図55-1の原点Oでの曲げモーメントMyの値が、公式集で求めた結果とFEMソルバーでのそれとが違っているという事実である。下の表は、b/aが1.5で、分布荷重強度がp=0.00735の場合のMyの結果で、明らかに違っている。これは第37話にある計算結果と同じものである。以下、この計算条件で話を進める。

55-a

本エッセイは、この指摘を、ある土木設計者から受けたことに始まり、その解明に乗り出した筆者の顚末記である。昨年(2008年)、業務の第一線を退いたものだから、少し時間的余裕ができ、本問題にチャレンジすることができた次第である。

 

FEMソルバーの結果にばらつきがあるというものの、どれも公式集の値とはずれているところから、筆者は若干の不安点をもちながらも、当初から公式集に疑いの目を向けていた。その疑惑が確信に変わったのは、元長野高専機械工学科の風間悦夫先生が開発された、統一エネルギ原理に基づく“ノードレス要素法”という従来のFEMソルバーとは流れを異にする新しいテクノロジーを利用させてもらったときである。このソルバーは、次の値を出した。

55-b

しかし、公式集が間違いであると断定するには、上の結果は傍証でしかない。それで、とうとう公式集の方を調査することにした。本問題点を指摘された設計者が参照された公式集はおそらく、文献1) だと思う。なんとなれば、3辺固定1辺自由板の具体的数値結果表を掲載しているハンドブック的な書籍は文献1) ぐらいだからである。

ところが、文献1) は、理論的展開式の紹介には、文献2) から引用しながら、設計者が利用できる数値表は、文献4) にあるそれらを転載しているのである。

ところで、文献2) は、元ネタともいうべき同一著者による論文が既に発表されており、それが、文献3) である。実に昭和33年の論文である。

文献2) は、FEM解析が普及した今となっては、役目を果たし終えた感のある著作物であるが、出版当時では非常に重宝されたのではないかと想像するに難くない書物である。当然、絶版本であるが、筆者が大阪の中央図書館まで足を運び、目を通した結果の感想である。

閑話休題。文献2) でも 文献3) でも、設計者が利用できる形式で数値表が掲載してある。それなのに、文献1) では、それらを掲載せず、あえて文献4) のそれを転載しているのである。これは、筆者の想像だが、その理由として、文献2)、3)での数値表は、ポアソン比がゼロのものであり、かつ、たわみ量の結果が記載されていなかったからではないかと思われる。

 

さて、筆者の関心は、文献2) の理論式で計算した結果が、果たして、文献1)で掲載されている数値表に一致するのかという点であった。

よく知られているように、2重曲率を持つ板の挙動を表現する偏微分方程式は重調和型の偏微分方程式になっている。このタイプの偏微分方程式の解法には“変数分離法”という常套手段がある(本エッセイ第24話参照)。2方向に分離した変数を個別に級数展開していくわけだが、境界条件により解式の難易度が俄然違ってくる。

たわみ(w)と曲げモーメント(w″相当)がゼロの単純支持の場合、たわみをSIN関数で展開しておけば自動的に境界条件が満足されるため、比較的容易に数値計算ができる。ところが、3辺固定1辺自由といったような境界条件では一筋縄ではいかない。文献2) でもご苦労な跡が見受けられ、展開結果の式の姿を見れば、その複雑さを物語っていることが容易に理解できる。

筆者は、級数展開における未定係数を決定する連立方程式を組み立てることから、プログラミングを始めてみた。境界条件の設定からくる3種類の方程式があり、これが級数展開時の展開項の数に比例する。すなわち、項数が1だと3元の連立方程式で済むが、10項採用すると、30元の連立方程式になってしまう。とにかく煩雑な式が展開されているものだから、慎重に事を進めてプログラミングして計算した結果は? なんと、次の値が出た。

55-c

これで決まりだ。理論的解法である級数解の結果がEEM解析のそれと一致を見たわけである。すなわち、公式集にある数値表内の数値に誤りがあったのである。計算のついでに、ポアソン=0とした、文献2) にある数値表と同様の計算をしたところ、全く同じ結果となっていた。理論と数値計算でそれぞれ引用した文献を別にしたことが祟りとなって現れたわけである。

文献1) の弁護を少しすれば、転載元の文献4) にこそ間違いがあったということになる。この古典的名著の著者、チモシェンコらに計算間違いがあったのかといえば、どうもそうでもなく、この書籍に記載している数値表もSmotrovという人が計算した結果を転載していることが脚注に書いてある。差分法で計算されたその文献の年代は1936年とのことである。

コンピューターのなかった当時のこと、おそらく収束判定を厳格にしなかった格子点不足の計算結果が源流となって、転載につぐ転載の連続で今日まで続いたのだろう。それにしても、数値の間違いはともかく、b/aとMyとの関係が直感に反していることに、今までに誰か、異を唱えた人がいなかったのだろうか。それが、筆者には不思議でならない。応用力学の神様のような人の著書にけちをつける勇気ある人がいなかったというわけだろうか。

ついでだから、最後に、文献4) にある数値表の一部を、今回計算し直した結果と比較した表を下に掲げておく。

55-d

[文献]

  • 土木学会:構造力学公式集、社団法人土木学会、1986
  • 東洋一・小森清司:建築構造学体系第11 平板構造、彰国社、1970
  • 東洋一・小森清司:日本建築学会論文報告集第60号、1958
  • Timoshenko,S.P. and Woinowsky-Krieger,S:Theory of Plates and Shells McGraw-Hill、1959
    長谷川節訳:板とシェルの理論、ブレイン図書出版、1973

2009年2月記

Advertisement

読者からの寄せられたコメント

  1. 増野希陸 より:

    エッセイ55話を拝読しました.
    自分も本公式を利用したことがあるにも関わらず,疑念を抱かなかったことを恥ずかしく思います.
    また甘えてしまい申し訳ございませんが,手元の公式集を訂正するためにb/a=0.7,0.8,0.9,1.25のβ5がお分かりでしたら教えて頂けないでしょうか.
    御返事をお待ちしております.

    • Yoshiaki Harada より:

      拙著エッセイの閲覧ありがとうございます。
      さて、ご希望の件ですが、公式的数値を弾き出すには、収束性の検証が必要なため一つの解に対していくつかの計算を必要とします。
      ここ数日、出かけることが多く、早急な回答は出来ません。少し時間をもらえるなら、ご回答出来ると思います。
      しばらくお待ちください。
      from原田@京都市

      • 増野希陸 より:

        ありがとうございます.
        お忙しいところご迷惑をおかけしてしまい申し訳ございません.
        宜しくお願い致します.

        • Yoshiaki Harada より:

          FEM計算を実施したところ、以下のようになっています。最後の桁の数値はゆらぐことはあるでしょうが、上位3桁は精度が高いはずです。
          ■b/a = 0.7 の時 β5 = -0.05527
          ■b/a = 0.8 の時 β5 = -0.05590
          ■b/a = 0.9 の時 β5 = -0.05624
          ■b/a =1.25 の時 β5 = -0.05670
          やはり、b/a>1.0 のケースで誤差がひどいようですね。

          • 増野希陸 より:

            ありがとうございました.b/aで全く異なる結果になることに驚きました.公式集をすぐに訂正します.
            こんなに早く対応して頂き,ありがとうございました.今後とも本エッセイの更新を楽しみにしております.

Yoshiaki Harada へ返信する コメントをキャンセル

ページ上部へ