3-10 解説
3-10-1 数値計算方法
問題を4変形で8節点をもつ要素に分割し、それぞれの要素の形状と変位を次の関数で近似して表し、要素に作用する力は節点力に置き換え、節点力がなす仕事量が要素の歪みエネルギ-に等しいという条件で、問題に与えられた条件を満足する変位と節点力を求めている。
形状関数
X座標=N1 X1+N2 X2+N3 X3+N4 X4+N5 X5+N6 X6+N7 X7+N8 X8
Y座標=N1 Y1+N2 Y2+N3 Y3+N4 Y4+N5 Y5+N6 Y6+N7 Y7+N8 Y8
変位関数(適合要素の場合)
X方向変位=N1 δx1+N2 δx2+N3 δx3+N4 δx4+N5 δx5+N6 δx6+N7 δx7+N8 δx8
Y方向変位=N1 δy1+N2 δy2+N3 δy3+N4 δy4+N5 δy5+N6 δy6+N7 δy7+N8 δy8
変位関数(非適合要素の場合)
X方向変位=N1 δx1+N2 δx2+N3 δx3+N4 δx4+N5 δx5+N6 δx6+N7 δx7+N8 δx8
+N9 Ux1+N10 Ux2
Y方向変位=N1 δy1+N2 δy2+N3 δy3+N4 δy4+N5 δy5+N6 δy6+N7 δy7+N8 δy8
+N9 Uy1+N10 Uy2
ただし、X1~X8:8個の節点のX座標
Y1~Y8:8個の節点のY座標
δx1~δx8:8個の節点のX方向変位
δy1~δy8:8個の節点のY方向変位
Ux1、Ux2:X方向変位の非適合成分
Uy1、Uy2:Y方向変位の非適合成分
N1=(1-ξ)(1-η)(-1-ξ-η)/4
N2=(1+ξ)(1-η)(-1+ξ-η)/4
N3=(1+ξ)(1+η)(-1+ξ+η)/4
N4=(1-ξ)(1+η)(-1-ξ+η)/4
N5=(1-ξ2)(1-η)/2
N6=(1+ξ)(1-η2)/2
N7=(1-ξ2)(1+η)/2
N8=(1-ξ)(1-η2)/2
N9=(1-ξ2)ξ
N10=(1-η2)η
ξ、η:要素の自然座標(図3-6参照)
適合要素の場合は形状関数と変位関数が同じであり、このような要素はアイソパラメトリック要素と呼ばれている。
規則性のないパタ-ンで要素分割を行ったモデルに、均一応力状態になる条件を与えて、計算結果が理論値の均一な応力になるかどうかを調べるパッチテストにおいて、本ソフトでは適合要素のモデルは合格し、非適合要素のモデルは合格しない。
節点の未知の変位を未知数とする連立方程式の解法には、スカイライン法を用いており、一般に連立方程式を解く過程でエラ-は起こらない。エラ-が起こるのは、剛体変位の生じるような条件を与えた場合や、単純な入力間違いなどモデル入力デ-タのミスが原因である。
なお、計算方法の詳細は1-1-10項に示す参考文献を参照下さい。
3-10-2 節点力と実際に作用する外力
一般に分布して作用する外力を数値計算のために節点に集中して作用する節点力に置き換えている。このため計算結果の節点力から外力を調べる場合は、計算結果の節点力のままではなく、3-8-5項で説明したように節点力と外力を互いに置き換える必要がある。
置き換えは、外力の作用する箇所の要素の表面に位置を指定して、3-6節の方法で応力度を計算し、それから表面に作用する応力を求める方法による。図3-26の左は、図3-23のモデル310~312の節点力の計算結果であり、図3-26の右は3-6節の方法で応力度を計算し、外力の作用する位置の表面に垂直な方向の垂直応力度を図示したものである。
図3-26の左のようなモデルの節点力の計算結果は、作用している外力の合力と作用位置を全体的に検討する場合に使えるが、外力の分布状態を検討するには図3-26の右のように表すことが必要である。なお図3-26には生じないが、一般に表面には垂直応力度のほか、表面に平行な方向のせん断応力度も外力として作用しており、これについても検討することを忘れてはならない。
図3-26 節点力と表面の応力度
3-10-3 外力の分布状態の誤差
図3-26の左と右を比べると、問題では外力の作用していない表面にも、計算結果には応力が分布している。これは有限要素法の数値計算方法に起因する誤差で、必然的に生じるものである。この誤差は、図3
-26のモデルを比べるとわかるように、要素分割を細かくすることにより小さくできる。
図3-27は、図3-26のモデルと同じ問題(図2-12)をより細かく要素を設定したモデル320で計算した、表面の位置での表面に垂直な応力度である。A点とB点で外力が集中して作用する状態に近くなり、外力の作用しない表面で応力度が0に近づいている。
3-10-4 要素内の応力度
1つの要素内で、応力度は連続した値で分布するが、隣り合った要素間で応力度は連続していない。図3-28は、図3-18のモデル308について、黒い要素とその右隣の要素の共有する辺上の同じ位置で、それぞれの要素内で計算した主応力度である。両者は大きく異なりかなりの誤差を含んでいる。これらは3-10-1項で説明したように、各々の要素ごとに変位を簡単な関数で近似し、その偏導関数として応力を求めているため出た誤差であり、数値計算上必然的に発生するものである。
誤差を小さくするには、その部分に細かく要素を設定したモデルで計算する。また隣合った要素の計算結果の平均値で表すのが適当である。
図3-27 モデル320
図3-28 隣合った要素で同じ位置の内点の応力度
3-10-5 要素内の変位
計算結果の変位は、一つの要素内で連続した値として分布している。適合要素を用いた場合(モデル入力デ-タの表3-1のN16を1としたとき)は、隣り合った要素の接合辺上の変位は同じになる。しかし、非適合要素を用いた場合(モデル入力デ-タの表3-1のN16を0としたとき)は、隣合った要素の辺上の変位は頂点に位置する節点で同じであるが、それ以外の位置では異なる。
3-10-6 理論解の得られる問題
図3-9の均一引張応力の問題、また図3-11の純曲げ応力の問題の計算結果は、モデルの要素分割にかかわらず理論解に一致する。3-10-1項に示した変位関数により、これらの問題の変位を正確に計算できるからである。残念ながら、一般の問題の変位の分布は簡単でなく、計算結果に誤差を生じる。
3-10-7 計算結果の妥当性の検討
計算結果は作成したモデルに関して妥当でなければならない。節点の変位と節点力の計算結果を対象にして次の項目が確認できれば、計算結果は妥当であり採用できる。
1)出力デ-タのIRRの値が0である。
2)出力デ-タの後の方にある RES.N./FOR.N. = の値が十分小さい。
3)最後に出る節点力の釣り合いの誤差を表す RES.F./NOD.F. の値が小さい。
3)の節点力の釣り合いの誤差は、3-9節の単精度計算の場合は数%程度になることがあるが、通常の場合は倍精度計算をしておりこの誤差は極めて小さい。
なお、モデル入力デ-タで節点力を与えたときに、有効数字の桁数に起因する誤差のため、計算結果の節点力が与えられた節点力と若干違うことがある。3-9節の単精度計算ではこの誤差がやや大きくなる。
3-10-8 計算結果の誤差の検討
一般に理論解は分らないので、計算結果に含まれる誤差を厳密に調べることは困難である。しかし、節点力の釣り合いと問題に与えられた条件との近似性から計算結果に含まれる誤差を次のように検討することができる。
1)計算結果の節点力の釣り合いに含まれる誤差を調べる。
X、Y方向の合計の0からの誤差
X軸、Y軸に関するモ-メントの0からの誤差
2)問題に与えられた条件の変位または外力に相当する要素内の変位または応力度を3-6節のように計算し、問題に与えられた変位または荷重の分布との差を調べる。ただし、変位や応力を計算しようとしている箇所から離れた部分は除く。
2)で変位や応力を計算しようとしている箇所から離れた部分を除くことができるのは、次の理由による。
サンブナンの定理より、同じ材質形状の領域に作用する合計とモ-メントが同じで分布状態の異なる外力により生じる変形と応力は、外力が分布する区域の大きさと同じ距離以上離れた所では差が小さい。このため、ある位置の計算結果の誤差を調べる場合は、その位置から離れた所に作用する外力の影響を、サンブナンの定理を参考にして判断できる。例えば図3-26のモデル312の左の図でA点の支承反力はある領域に分布しているが、合計とモ-メントが問題における集中した反力(上向きに50)と同じであり、その分布する領域の大きさより離れている中央の荷重が作用する付近は、A点の支承反力による影響が小さい。