嶽誓㌫篇㌫篇況㌻鑛　　∠乙く＿）　　　「

有限要素法による3次元弾性体の応力集中係数

（昭和48年10月20日 原稿受理）

大学院機械工学科有 吉 一 一 第二部機械工学科光永公一

Stress Concentration Factor of Three−Dimensional

by Kazuichi ARIYOSHI

  Koichi MITSUNAGA

  This paper reports one of ways of making the analytical program for three・

dimensional problems（Axisymmetric problems）by F． E． M．．

  The stress concentration factors of the bar having a circumference notch of the

semicircle were calculated by using F． E． M．．

  The stress distributions and the strees concentration factors were compared with

the values obtained by Nishitani， Ishida and Neuber．2）

      に述べる。

 1・ は じ め に       有限要素としては3角形リング要素（図1）を  前報の2次元弾性体（有孔板）の解析1）に引き  用い，その要素内のrおよびz方向の変位μ，ρを

−・

を有限要素法によって数値計算したので・その結     （α、，α、，……α、は未定定数）

果を報告する。       と仮定する。2次元弾性体の場合と同様にして，

 また・数値計算する為のプログラムを改良した  ひずみ＿変位マトリックスを求めるわけである ので・それについても述べる。      が，ひずみは要素内の座標（r，z）に依存し要素内  この3次元弾性体の数値計算に用いたプログラ

ムは ボルト・ナットによる板締付けにおける板         Z

 その前段階として解析解の分っている半円の環 状切欠を持った丸棒の一軸引張におけを切欠底の 応力集中係数および最小断面における応力分布を

嶽誓㌫篇㌫篇況㌻鑛  ∠乙く＿）   「

当性を検討した。

 有限要素法による軸対称問題の解析方法3）4）5）は

2次元弾性体の場合とほとんど同じで，しかも数 値計算に用いるプログラムは手元にある2次元弾 性体のプログラムに少し手を加えることによって

作成できる。解析過程において異るところを簡単        図1三角形リング要素

一定とならない。そこで一番簡単な近似法とし  小断面における応力分布を求める。

て，要素の重心（ア，あを用いて要素内のひずみを   図2−1において形状・荷重および変形がr軸 一定と近似する。したがって，そのように近似す  に対して対称であるから図2−2の如くモデル化

ることによって，要素内の応力は一定となり，簡  できる。すなわち・r軸上では2方向を・2軸上 単に各要素における節点カー節点変位マトリック  ではr方向の変位を拘束する。切欠の種類として

ス （要素の剛性マトリックス）が求まる。そし  は・λ（一ρ／のを1／2・1／3・1／4・1／5の4種類と

て，各節点における外力と変位の関係が多元一次  した。なお・弾性定数Eは2・1×104kg／mm2お 方程式となる。       よびボアソン比レは1／3を用いた。

 プログラムの改良は次のようにした。現在手元   図2−2の如くモデルが決定したので要素分割 にある外カー節点変位の多元一次方程式を解くプ  を行う。ここで・モデルはz軸に対称で・しかも

ログラムは計算機内部の記憶場所を非常に大きく  有限要素を三角形リング要素としたので図2−1 占有している。外カー節点変位マトリックスが  の四半分について三角形の有限要素に分割すれば 200×200以上（節点数が100以上）になると，  よい・一例を図3に示す（図3はλ一1／3の場合 九大の計算機では処理することが不可能である  を示す）。分割時の注意すべきことは前報の平板

（ただし，内部記憶装置のみによる処理のことで  の場合1）と全く同じであるから・ここでは省略す

もっと大きなマトリックスの処理が可能なように   要素分割が終ったら一様な引張荷重をそれと等 プログラムを改良した。それについても合わせて  価な節点荷重に置き直す・この方法として・次の 報告する。       ようにした。

2． モデルおよび要素分割       r．

 図2−1に示すような直径（2の，長さ（2L），

切欠半径ρなる丸棒に一軸引張力が作用した場合 の切欠底（A点）における応力集中係数および最

  Z       Z

日！｝川日       ｝｛日 ｝1｝｝

☆ A

2■

1日い日日

  σ

r

■

0

200

150

50

0

・ ・ニゴー一→ 一一］ ・

2一ユ      2−2      

−一一【一一一司」r

    図2 モ デ ル      図3 要素分割の一例（λ＝1／3）

 図3において，載荷節点数を5コとし，その点       Start における荷重の大きさ孔五……九を      1．

      Read data

五一π（r22）2・・      「…￡…一

       ロ       わ む    の け   は      

       ロ

カーπ〔巳ち）2−（写ヨ2〕・・ 蒜、3㎜ 

    匡一2，3，4       1  ：：ぽ1y°f enti「e stif血ess

九一π〔（㌃）2一傷九）2〕・・   L−一τ一一誌＿、−en，。f，晒eS，

    但し，σは載荷面の平均応力       matrix by restraint c・・disi・n・

とした。       5． Calculation of invert matrix．

 以上により拘束条件，要素分割および節点荷重     「＿．一＿一．＿

が決ったので，これらのデータをもとにして数値    1  6 Calcul。ti。n。f displacemenL

計算を行い・応燥中徽および応力分布を求め 歴誌、乳

       Calculation of reaction and

た。      ・

 本論文では前報のように切欠底の近傍（部分要     L．一＿＿一一＿

素）を取り出しての繰返し計算は行なわなかっ       ＆  print。f calculated た．しかし，プ。グラムは部頒素を殼ての繰    「esults（δ・Rσetc）・

返し計算ができるようになっている。       St。p

 3 フローチャート      図4 フローチャート

 プログラムの主な流れは前報の2次元弾性体の

㌘籔ξ』燃漂三鯵麓慧㌢『 ・・y   ｛≡曲三ヨー

       、、

 4．計算結果および考察       ・        △は辞結果  最小断面における2方向の応力σ。の分布を求め

ておいて，そのグラフから切欠底の最大応力を求   。2。   、

める。そして，その最大応力を最小断面における

2方向の平均応力σ．で割ったものを応力集中係 数αとした。その値と西谷，石田およびNeuber らの得たα2）と比較したものを図5に示す。図

中，横軸λは切欠半径ρと丸棒の半径αとの比    t・。

（ρ／のである。この図から分るように西谷の解に ほぼ一致した。

、

、

、、

        なロ

西谷 、＼、

    Neuber   、、、、〜

   O．5      1・0   λ＝ρ／a

図5 応力集中係数

 なお，λ一1／2，1／5について，最小断面におけ   この図から，z方向の応力分布は西谷の解に近 る応力分布を西谷，Neuber2）のそれと比較し，  い結果が得られ（λ一1／2の場合はほぼ一致したと その結果を図6に示す。       みてよい），そして㊨σθの分布についてはλ一1／2

をそれぞれ最小断面における2方向の平均応力σ．  の得た分布に近い結果が得られたことが分る。

で割って無次元化し，横軸は丸棒の中心から距離   これらの計算において，計算機の記憶容量の制

rと丸棒の半径αとの比（r／のを取った。     限から総節点数を100以下とし，要素数をできる

だけ多くした場合（要素数130〜140）と少なく  が，結果はコンマ数％程の変動しか見られなかっ した場合（要素数 70〜80）の計算をしてみた  た。図5，6に示した計算結果は全要素数が70〜

2．o

ξ

、1・o

一△一本計算

一西 谷

一・一一， meuber

100の場合である。

 また，2次元（有孔板）の場合と同様にして，

部分要素を取り応力計算したが西谷らの解を大き く越へてしまってたので部分要素を取っての計算

5． む  す  び

        3次元軸対称問題の有限要素法によるプログラ

       棒の一軸引張における応力集中係数および最小断 入情     面における応力分布を計算した結果，つぎのこと

゜        ・・5 ，／。     1°     1） 部分要素を考えず，応力分布を求めておい

   図6 最小断面における応力分布      て，グラフから取って来た応力集中係数は西谷の

SUBROUTINE ASMAT

（KAKOM， SM， IS， A， NE， NOD NHEN， KZ， KZ1）

 1＝1，NOD      I＝1

い一1，NOD      J−1

lK−1，NHEN     K−1

1 L＝・1，NHEN       L＝1

KIK＝KI十K KJL＝KJ十L ISK＝IM十K JSL＝JS十L

図7  ［TSM］ の  組  立

餓かなりよ一致をみた。     であるがN・uberの働こ近い結果が与えられた・

2）劇断面における応力分布を西谷および 3）このことから・本計算に用いたプ゜グラム Neuberの得た分布と比較した結果2方向の応  が妥当なものであると考えられる。

力分布につい嘩西谷の解とよし・一致をみた・と なお・本計算は九州大学のFACOM 23°−6°お くに切欠底付近の応力勾配はよく一致している。  よび九州工大のOKITAC 4500を使用した。

そしてθおよびr方向の応力分布はλ一1／2のみ   計算に御協力頂いた，九州大学計算機センタ

SUBROUTINE RSM

（IS， A， KZ， KZ1， NOKR， NOKZ， KOZ，

NT， F， MM）

   ＿＿ ＿＿ 1ニ1，NT INDEX I＝I   l

＿＿＿＿＿＿＿＿＿1＝1，KOR

       一一一一一一一一一11＝1，KOZ

1＝4NT   I＝1

 1−1，NT     I＝1

「     IA−INDEX（1）

l       F（1）＝F（IA）

1

1       J＝1

！rl三型L−一一…一一一一一一一

1 1       JA＝INDEX J

l 8         」     K1＝1000＊1十J l l       K2＝1000＊IA十J

い   K−1，Kz K−・

U    l IS（K）−Kz Yes

ll        一一一一一

ll    L＝』呉埠．占1

1 1       1S（L）＝KI  ll  ISL＝Kl

l ｜

U   −・     1（K）−KI IS（K）−K1

L」

図 8

一および九州工大の計算機センターの方々に深く  クスの組立，拘束条件による再編成および反復法 感謝します。       による節点変位を求める計算のサプルーチンのフ       ローチャートを示す。記号の説明を省き簡単に内

      容の説明をする。      ，  前報の2次元弾性体および今回の3次元軸対称   フローチャートの説明……

問題における外カー節点変位の多元一次方程式を   ① 外カー節点変位マトリックスの組立て 解くプログラムは外カー節点変位マトリックスを      （図7，ASMAT）

そのまま2次元で記憶させている。2次元および   各要素における節点カー節点変位マトリックス 3次元問題の計算に必要なメモリー数を計算して  ［SM］を外カー節点変位マトリックス［TSM】に みる。たとえば，総節点数が100とすると，この  組み入れる。［TSM］内の0でな要素αゴ．∫をA マトリックスの大きさが200×200となり必要な  （K）に記憶させ，Z，ノをIS（K）に記憶させた。

メモリー数は40KWである。そして，応力計算  一例を表1に示す。

など前後の計算に必要なメモリー数は15KWで

あるから合わせて55KW必要である。    表1マトリックス〔TSM〕の記勧法・

九大の計撒を佑て㈱ぽ骸置のみ使 κレ5（κ）一・…×⇒ 継）⇒・

用），記憶容量が64KWであるから，節点数が   1 100前後まで計算が可能である。      2  そこで，外カー節点変位マトリックスの内部の   3 大半が◎であることから，Oでないもののみ記憶   1 させ，反復法によって多元一次方程式の解を求め   

100i ユ2055 5013 4026

 0．21388

−0．70042

−0．02971  1．25793

るようにプログラムを改良した。     ②拘束条件による［TSM］の再編成

図7・8および図9に外カー節点変位マトリッ      （図8，RSM）

       ［TSM］の再編成というのは，［TSM］マトリ       1      ックスのなかで，拘束された節点変位に対応する       行と列を取り除き，他のマトリックスの項を順次

酬＝◎      番号の若いほうへつめて行く操作のことをいう。

1＝1，NT      I＝1       この操作をIS（K）， A（K）に施す為のサブル l J＝1，NT     J＝1      一チンである。

1「一備一 }晶一一一 ﾟ蒲三．dJ≦MM         ③ 変位の計算（図9， MAT）

r− 一一西一西一西一一゜°一一

l I

幽 8  1  8  1

8 8 81    一＿＿＿＿

‖ 8

し⊥＿＿一一一＿．．一＿一一＿

民KzK1−1・蹴     たとえばタト加は願変甑の関係が次の

 「吟＝一一一    Ye§       多元一次方程式とする。

l   JS（K）＝K1

ISK＝0

1＝1，Kz     I・＝1

r香一一ひ一゜°       Yes l      IS〈D＝0

1    NN＝N 十1 I    ISNN＝ISI

1   ANN＝AI

L鞠＿● ．＿＿

KZ＝NN RETURN

01，、X1＋α、，2×2＋α1．3×3＋・・一＋α1．。X。−b1

α2，1x，＋α23x2＋02、3ズ3＋……＋α勧x。−b2

 α。．、x、＋α。．，x、＋……・・一＋α脚x。−bκ

この式を書き換えて，

x㌍一ち／輪一

ﾃ（嚥＋嚥＋一・情▲

x戸L砺／晦一

H・旦＿．＿

K＝1，NT    K＝1  1  1＝IS K）／1000  1  J＝IS（K）−1＊1000

 1     1＝J   no

l  B（1）＝B（1）／A（K）

l  X2（1）＝B（1），AA（1）＝A（K）

K＝1，NT      K＝1

1         1＝IS（K）／1000 1         J＝IS（K）−1＊1000

1＝1，KZ  I＝1      1LL＝M−l   ILL＝O  l   Xl I＝X21      PRINT

 ；  X2（1）＝B（1）      X2（1）

図9変位の計算

・………・…・ @      表2 プログラムの改良前，後のメモリー数と

．．．．．．．．．．．．．．．．＿．．      計算時間の比較

x汽／輪一

＋α励一1x。．1）

節点数100の場合

改良前1改良後

節点数200の

黶@ 合       メモリ数（kw）  55  17   ⁴⁰

とする。X、に適当な初期値を入れxl1）を求め，   計算時間（分）  1．5   7 そして，xPを川に入れxi2）を求めるという操

作を繰返して，収束したx｛K）が求める変位x、で  ぐらいに減少したからである。また，表2には節 ある。このような解法を反復法と呼び6），サブル  点数が200の場合の改良プログラムでのおよその 一チンMATでこの操作を行っている。     メモリー数を合わせて示した。これから，分るよ  この改良プログラムと改良前のプログラム必の  うに改良前に較べて節点数が200になっても楽に 要なメモリー数および計算時間の比較したものを  計算が可能である。

表2に示す。この表から，必要なメモリー数が少   ここで欠点は計算時と解の収束性の問題であ なくなったことが分る。これは外カー節点変位マ  る。計算時間はまだまだ短くすることができる。

トリックス［TSM］の記憶場所が改良前の25％  収束性の問題は現在検討中である。

        参考文献   4）会獣諸纂㌶驚素法の応用旧撒学

1）有吉，光永ほか，九州工大研究報告，26号（昭   5）有限要素法による構造解析プログラム・日本鋼構  48）       造協会編（培風館・昭45）

2） 西谷，機械学会論文集26−167（昭35）985     6）FORTRN数値計算とプログラミング・中村ほか

3） 0．C．ツィエンキヴイッツ，マトリックス有限要    （共立出版・昭47）

 素法（培風館，昭45）