計算工学講演会論文集

Vol.9 (2004年5月)

計算工学会

ASI-Gauss有限要素法における破断･接触アルゴリズムの検討
On Member-Fracture and Contact Algorithms for ASI-Gauss Finite Element Method○ 磯部大吾郎1)， チョウ ミョウ リン2)
Daigoro Isobe and Kyaw Myo Lynn1) 工博

筑波大学助教授 機能工学系（〒305-8573 茨城県つくば市天王台1-1-1, isobe@kz.tsukuba.ac.jp）
2)筑波大学大学院生（〒305-8573 茨城県つくば市天王台1-1-1, kyaw@kz.tsukuba.ac.jp）

Recently, a new finite element code that can be efficiently applied to structural collapse analyses of framed structures underimpact loads has been developed. The code is developed by using the ASI-Gauss technique, a modified version of the formerly
developed Adaptively Shifted Integration (ASI) technique for the linear Timoshenko beam element, which computes highly
accurate elasto-plastic solutions even with the minimum number of elements per member. The ASI-Gauss technique gains still
higher accuracy especially in elastic range, by placing the numerical integration points of the two consecutive elements forming anelastically deformed member in such a way that stresses and strains are evaluated at the Gaussian integration points of thetwo-element member. Moreover, the technique can be used to express member fracture, by shifting the numerical integration point
and by releasing the resultant forces in the element simultaneously. An elemental contact algorithm, which uses geometric relation
for determining contact and gap elements for simulating impact phenomena, is also implemented to the code. Practical results canbe obtained in impact collapse analyses, however, further estimations should be made to the member-fracture and contact
algorithms to improve the accuracy as well as reality in collapse stage. In this paper, some efforts are taken to improve the validityof the algorithms, and it is verified by carrying out simple numerical tests.
Key Words: ASI-Gauss Technique, ASI Technique, Member Fracture, Contact, Impact Loads, FEM

1．はじめに．はじめに
先の世界貿易センタービル崩壊を受け，飛行物体との衝突過程をシミュレートすることを目的とし，衝突荷重下の骨組構造体に対して有効な有限要素解析手法が開発された[1]．本手法には，計算コストを低く抑えることが可能な ASI 法（順応型 Shifted Integration 法）[2]をさらに改良した，ASI-Gauss 法を適用した．ASI-Gauss 法では，２つの要素をサブセット要素として考え，そのガウス積分点に相当する位置に応力評価点を配するように数値積分点をシフトすることで，弾性変位解の精度を向上させている．また，積分点のシフトと同時に断面力を解放することで破断を表現し，幾何学的な位置関係に基づいて要素間をギャップ要素で拘束することで接触を表現可能とした．この方法は，衝突荷重を節点力として加える方法に比べ，衝突物体のパラメータの差異を考慮でき，崩壊過程の解明に有効であると考えられる（図１参照）．
本稿では，これらの破断および接触アルゴリズムに焦点を当て，簡単な 2 部材モデルに対する衝突解析を通じ，要素間の拘束時間を設定するためのパラメータなどについて検討した結果について報告する．

２．ASI-Gauss法
法
２．
線形チモシェンコはり要素を使用する ASI 法では，弾性域での数値積分点位置は線形解析に対する最適位置とする中央点に置き，全塑性断面の発生直後に，その点に正確に塑性ヒンジが形成されるように数値積分点をシフトする[2]．ASI-Gauss 法での変更点は，この弾性域での数値積分点位置のみである．
線形チモシェンコはり要素と，応力評価点あるいは塑性ヒンジ形成位置が陽に与えられている物理モデルである剛体・ばねモデルのひずみエネルギ近似式の等価性を考慮することにより，数値積分点位置と塑性ヒンジ発生Fig. 1 : An example of impact collapse analysis using
member-fracture and contact algorithms
点（または応力評価点）位置の関係は
s1 = Är1 or r1 = Äs1(1)と表現される[3]．ここに，s1 および r1 はそれぞれ，数値積分点位置および塑性ヒンジの位置である．要素全体が弾性域にある場合は，ASI 法ではその最適な積分点位置として要素の中央点を使用するが，ASI-Gauss 法では図 2に示すように 2 つの要素を 1 つのサブセット要素としてまとめ， 3 次はり要素と同様にガウス積分点に相当する位置に応力評価点を配置する．すなわち，2 つの線形チモシェンコはり要素を用い，2 点積分法に基づいた１つの要素を擬似的に構成するのである．これは，ガウス積分点

-11Member0-1/√3rg 0 sg-1 -11/√3sg 0 rg

L(s + sg )2LÅ My (s) = Å My (sg ) Ä Å Vx (sg ) (s + sg )2

1

Stress evaluation pointfy ëFig. 2 : Locations of numerical integration and stressevaluation points in ASI-Gauss techniqueで評価される曲げ変形精度が数学的に保証され，その点が 2 点積分法における最適な積分点位置であることを利用している．この時の数値積分点位置 sg と応力評価点位
置 rg との関係は，(1)式より(2)sg = Ä rg or rg = Ä sgとなる．よって，ガウス積分点に相当する応力評価点位置，およびその位置を表現するための数値積分点位置は，
図 2 に示す通りとなる．この位置が，ASI-Gauss 法の弾性域における数値積分点位置となる．このときの要素剛性マトリックス，一般化ひずみベクトルと断面力増分ベク
トルは次式で表される．
T(3)(4)(5)[K] = L[B(sg )] [De (rg )][B(sg )]fÅ è(rg )g = [B(sg )]fÅ ugfÅ õ(rg )g = [De (rg )]fÅ è(rg )g

ここで L は要素長，[B] はひずみ－変位マトリックス，[De ]は弾性の応力－ひずみマトリックス， fÅ ug は変位増分ベクトルである．
はり理論に従えば，曲げモ－メント Mx ， My とせん断
力 Vx ， Vy との関係は
Vx = ÄdMy;dz

Vy = ÄdMxdz(6)と与えられる．よって，曲げモ－メント増分 Å Mx (s) および Å My (s) の要素長方向分布は， s = sg における曲げモ－メ ン ト 増 分 Å Mx (sg ) ， Å My (sg ) お よ び せ ん 断 力 増 分
Å Vx (sg ) ， Å Vy (sg ) を用いて，次式により近似することができる．

êMx
Mx0

ë2+êMyMy0ë2êNN0

ë2+

êMz
Mz 0ë2T[K] = L[B(s1 )] [Dp (r1 )][B(s1 )]fÅ è(r1 )g = [B(s1 )]fÅ ugfÅ õ(r1 )g = [Dp (r1 )]fÅ è(r1 )g=1(8)(9)(10)(11)

ここで，[Dp ] は塑性の応力－ひずみマトリックスである．
上記のシフト操作を材料性状に合わせて順応的に行うことにより，最小限の要素数で解を求める．メモリ消費量を抑えるためにソルバーには共役傾斜法（CG 法）を採用した．また，動的な問題では振動周期の誤差低減のために分布質量マトリックスを採用し，時間積分には陰解法を用いた．

３．部材破断および要素接触アルゴリズム
３．部材破断および要素接触アルゴリズム
本研究で使用する部材破断アルゴリズムでは，2 要素で構成される部材の曲率と軸ひずみの値によって部材の破断を判定する．線形チモシェンコはり要素は，数値積分点を１つだけ持ち，横たわみの不連続性を許容する剛体・ばねモデルと等価である．そこで，破断面をどちら

s=0
Plastic hinge

s=0

Node+

ここで， Mx ， My ， N ， Mz はそれぞれ x, y 軸回りの曲げモーメント，軸力，ねじりモーメントである．右下添え字の”0”は全断面塑性値であることを示す．降伏関数では，せん断力の影響を無視した．
要素両端のどちらか一端が塑性化した直後の増分ステ
ップにおいては，要素内の数値積分点を(1)式に従ってシフトする．例えば，全塑性断面が要素の右端( s = 1 )に発生したならば，数値積分点は要素の左端( s1 = Ä 1 ，よって
r1 = 1 ) にシフトされ，全塑性断面が要素左端に発生したならば，積分点は要素右端にシフトされる．このときの要素剛性マトリックス，一般化ひずみベクトルと断面力増分ベクトルは，ASI 法の場合[2]と同様に次式で表される．

s=1s=sg

s=sg

(7b)

上式は，曲げモ－メントが要素内で線形に変化し，2 つの端点 (s = Ü 1) のどちらかで，最大絶対値をとることを示している．曲げモ－メントを除く諸量は要素内で定数値をとるので，要素両端の断面力を上式等で算定し，それらの値を降伏関数に代入することにより塑性化を判定す
る．本研究では，以下の降伏関数を用いた．

Element 2
Element 1
sg = − rg = 1- (2/√3)
Numerical integration point(7a)Å Mx (s) = Å Mx (sg ) Ä Å Vy (sg )1Fractures=0Numerical integration point(b) Plastic stage(a) Elastic stage(c) Member fracture
Fig. 3 : Locations of numerical integration points during each stage

A1

Bi1(xi1,yi1,zi1)および Bi2(xi2,yi2,zi2)とおくと，これらの 4 つの節点が同一平面上に存在する条件は，次式のように表される．

Bi1

Colliding element Ecf (x; y; z) ëElement Ei
A2f(yi1 Ä yf 2 )(zi2 Ä zf 2 ) Ä (yi2 Ä yf 2 )(zi1 Ä zf 2 )g(xf 1 Ä xf 2 )+f(xi2 Ä xf 2 )(zi1 Ä zf 2 ) Ä (xi1 Ä xf 2 )(zi2 Ä zf 2 )g(yf 1 Ä yf 2 )Bi2+f(xi1 Ä xf 2 )(yi2 Ä yf 2 ) Ä (xi2 Ä xf 2 )(yi1 Ä yf 2 )g(zf 1 Ä zf 2 ) = 0Fig. 4 : Geometric relation in the determination ofelemental contact(13)

接近する要素どうしが上式を満たし，かつ特定の距離内に存在する条件式，すなわち

Bi1
Colliding element Ec

jA1 Bi1 j + jA1 Bi2 j + jA2 Bi1 j + jA2 Bi2 j î Cl (Lc + Li )A1
A2

Gap elements

Element EiBi2Fig. 5 : Four gap elements binding the two elements
in contact
か一端に設定し断面力を解放すれば，破断の表現が可能となる．このことを利用し，図３に示すように，まず要素の降伏が判定されたら，前節で示したアルゴリズムに従って要素内の数値積分点をシフトし，塑性ヒンジを要素端に発生させる．次に，その要素が破断臨界値を超えた場合に断面力を瞬時に解放し，それまで作用していた内力を解放力として作用させる．いずれの段階でも，同一部材内の隣接要素の積分点位置は中央点にシフトされる．また，破断した要素もその積分点位置は中央点となる．これは，各々の場合で１点積分が最適であるためである．本研究では，破断の判定パラメータとして部材の曲率 îx ， îy および軸ひずみ èz を使用し，これらのうち一つでも以下に示す条件式を満たした場合に，その要素が破断すると判定した．

êîx
îf xë

Ä 1 ï 0 or

ê

îy
îf yë

Ä 1 ï 0 orêèzèf zë

Ä1ï 0(12)ただし，îf x ，îf y は x，y 軸回りの破断曲率，èf z は破断軸ひずみである．部材が破断する際には破断面に新しい節点が設けられ，その２節点で要素質量が等分割されるように設定した．
次に，本研究で使用した要素接触アルゴリズムについて述べる．部材破断のみを考慮した解析では，破断した部材が他の部材をすり抜ける等の挙動を示すため，衝突エネルギが構造物に適切に伝達されない可能性がある．
また，衝突時の力を節点力のみで入力すると，構造物に対する局所的な作用となってしまい，大規模な構造物間の衝突現象を表現するには無理が生じる．そこで，要素間の空間内位置関係により接触を判定する，以下のような要素接触アルゴリズムを採用した．
要素間の接触は，(i)接近する要素どうしが特定の距離内に在ること，(ii)それらの 4 つの節点が同一平面上，またはそれに近い状態に在ること，の 2 つの条件によって判定する．図 4 に，衝突要素と被衝突要素の節点位置関係を示す．衝突要素 Ec の節点をそれぞれ A1(xc1,yc1,zc1)およ び A2(xc2,yc2,zc2) ， 他 の 要 素 Ei(i=1,2,3, … ,n) の 節 点 を(14)を満たした場合には接触と判定する．ここで，Lc およびLi(i=1,2,3,…,n)はそれぞれ衝突要素 Ec の要素長，被衝突要素 Ei(i=1,2,3,…,n)の要素長である．また Cl は接触距離に関する定数で，本来，2 要素の接触時の位置関係により変化するものだが，その平均的な値として Cl=1.8 と設定した．
(14)式を満たすが同一平面上に存在しない要素については，2 つの要素の 4 節点が平面に近い形状を成す条件，すなわち(15)f (x; y; z) î Cfを使用し，この条件を満たした場合に接触と判定した．
ここで，Cf は平面形状の度合いを示す正の定数で，部材の厚み等を考慮して本研究では Cf=3.0×10-4 m3 と設定した．接触を判定された 2 つの要素については，図 5 に示すように節点間に計 4 本のギャップ要素を結合し，要素間を拘束する．ギャップ要素には，他の要素と同一の材料定数と幾何形状パラメータを使用した．

４．数値例
破断・接触アルゴリズムの有効性を検証するため，図 6に示すような簡単な 2 部材モデルの衝突解析を行った．
まず，自由運動をする部材に初速度 50m/s を与え，両端を固定した部材に水平に衝突させた．ここでは，部材が弾性挙動をすると仮定し，ギャップ要素の拘束を外す判定条件について検討した．図 7 に拘束を外さない場合の接触部材中央の節点における変形量の絶対値，およびギャップ要素の曲げひずみエネルギの時刻歴を示す．変
形が戻った点（図中×印）でギャップ要素を削除し拘束を解放すると，図 8 に示すように 2 部材の運動エネルギ(KE)，位置エネルギ(PE)，ひずみエネルギ(SE)の総和は推移する．図から，接触中，拘束解放後においても総エネ

L = 0.5 mLLLL−

I x = I y = 5.21×10 7 m 4
ρ = 7900 kg / m3 , E = 206 GPaυ = 0.3, σ y = 245 MPa

LZL−A = 2.50 ×10 3 m 2

M x 0 ≡ M y 0 = 7.66 ×103 Nm

YN 0 = 6.13 ×105 NX

M z 0 = 1.12 ×105 Nm

Fig. 6 : Two-beam model and the propertiesAbsolute value of deformation
Bending strain energy
Release point

5500.44400.3

30.2

20.1

1

10000
1015
Time (ms)Energy (×103J)

Energy (×103J)

Deformation value (m)0.5ZZ20ZYYXX

X

12 ms

14 ms

16 ms

X

Before impact On impact (10 ms)

10Time (ms)Fig. 8 : Transition of energies in two-beam modelZYYX200Z

Z3020Fig. 7 : Time histories used for determination ofrelease pointY

KE+PE+SE
KE+PE
SE

YX

20 ms

Fig. 9 : Elastic behavior of the two-beam model before and after impact
ルギが保存されることが分かる．変形量の絶対値に基づいて拘束の解放を自動的に判定すると，この場合には 2部材モデルは図 9 のような挙動を示した．
次に，部材の降伏と破断を考慮した衝突解析を行った．
この場合には変形量が戻らず，その値による拘束解放の判定は困難である．また，実際には部材が固まりとなって動き，周囲に与える力学的影響が少ないことが考えられるため，被接触要素が降伏または破断した場合には，拘束の解放は行わないことにした．例として，部材の初速度を 100m/s，破断曲率を îf x = îf y = 0.04, 破断軸ひずみを èf z = 0.3 とした場合の解析結果を図 10 に示す．

５．結論
本研究では，ASI-Gauss 有限要素法における部材破断および要素接触アルゴリズム対し，実際の現象との適合性を高めるための検討を行った．その結果，接触した要素間の拘束時間を変形量で判定すると，エネルギを保存した挙動が得られることが分かった．逆に言えば，エネルギの消失を考慮したい場合にもこの方法は有効であると[2] Y. Toi and D. Isobe. Adaptively Shifted Integration
Technique for Finite Element Collapse Analysis ofFramed Structures. Int. J. Numer. Methods Eng., 36,
1993, 2323-2339.
[3] 都井 裕：骨組構造および回転対称シェル構造の有限要素解析における Shifted Integration 法について,
日本造船学会論文集, 168, 1990, 357-369.

X

Before impact On impact (10 ms)

Z

ZYYX

参考文献
[1] K.M. Lynn and D.Isobe: Structural Collapse
Analysis of Framed Structures under Impact Loads
Using ASI-Gauss Finite Element Method, Proc. 1st Int.
Conf. Design and Analysis of Protective Structures
against Impact/Impulsive/Shock Loads, Tokyo, 2003,
371-387.

Z

Z

Z
Y

考えられる．今後，減衰効果等も導入し，大規模骨組構造に対する解析を実施する予定である．

X

12 ms

Y

14 ms

Z

Y
X

X

16 ms

Y
X

18 msFig. 10 : Behavior of the two-beam model considering yielding and member fracture