弾性学の基本式を円筒座標系や球座標系で表現すると、軸対称問題の理論展開がすっきりするので、いろいろ興味深い問題を提供してくれることになる。そんな中の一つに、対向集中荷重を受ける中実球の弾性問題がある。

軸（上図ではＺ軸）対称の球体応力場問題では、変数ψに関係なくなるので一般解は結局、次の(球座標系形式での)ラプラス方程式の解となることが知られている。

ここで、φは応力関数と呼ばれているものだ¹ このタイプの偏微分方程式の解法は、定番の変数分離法を使用する。すなわち、φ (r,θ)を次のように r の関数と θ の関数の積の形式に仮定する。

式(2)を式(1)に代入して、ｒと θ の変数、関数を左右に分離して整理すると、

を得る。

式(3)の左辺はｒに無関係であり、右辺は θ に無関係である。ということは、両辺の式はある定数とならざるをえないわけだから、この定数を仮に λ とすれば、式(3)左辺から、次の方程式が生まれる。

また、式(3)右辺からは、次の方程式が生まれる。

つまり、一つの偏微分方程式が二つの常微分方程式に分解できたわけである。

 

式(4)の方で、 χ=cosθ , Θ(θ)=y(χ) なる変数置換を施し、λも改めて、-n(n+1)と置き直せば、ルジャンドル方程式の標準系と呼ばれている有名な変数係数の2階微分方程式の形になっていることが分かる。

この方程式の解が、数理物理方面で非常に有用なルジャンドル多項式である。既に本エッセイ第24話で紹介したように、熱伝導方程式からフーリエ級数が生まれ、ここの話のように球体の物理量を表す物理方程式(球体での静電ポテンシャル問題もしかり)からルジャンドル多項式が生まれるという、応用数学の興味深いところである。
もう一方の微分方程式(5)の方もコーシー方程式と呼ばれているそうで、その基本解は、rⁿ と r^-(n+1) であることが分かっている。それで、式(2)は次のようになる。

まあ、ここまでは一般論であり、比較的展開式がシンプルであるが、荷重を与えた具体的な境界値問題となると、煩雑な展開式が待ち構えている。一つには、荷重項による微分方程式の特解の追加、一つには、物理成分である応力成分と直結する球座標表現の応力関数が見つかっていない(筆者の知る範囲であるが)事情からある。
後者については、円筒座標系表示の応力関数が文献(1)でラブにより導かれているので、球座標系→円筒座標系の座標変換を利用して、 φを円筒座標系表示してから、各応力値を求める方法が文献(2),(3)で紹介されている。

さて、冒頭の問題の求解についてだが、有限要素法(FEM)で解析した結果が理論解とずれているという指摘をある顧客から受けた。理論解というのは、もちろん微分方程式を直接解いた方法のことである。比較対象となった文献は(4)であった。具体的には、θ=0の場合の応力 σ_θ 、すなわち、Ｚ軸直角方向の軸応力の分布のことである。
図57-2にあるように、文献(4)では、球中心から載荷点に向けてほぼ一定に近い引張りの分布を見せていた応力がr/a=0.6を過ぎるあたりから圧縮へ反転する様子が見て取れる。ところが、FEMでの結果は、圧縮／引張の反転位置がr/a=0.8以上であると言われるのである。

ところが、筆者が文献(3)にある式をプログラミングして計算したところ、下の表の結果となった。数値は無次元化した σ_θの値だが、これを見ると、文献(4)の数値計算は式(7)での多項式展開をN=3までを採用した計算のようである。

採用項数を上げていくと、引張／圧縮分布のクロス点位置はどんどん載荷点に近づく様子を示してくる。この状況はFEMの解析結果と一致していて、FEMでのメッシュ密度と理論解での多項式展開項数の多寡がよく対応していることを物語っている。

ことのついでに物理的なことを言っておけば、当問題は、載荷点付近(載荷点は応力場の特異点)の一部で大きな圧縮応力場が発生し、それと平衡をとるようにＺ軸のほとんどの領域で引張応力場となる問題のようである。

前の第55話の公式集の間違い話でもそうであったが、コンピューターが存在しなかった時代、もしくは存在していても自由に使用できなかった時代での数値計算例を参照する場合は注意が必要である。多くの文献では、計算条件の詳細を記すことは稀なので、これらの計算結果を引用するときは慎重に見なければならない。

話をここで終えてもいいのだが、せっかくだから、上で扱った問題を別のアプローチで考察する方法を考えてみることにする。

上で紹介した方法は、数学的に進める数理弾性学そのものであり、この方面が苦手な設計者には、少し難解な話だったかもしれないので、以下、工学的手段で判断することを考えてみる。

図57-3にあるような２つの弾性球の接触問題を考察材料とする。この問題は、ヘルツの接触理論として広く知られているものの一つで、この理論結果を応用するのである。

図57-3で、R₂ → ∞ および E₂ ≫ E₁ を仮定すれば、剛な平面板に載せた弾性球を押しつける問題となり、この問題は近似的に図57-1の問題と同じになるというものだ。

ところで、z軸上の σ r の分布を表す式は下の通りである(文献(3)より)。

式(8)においては、a とF_maxという二つのパラメータがあり、どちらもR1、R2、E1、E2などの関数となっている。aの方はヘルツの接触理論で最初に仮定される接触円の半径を表し(図57-3、式(9))、応力分布の重要なパラメータである。

F_maxの方は接触円での圧力分布のピーク値であり、応力値の絶対量にはもちろん関係するが、応力分布模様には関与しない。そこで、式(8)の[ ]内を計算してみれば、引張／圧縮分布のクロス点位置が a に依存していることが分かることになる。しかも、圧縮領域が載荷点より、a の値程度の小さな範囲であることが分かる。実際、式(9)に妥当な数値を代入してみると理解できるが、a の値は小さいものだ(これは、直観的にも理解できるが)。

FEMでの高密度メッシュ化、ルジャンドル多項式での高次項採用の結果を別の面から検証したことを物語っていると言える。

 

[文献]

• A.E.H.Love：A Treatise on the Mathematical Theory of Elasticity, 4th Cambridge、1927
• S.P.Timoshenko・J.N.Goodier：Theory of Elasticity,3edition, McGraw-Hill, 1934
• 中原一郎他：弾性学ハンドブック、朝倉書店、2001
• 平松良雄・岡行俊・木山英郎：非整形試験片による岩石の引張り強さの迅速試験、日本鉱業会誌、昭和40年６月

2009年初夏記

1. 数理弾性学の泰斗ラブは、彼の有名な著書[文献(1)]の中で、“応力関数”という用語を使用しているが、それを引用する他者の文献では、“歪関数”と言ったり、“変位関数”と名付けたりしている。 [↩]