M523 初期不整を有する二軸力を受ける球殻の非線形挙動に関する研究

平松宗也（指導教員：正岡，柴原）
９節点要素は Fig. 2 (a)，16 節点要素は Fig. 2 (b)に示す
Tying Point を使用し，歪を内挿した．内挿関数は一般に式（2）に示されるようなラグランジェ補間を用いて導いた．また面内せん断歪γξηに関しては Tying Point が少なくなっており，内挿関数の次数が低減される．

１．緒言
薄板溶接構造物には溶接時の熱影響等により初期不整が発生する．初期不整を有する球殻の座屈強度は古典座屈理論によって求められる座屈強度より極端に低い値を示すことが知られている．例えば LNG 船球殻タンクには二軸力が作用することにより座屈などの非線形挙動が発生し，初期不整の影響による強度低下が生じるが，現在座屈強度低下率は実験から得られた知見による式により求められている．本研究では，より一般的で合理的に初期不整を有する球殻の座屈強度低下率を解析するために非線形有限要素法による解析手法を開発した．

３．ベンチマークテスト
本研究で作成した MITC 曲面シェル要素の精度を確かめるため２種類のベンチマークテストを行った．なお，ベンチマークテストは線形弾性範囲において行った．

２．有限要素法解析
現在研究室で開発されている有限要素法プログラム
ULTSTRUCT は MITC 平面シェル要素を使用している．
本研究では曲面構造物である球殻を解析するため現在の MITC 平面シェル要素に代わる MITC 曲面シェル要素を導入した．MITC 曲面シェル要素は埋め込み座標を用いることにより１要素内の節点を１平面上に存在する必要が無い．これにより要素自体が平面ではなく，曲面をもつ要素であり，平面を表すことはもちろん，曲面を表すにも長けた要素であると考えられる．
また shear locking の一因が面外せん断歪の内挿関数にあることが明らかになり，近年，面外せん断歪の内挿方法を見直す要素が種々開発されている．これらの要素における面外せん断歪は，従来の方法に従って変位から内挿して求めるのではなく，あるサンプリング点の面外せん断歪を再定義した内挿関数から新たに求められる．今回の４節点 MITC 曲面シェル要素では Fig. 1 のように４つのサンプリング点（Tying Point）より面外せん断歪を内挿した．このときの内挿関数は式(1)に示す．
また４節点 MITC 曲面シェル要素に続いて９節点，１６節点の MITC 曲面シェル要素を導入した．多節点の要素を導入することにより内挿関数が多項式になるので，少ない要素で構造物の詳細な変形，応力を表現できると考えられる．多節点要素においては５つの歪成分εξξ，εηη，γξη，γξζ，γηζを Tying Point を使い内挿する．

hhg A xz

g Bhz

xxh1(1 + h )g xzA + 1 (1 - h )g xzC2211D()
g hz = 1 + x g hz + (1 - x )g hzB
22g xz =hg Dhz
g C xz・・・(1)xx

Fig. 1 Sampling points of transverse shear strain of MITC4
a aa aa abba ca c(ⅰ)
b(ⅱ)

b=(ⅲ)

aba(ⅱ)b35b = 0 .861
c = 0 .339
a=a = 1/ 3
aaab(ⅰ)
bbaaa

35(ⅲ)
cc(b)MITC16(a)MITC9Components (ⅰ) e ， g xz (ⅱ) g (ⅲ) e ， g hz
xxxhhhPn ( x ) =M523-1nåk = 0Õ (x - xi )´ ykÕ (xk - xi )・・・(2)

Fig. 2 Strain interpolations and tying points of the MITC9 and MITC16

3.1 Pinched Cylinder
Fig. 3 に示すような円筒のモデルを使用した．円筒の４．球殻の解析
本章ではMITCシェル要素を用いて曲面構造物である中央の上下方向より集中荷重 P=1 を作用させ，荷重点の変位を理論値 1)と比較した．境界条件は円筒の両端においては長さ方向変位を除いて固定とし，両端以外においてはすべて自由とした．節点数は周方向に 48 節点，長さ方向に 49 節点である．また図中の R は半径で 300，L は全長で 600 である．
また板厚は3.0，ヤング率は3.0×106，ポアソン比は 0.3 として解析を行った．結果を Table 1 に示す．Table 1 は厳密解を１とした時の有限要素法解である．

球殻の解析を行う．初期不整の影響により球殻の座屈強度は低下する．また球殻は座屈後の荷重保持能力が低下する．そこで後座屈の挙動を現す b-factor について詳しく考える．b-factor は球殻や円筒等の座屈後の挙動を示す係数であり b-factor がプラスの時は座屈後も荷重保持能力は増加し，b-factor がマイナスの場合は座屈後の荷重保持能力は低下する．円筒，球殻の場合には一般に b-factorはマイナスである．この初期不整による座屈強度の低下率と座屈後の荷重保持能力の低下を精度良く推定することは重要なことである．そこでここでは有限要素法によって解析を行うことにより球殻の座屈強度の低下率を求
め過去の実験的研究（1981）など 3) 4)と比較検討する．

3.2 上部 18°をカットした Hemispherical shell
Fig. 4 に示すような上部 18°をカットした半径 10 の1/8 球モデルを使用してベンチマークテストを行った．
このモデルは双曲面を有し曲面シェル要素のベンチマークテストとして有効である．荷重を点Ａ，点Ｂに P=1 作用させ荷重点の変位を理論値 2)と比較した．モデルの板厚は 0.04，ヤング率は 6.8×107，ポアソン比を 0.3 として解析を行った．節点数は辺 AB 辺 AD 共に 49 節点を持つように分割した．結果を Table 2 に示す．Table 2 は厳密解を１とした時の有限要素法解である．
以上のベンチマークテストをおこなったことにより以下の知見を得た．
１．曲面シェル要素はどの要素も正確な値を得ることができる．
２．曲面構造物では多節点要素を用いると少ない要素分割でも精度の良い解を得ることができる．
P
Ｌ/ 2Ｌ/ 2
Table 1 Accuracy of MITC elements

4.1 球殻の理論
初期不整が強度低下におよぼす影響は以下のKoiterの座屈後の挙動を示す理論式（3）5）を展開することにより推測される．

(1 - r )x + ax + bx +  = rx + order(x ,x ) ・・・（3）2(3)式でρはノックダウンファクターと呼ばれるもので初期たわみを有することでの強度低下を示すものであり式中の b が b-factor である．
（3）式において球殻の場合 aは 0 となり，また高次の項を無視すると（4）式となる．
æl öl3ç1 ÷ç l ÷x + bx = l xèCL øCL

ＭＩＴＣ4
ＭＩＴＣ9ＭＩＴＣ16・・・(4)（4）式を両辺ξで微分してξを求め，（4）式に再度代入することで（5）式が得られる．
（5）式をグラフ化するとFig. 5 ようになる．Fig. 5 より b-factor と初期不整 x が強度低下におよぼす影響が大きいことが分かる．

Exact = 0.18248×10-4R33æl ö2 3 3çç1 ÷÷ =2è lCL øw FEM / w Exact.
0.87
1.03
1.04

( - b )x ll・・・(5)CLl / lCL = r
1.25

1PFig. 3.Pinched cylinder with rigid end diaphragms

0.75

0. 5DC

Table 2 Accuracy of MITC elements

0.25

Exact = 0.0924

APB
P

ＭＩＴＣ4
ＭＩＴＣ9
ＭＩＴＣ16-bx

0w FEM / w Exact.
0.89
1.06
1.04

00.2

0.4

0.6

0.8

1Fig. 5 Effect of initial deflection and b-factor on
knockdown factor

Fig. 4Hemispherical shell with 18°holl

M523-2

4.2 解析対象Odland の実験的および理論的検討 3)に用いられた球殻ρ，後座屈挙動係数ｂと解析解によるρ，ｂとの比較を行う．ここでは Odland，Dnv の Classification Notes と比モデルである Toroidal Shell Segment を解析対象とした．
解析対象を Fig. 6 に示す．詳しい寸法，材料定数は Table
3 に示す．Table3 中における R は半径，t は板厚，νはポアソン比，E はヤング率である．

較する．Dnv の Class-Notes における球殻の座屈に関するルールは以下のように定められている．
ノックダウンファクターを考慮に入れた座屈荷重は以下の(8)式で求められる．
・・・(8)
N CR = Fs 0 t

4.3 荷重条件
荷重は子午線方向の圧縮または引っ張り P，外圧 p の２つを考えた．外圧のモデル化は各節点での法線方向にその点の要素の面積に対して力を作用させる．P と p の合力によって子午線φ方向にσφが作用し，面内θ方向にσθが作用する．この２つの外力を変化させてσφ/σθを変化させる． ここで，σφ，σθは(6)式で表される．
(6)式中のφ0 は Fig. 6 に示す．

ここでΦは塑性修正式，σ0 は降伏応力，t は板厚である．
塑性修正式は(9)式で表される．(9)式中の l は(10)式で表され，降伏応力と弾性状態での座屈応力の比である．

1 pR
P
cos 2 f0 2pR
t 21æ
1 pR
P ös q = - ç pR cos 2 f0 ÷tèt 22pR øsf = -1F=
l=s0
se・・・(10)

ここで(10)式中のσe は(11)式である．
se = r・・・(6)Eh = (1 - r ) 2 3n(q + f 0 ) ・・・(7)a4.5 有限要素法解と解析解の比較
本節では以上の解析条件において非線形有限要素法解析を行い，解析結果より得たノックダウンファクター

・・・(11)

f0Fig. 6 Toroidal Shell Segment model
Table 3 Dimensions and material properties of the model
R[mm]

t[mm]

L[mm]

ν

E[MPa]

1000

2800

0.370000・・・(12)1.15(sf / sq )・・・(13)b = -0.5e

解析的には(13)式で b を求め，(12)式を解くことでρを得る．
また，(8)式に(9)，(10)，(11)式を代入すると(14)式となる．
ここで弾性解析を行う時(14)式中のσ0 は無限大と考えられるので(14)式は(15)式となる．
N =tCR・・・(14)11+ss220e

N =s t = rtL

d3 3-b r2t(12)式は Koiter の後座屈理論の式であり，ここで t は板厚であり b-factor， b は(13)式のように定められている．
またd は初期撓み量である．

CRR

t3(1 -n 2 ) R

ここでノックダウンファクターρは以下の(12)式を h=0として解くことで求められる．

4.4 初期撓み
初期撓みは，1/4 モデルに対して周方向に 10 半波，子午線方向に１または２半波となるように，式(7)で与えた．
ただし，ここで m=10，n=1 または 2 である．
また解析時，初期撓み小という場合は(7)式の 1/100 の初期たわみを与えたものである．最大初期撓み量δ 0 は DNV のClassification Notes, No.30.34)を参考に与えた．
d = d 0 cos (4 m f )sin・・・(9)1 + l4eE

t3(1 - n ) R2t = rs tCL・・・(15)(15)式中のσCL は球殻に一様圧縮が作用した時の古典座屈応力であるので有限要素法解析による b-factor，ρは解析によって得られたσθをσCL で除することでρを求め，ρを(12)式に代入することで b-factor が求まる．
ここで有限要素法解析を行ったσφ /σθ＝－0.296，－0.145，0.096，0.49，1.38 の 5 パターンについてのルール式と有限要素法解析によるρを求め，Table 4 に示す．
Table 4 より，σφ/σθ＝1.38 の最大初期撓みが 0.0067 の時を除いてルール式のノックダウンファクターと有限要M523-3

1.38

また Table 4 で求めたノックダウンファクターより有限要素法解析による b-factor を求め Fig. 7 に示す．縦軸は荷重条件σφ/σθで，横軸は b-factor である．また図中の黒丸破線は(13)式で示される Dnv Classification Notes における定式であり，二重丸破線は Odland による定式である．
Fig. 7 より初期撓みが大きい場合，有限要素法解析による b-factor は Dnv 式，Odland 式と非常に近い値を示すことが分かる．

1

4.6 塑性修正についての考察
4.5では有限要素法解析における弾性解析を行い(14)式における降伏応力を無限大と考えることでルール式との比較を行った．本節では材料的非線形性を考慮した有限要素法解析を行うことでルール式の塑性修正を考慮した解と有限要素法解について比較する．
荷重条件はσφ/σθ＝－0.296 において降伏応力を127[MPa]として解析を行った．その他の解析条件は弾性解析時と同じである．解析解において塑性修正を考慮する時，(10)式のσ0 を 127[MPa]とすることで(15)式は(16)式となる． ここでρ’は初期不整による強度低下率と塑性化による強度低下率の両方を考慮に入れた強度低下率である．Table 5 はρ’をルール式と有限要素法解において比較したものである．
塑性化を考慮に入れた場合もルール式による強度低下率と有限要素法解析による強度低下率は非常に近い値を示す．

N CR = F s 0 t = 0.355s 0t = r 's 0t ・・・(16)σφ/ σθ

素法解析によるノックダウンファクターが良く一致していることが分かる．

1.15(σ /σ )φθDnv– 0.5e
2.19(σφ /σ θ)Odland – 0.225em=10，n=2(Large init.def.)
m=10，n=2(Small init.def.)
m=10，n=1(Large init.def.)
m=10，n=1(Small init.def.)

0–1
0.01 0.02

0.06 0.1

0.2

0.6 1

52–bFig. 7 b-factor for the spherical shell segments
Table5 Knock down factor ρ’ of analytical solution and FEM
σφ/σθδ0-0.300.41

ρ'Dnv

FEM(n=2)

FEM(n=1)0.360.32

0.31

５．結言
本研究においては多節点MITC 曲面シェル要素を導入し，曲面構造物の解析を非常に精度よく解析できることを確認した．さらに MITC シェル要素を用いて初期不整形状，初期不整量，応力状態を変化させて球殻の解析を行い，強度低下率を算出し，過去の実験と良く一致することを確認した．この結果より非線形有限要素法解析による初期不整を有する球殻の座屈強度低下の解析手法の有用性を示すことができた．

Table4 Knock down factor ρ of analytical solution and FEM
σφ/σθ
-0.30
-0.15
0.08690.49
1.38

δ0参考文献ρDnv

FEM(n=2)

FEM(n=1)

0.410.52

0.52

0.55

0.0041

0.97

1

0.99

0.59

0.43

0.43

0.41

0.0059

0.96

0.98

0.98

0.75

0.36

0.37

0.34

0.0075

0.94

0.98

0.98

0.73

0.32

0.31

0.35

0.0073

0.94

0.91

0.92

0.67
0.0067

0.24
0.92

0.24
0.71

0.24
0.72

1）
Ted Belytschko and Bak Leong Wong : Assumed Strain Stabilization Procedure
for the 9-Node Lagrange Shell Elements, International Journal for Numerical
Methods in Enginerring, Vol28, 385-414 (1989)
2）M.L. Buclem and K.J. Bathe: Higher-Order MITC General Shell Elements,
International Journal for Numerical Methods in Engineering, Vol36, 3729-3754(1993)
3）Odland, J. : Theoretical and Experimental Buckling Loads of Imperfect
Spherical Shell Segments, Journal of Ship Research, Vol.25, No. 3, Sept. 1981
4）Dnv : Classification Notes“Buckling Criteria of LNG Spherical Cargo Tank
Containment Systems-Skirt and Sphere” Note No.30.3 December 1997
5）Hutchinson, W.T. : Imperfection Sensitivity of Externally Pressurized Spherical
Shells, J. Apple. Mech. 34, 1967

M523-4