有限要素法における 2 次元弾性体要素の精度に関する研究

建設構造研究室 渡邊 章吾
指導教官 岩崎 英治
長井 正嗣

1.諸言
有限要素法では構造物を解析する際、部材が非圧縮状態に近くなったり、非常に細長くなると、極度に精度が落ちることがある。
これはロッキングと呼ばれる現象である。

ルギーモードの数と形状を導く。また、低減積分時における、ゼロエネルギーモードの発生状況を調べる．
2.ロッキングと剛性行列の特異性
2.1 体積ロッキング

ロッキングは有限要素のつりあい式となる剛性方程式の剛性を過大評価してしまうことによって生じる精度低下である。このロッキングについても大きく 2 つに分けることができ、一つは部材の非圧縮性によって起こる体積ロッキング、もう一つは、非常に細長い梁を解析した時にせん断ひずみが拘束され発生するせん断ロッキングである。
このような精度低下を避けるために低減積分法という積分方法がある。これは、解析の際、剛性行列の積分点よりも一つ少ない積分点で積分する方法である。低減積分法
によって、うまく非圧縮による拘束や、せ図 1 解析対象とする片持ちはり

最初に、図 1 のような片持ち梁を対象に，非圧縮性によるロッキングを示す。ここでは，4 節点四辺形要素、8 節点四辺形要素、9 節点四辺形要素を使用し、ポアソン比を 0〜0.49999 まで変える。要素分割を図 2 に示す。なお、ヤング係数は 2×103kN/cm2とした。図 3 にその結果を示す。横軸はポアソン比、縦軸は(計算値/正解)である。

ん断ひずみの拘束を回避できるが、このような低減積分を行ったときには、ゼロエネ
ルギーモードと呼ばれるひずみエネルギー(a) 4 節点四辺形要素に含まれない変形モードを生じることがある。このようなゼロエネルギーモードを含んだ剛性行列は、解が得られなかったり、(b) 8 節点四辺形要素解が得られても、特異な変形モードが生じたりする。そこで，本研究では、4 節点、8節点、9 節点四辺形要素において、各要素の変形モード毎のひずみを調べ、ゼロエネ(c) 9 節点四辺形要素
図 2 要素分割ここで、ポアソン比が 0.5 に近づくという1.2(計算値/正解)

1

ことは、ラメの定数λが限りなく大きくな

0.84節点要素完全積分
8節点要素完全積分
9節点要素完全積分

0.60.4

るということにつながる。
これにより、ひずみエネルギーの停留条
件に、(εxx+εyy)=0 という非圧縮条件が付

0.2000.1

0.2

0.30.4

け加えられることとなり、誤差が生じる。

0.5

ポアソン比

これが体積ロッキングである。

図 3 モーメント荷重が作用した

2.2 せん断ロッキング

梁のたわみ＜完全積分＞

非常に細長い部材に、2 次元弾性体要素
なお、剛性行列は、Gauss-Legendre の数を用いた場合などに、著しく精度が悪化す値積分法を用いて、正確に積分している。

ることがある。体積ロッキングの時と同様

また、右端に 5kN の鉛直荷重を作用させたの 片 持 ち 梁 に お い て ， 部 材 長 を 12 〜

ときの結果を図 4 に示す。

1200cm まで変化させてモーメント荷重、及び鉛直荷重を加えたときのたわみを調べ1.2てみた。ポアソン比には，0.3 を用いる．結(計算値/正解)

1

果を図 5,6 に示す。

0.8
4節点要素完全積分
8節点要素完全積分
9節点要素完全積分

0.60.4

荷1.2

わ

完 積分

算

わ

10.2

00.1

0.2

0.30.4(計算値/正解)

00.5

ポアソン比

図 4 鉛直荷重が作用した

0.8
4節点要素完全積分
8節点要素完全積分
9節点要素完全積分

0.60.40.2
0

梁のたわみ＜完全積分＞

10100

1000
部材長(L)

10000100000

図 5 モーメント荷重を作用したモーメント荷重を作用したモーメント荷重、鉛直荷重を作用させたもの、どちらにおいてもポアソン比が、0.5

鉛直荷1.2に近づくほど精度が悪化している。また、4梁のたわみ＜完全積分＞
わ完 積分
算
わ2 次元弾性体のひずみエネルギーは次のような式で与えられる。
U={}1(λ + G )(ε xx + ε yy ) 2 + G (ε xx − ε yy ) 2 + 4Gε xy2 dA2 ∫A(計算値/正解)

1

節点要素の場合は非常に誤差が大きい。

0.8
4節点要素完全積分
8節点要素完全積分
9節点要素完全積分

0.6
0.40.2
0
10100

1000
部材長(L)

10000100000

ラメの定数λは、次のように表される。

λ=νE(1 + ν )(1 − 2ν )図 6 鉛直荷重を作用した梁のたわみ＜完全積分＞

この場合もまた、4 節点四辺形要素のたわみにおいては、極度の精度低下が起きていε xx =る．

UU1+η 3aa2ε xy =ε yy =

VV2+ξ 3bbU
V
U 2 V1+ +ξ 3 +η 3babaここで、2 次元の細長い部材ではせん断ひずみがゼロという拘束条件

2ε xy ≅ 0が課され、この条件により発生可能なひずみは
図 7 長方形要素

ε xx =

U1aε yy =V2bこのような精度低下について、図 7 のよのみとなる。図 8 の変形モードより、Ｕ1,V2

うな、長方形要素を用いて考察する。x と yは、部材の伸縮、U3,V3 は、曲げ変形を表し方向の長さが、それぞれ、2a、2b の 4 節点ていることが分かる。上式により、梁のよ長方形要素のひずみを、陽な形式で表す．

うな細長い部材を，長方形断面要素を用いこのとき、座標(x,y)と向きが一致し，要素て解析する場合には、部材の伸縮変形は、中央に原点を持つ−1〜1 の範囲の値をと可能であるが、曲げ変形が拘束され、支持る正規座標 (ξ ,η ) を導入する。

条件や、荷重条件によっては、精度が非常

ここで、各節点の x 及び y 方向変位 u,v は、に悪くなることがある。これが、せん断ロベクトル表示で表すと、正規座標 ξ ,η から、ッキングの原因である。

次のようにして求められる。

u = sU 0 + ξ U 1 + η U 2 + h U 3v = sV + ξV1 + ηV 2 + hV3

2.3 低減積分法
低減積分法とは、剛性行列を求める際に，また、これらのＵi,Vi(I=0,1,2,3)は、変形全てのひずみ成分について、1 点少ない積量(剛体回転を含む)を表している。これらの分点数を用いて積分を行うことである。一諸量に対応するモードを図 8 に示す。

方，選択的低減積分は、非圧縮やせん断ひずみがゼロとなるような拘束の課せられるひずみ成分だけを低減積分し、残りの項は、正確な積分を行う方法である。せん断ロッキングが起こる問題では、せん断ひずみのみ低減積分し、体積ロッキングが起こる問図 8 剛体モードと変形モード

題では、体積ひずみのみ低減積分する。先の計算例と同じ片持ち梁に，低減積分法と2 次元弾性体要素のひずみ式は次のように選択的低減積分法を用いた結果を、図 9〜表される。

16 に示す。

1.61.31.5

1.21.41.1

4節点要素低減積分
8節点要素低減積分
9節点要素低減積分

10.9(計算値/正解)(計算値/正解)

1.4

1.10.90.70.80.6

1000.10.20.3
ポアソン比

0.4

1000.51000
部材長(L)

10000100000

図 9 モーメント荷重作用

図 13 モーメント荷重作用(体積ロッキング問題)
体積ロッキング問題)＜低減積分＞(せん断ロッキング問題)
せん断ロッキング問題)＜低減積分＞
1.6

1.31.51.2
1.1

4節点要素低減積分

18節点要素低減積分

0.9

9節点要素低減積分(計算値/正解)

1.4

低減積分 (計算値解)

1.4
1.3
4節点要素低減積分
8節点要素低減積分
9節点要素低減積分

1.2
1.1
10.8

0.9

0.7

0.8

0.6
00.1

0.2
0.3
ポアソン比

0.4

100.5

100

1000
部材長(L)

10000

100000

図 10 鉛直荷重作用

図 14 鉛直荷重作用(体積ロッキング問題)
体積ロッキング問題)＜低減積分＞(せん断ロッキング問題)
せん断ロッキング問題)＜低減積分＞

1.4

1.2

1.3

1.15

1.2

4節点要素選択的低減積分
8節点要素選択的
低減積分
9節点要素選択的
低減積分

1.1
10.9
0.8

(計算値/正解)(計算値/正解)

4節点要素低減積分
8節点要素低減積分
9節点要素低減積分

1.2

10.8(計算値/正解)

1.3

選択減積

算

解

1.1
4節点要素選択的低減積分
8節点要素選択的低減積分
9節点要素選択的低減積分

1.05
10.95
0.9
0.85

0.7

0.80.6
00.1

0.2
0.3
ポアソン比

0.4

100.5

100

1000
部材長(L)

10000

100000

図 15 モーメント荷重作用(体積ロッキング問題)
体積ロッキング問題)＜選択的低減積分＞(せん断ロッキング問題)<
せん断ロッキング問題)<選択的低減積分)<選択的低減積分>選択的低減積分>

1.4

1.2

1.3

1.15

1.2
4節点要素選択的低減積分
8節点要素選択的
低減積分
9節点要素選択的
低減積分

1.1
1
0.9
0.8

(計算値/正解)(計算値/正解)

図 11 モーメント荷重作用

選択減積

算

解

1.1
4節点要素選択的低減積分
8節点要素選択的低減積分
9節点要素選択的低減積分

1.0510.950.90.85

0.70.80.600.10.20.3
ポアソン比

0.40.5101001000
部材長(L)

10000100000

図 12 鉛直荷重作用

図 16 鉛直荷重作用(体積ロッキング問題)
体積ロッキング問題)＜選択的低減積分＞(せん断ロッキング問題)<)<選択的低減積分>
せん断ロッキング問題)<選択的低減積分選択的低減積分>ε xx − ε yy = D2 + xD5 + yD4

2.4 ゼロエネルギーモード
先の結果のように低減積分法、選択的低2ε xy = D 3 + xD 4 + yD 5

減積分法共に、完全積分の時に見られたよ

これらの Di は、要素内の変位式

うな極度の精度低下は見られない。しかし、

u = N 1u1 + N 2 u 2 + N 3 u 3 + N 4 u 4 = C1 + C 2 x + C 3 y + C 4 xy低減積分法においては、ひずみエネルギーv = N1v1 + N 2 v 2 + N 3 v3 + N 4 v 4 = F1 + F2 x + F3 y + F4 xyに含まれない変形モード(ゼロエネルギーから得られる。

モード)が剛性行列に含まれている．このよ低減積分時は 4 節点四辺形要素の積分点

うな剛性行列は、支持条件や荷重条件によは、x=y=0 なので、M4 と M5 のひずみはゼっては解が得られない場合があり，また，ロとなり、ひずみエネルギーに含まれない解が得られても、通常では起こりえない変ことになる。よって、これらのモードはゼ形状態が表れることがある。それぞれの要ロエネルギーモードとなる。8 節点、9 節点

素のゼロエネルギーモードは図 17 のよう四辺形要素においても同様に調べると図になる。この低減積分時のゼロエネルギー

17 のようなゼロエネルギーモードが発生

モードの数と形は，各要素について、ひずしていることが分かる。

みを変形することにより、理論的に導き出すことが出来る。ここでは，4 節点四辺形

表 1 4 節点四辺形要素の要素のゼロエネルギーモードを具体的に示

節点変位に対するひずみ

す．変形モードは図 18 のようになる。また、それぞれのモード毎の節点変位とひずみの関係は、表 1 のようになる．

節点変位ひずみ

モード U1 U2 u3 U4 V1 V2 V3 v4 εxx+εyy εxx-εyy 2εxyM1-1 1 -1 1 -1 -1 1 1400M2-1 1 -1 1 1 1 -1 -1040M3-1 -1 1 1 -1 1 -1 1004M41 -1 -1 1 0 0 0 02y2y2xM50 0 0 0 0 -1 -1 12x-2x2yまた、極端に短い梁に低減積分法を適用
図 17 ゼロエネルギーモード

すると大きな誤差を生む場合がある。4 節点四辺形要素の部材長が L=25.9 の時の変形状態を図 19、また、剛性行列の最小固有値に対応する固有ベクトルを図 20 に示す．

図 18 4 節点四辺形要素の変形モード

要素内のひずみεxx+εyy,εxx-εyy,2εxy は、節点変位から次のように表される．
ε xx + ε yy = D1 + xD5 + yD4

図 19 4 節点要素の変形モード(L=2(L=25.9)
節点要素の変形モード(L=25.9)

図 20 4 節点要素の固有ベクトル(L=29.7)
節点要素の固有ベクトル(L=29.7)

これより，先に述べたゼロエネルギーモ

Reduced Integration and Theードが発生していることは明らかである。

Shear-Flexible Beam Element,

図 21 に 4 節点四辺形要素の部材長を変えた

Int .j .numer, methods engng., Vol.18,ときの剛性行列の固有値を示す。

pp195-210, 1982.

このように、部材長が短いときに，固有(3) 土木学会構造工学委員会

計算力学とベクトルにゼロエネルギーモードを含んだ

その応用に関する研究小委員会：構造固有値が，梁の変形モードを表す固有値よ工学シリーズ 7 構造工学における計算り小さくなり，本来は，発生しないはずの力学の基礎と応用, 土木学会, 1996.

変形が顕在化することが分かる。このような現象は、9 節点四辺形要素においても見られた。

100000

固有値(λ)

100001000
ゼロエネルギーモード発生
標準モード発生1001011010010000.1
部材長(L)

図 21 部材長に対する固有値(4 節点四辺形要素
節点四辺形要素)

3. 結論
本研究では、陽な形での式の変形により，ゼロエネルギーモードの数と形状を導くことが出来た。
また，極めて短い梁において低減積分法を用いると，特異な変形状態が現れ，大きな誤差を生じる場合があるということが分かった。
4. 参考文献(1) N.Bic’anic’ and E.Hinton：Spurious
Modes in Two-Dimensional
Isoparametric Elements, Int. j. numer,
methods engng., Vol14, pp1545-1557, 1979.
(2)G.Prathap and G.R.Bhashyam