長崎大学工学部研究報告 第27巻 第48号 平成9年1月
一面せん断試験 にお ける粒状体 の進行性破壊 と ダイレイタンシー特性 の個別要素法 による把握
****彦強
由
橋村
棚西
151
*****一郎
正英
崎山
潰木
A St udyonPr ogr e s s i veFai l ur eandDi l at anc yChar ac t e r i s t i c sof Gr anul arMat e r i alunde rDi r e c tShe arUs i ng
byDi s t i nc tEl e me ntMe t hod
by
YoshihikoTANABASHI*,SeiichiHAMASAKI** TsuyoshiNISHIMURA*** andHideoKIYAMA***
Itisfairlyimportanttoestimateafeatureofprogressivefailureandadilatancycharacteristicsof granularmaterialsuchassand,andadevelopmentprocessofshearbandhasbeenrecentlynoteworthy concerningwithprogressivefailureinthemicro‑mechanics. Thispaperfirstdescribesthetestresults ofdirectsheartestforgrassbeadswhichhaveequivalentdiameters. Thepapernextdealswitha numericalsimulationfortheabovetestsusingbyDistinctElementMethod(hereafterreferstoasDEM) developedbyCundall,P.A.(1971)atfirstasacomputermodelforsimulatingprogressivelargescale movementinblockyrocksystems.IthasbeenclarifiedthatDEM modifiedbyoneoftheauthorscan estimatesufficientlythefeatureofprogressivefailureandthedilatancycharacteristicsofgranular materialfrom thecomparisonwithobservedandcalculatedbehaviorofgranularmaterialunderdirect shear.
1.は じ め に
不連続体 を取 り扱 う1手法 として, カ ン ドル(Cun‑
dall,P.A.,1971)の提唱 した個別要素法 (DistinctEle一 mentMethod,以下DEM と略称)が ある.これ は,砂 の ような粒状体や亀裂性岩盤 な どの動的挙動 を取 り扱 うのに適 している. この方法 を用いて,一面せ ん断試 験や平面 ひずみ試験等 をモデル とし,粒状体 の進行性 破壊 をシ ミュレー トす ることは可能である. それ と同 時 に,粒状体特有 のダイ レイタ ンシー特性 (せん断過程 に生 じる体積変化)や,せん断帯形成過程 を観察す るこ
とも容易 にで きる.なかで も,粒状体 を取 り扱 う上 で ダイレイタンシー特性 を把握す ることは必要不可欠で ある.本研究で は, このダイ レイタンシー特性及び進 行性破壊 を,一面せん断試験 をモデル とし,粒子配列 等 の視点か ら検討 を行 ったので, ここに報告す る.
2.個別要素法(DEN) 2. 1 DEM:の概要
DEMで は,岩盤や粒状体 を完全 に切 り離 された ブ ロ ックの集合体 として扱 うので,通常 の有 限 要素 法 平成8年10月25日受理
*社会開発工学科 (CivilEngineeringDepartment)
**大学院修士過程社会開発工学専攻 (GraduateStudent,CivilEngineeringSpecialty)
***鳥取大学工学部 (FacultyofEngineering,TottoriUniversity)
152 棚橋 由彦 ･潰崎正一 ･西村 強 ･木山英郎 (Finite Element Method)や剛体 ばねモデル(Rigid
BodiesSpringModel)で は表現で きない,動的な破壊 過程 の解析 が可能 である. この手法 は,剛体 で表 され る各要素の境界部分 の接触力,及 びそれ によって引 き 起 こされ る要素変位 を,運動方程式の時間差分 によっ て各時間 ご とに追跡 してい くもので,時間依存 の問題 等 に有効 で,物理量 の離散化 によ り大変形問題 に対 し て も適用可能である.
2.2 DEn4の原理1)
2. 2. 1 運動方程式
剛体 ブロックの運動 には,重心 の並進運動 と重心 ま わ りの回転運動 を伴 う. また,2つの粒子が接触 ある い は衝突す る時,粒子 は完全 な弾性体で はない し, ま た接触点近傍 は,局所的 な塑性変形が生 じ,完全弾性 衝突 とはな らない.粒子 の もつ弾性的,非弾性的性質 は,接触点間 に挿入 した フォーク ト(Voigt)モデル (弾 性 スプ リング と粘性 ダ ッシュポ ッ トの並列配列)で表 現す る(Fig.1).
∩エ1yy不し
(a) Normaldirection (b) Tangentialdirection Fig.1 Modelofcontactpoint
(1)重心 の並進運動 に関す る運動方程式
粒子の質量m,並進 に伴 う未知変位 u,弾性 スプ リ ングの剛性定数K,粘性 ダ ッシュポ ッ トの粘性定数甲 とすれば,並進 に伴 う運動方程式 は次の様 になる.
m留 +揺 +Ku‑F (1) ここに,F は外力 を表す.
(2)垂心 まわ りの回転 に関す る運動方程式
粒子の質量 mお よび慣性 モーメ ン トI,回転 に伴 う 未知 回転 量 ￠,弾性 ス プ リングの剛性定 数 且,粘性 ダ ッシュポ ッ トの粘性定数 符とすれ ば,回転 に関す る 運動方程式 は次 の様 になる.
･掌 ･W2普 +Kr24‑M (2) ここに,〟 は外力 によるモーメン トを,γは重心か ら の尿巨離 を表す.
2.2. 2 接触力
2要素i,jの接触面 に作用す る力 を,法線方向に作 用 す る圧縮力 fnと,接線 方 向 に作 用 す るせ ん断力fs
(要素 tに関 して時計 回 りを正)に分 けて考 える.
(1) 法線方向の作用力
微小時間増分△t間の法線方向の相対変位増 分△un に比例 した坑力増分△enを生 じる弾性 スプ リング(剛 性定数Kn)と,相対変位速度△un/△tに比例 した坑力 増分△dnを生 じる粘性 ダ ッシュポ ッ ト(粘性定数 恥)
の並列配列 を考 える(Fig.1(a)).△en,△dnは次式で 表せ る.
Aen‑Kn･Au,I
(3)
△dn‑ qn
･
但 し,圧縮力 を正 とす る.
したが って,時刻 tにおいて接線方向に作用す る弾 性抗力[en]tと粘性抗力[dnL は次式の様 にな る.
[en]t‑[en]t̲At十△en ldn]t‑Adn
また,上式 には次式 の条件が付 され る.
[en]t<0 の とき len]t‑ldn]t‑0
(4)
(5) 以上 よ り,時刻 tにおける2要素間の法線方向圧縮力 lfn]士は次式で計算 され る.
lfn]f‑[en]t+[dn]t (6) (2)接線方向の作用力
接線 方向の相対変位増分△us(要素 iに関 して時計 回 りを正)に対 して も同様 にせ ん 断 抗 力 を与 え る フォー ク ト(Voigt)モデル (せ ん断剛性 定数Ks,せ ん 断粘性定数 恥)を考 える(Fig.1(b)).△es,△dsは次式 で表 せ る.
Aes‑Ks･Aus
(7)
△ds‑qs
･
したが って,時刻 tにお ける接線方向の弾性抗力【es]t と粘性抗力[ds]t(いず れ も要素iに閑 し時計 回 りを 正)は次式 の様 になる.
一面せ ん断試験 にお ける粒状体 の進行性破壊 と ダイレイタンシー特性 の個別要素法 による把握 les]t‑[es]卜At+△es
lds]t‑△ds
上式 には,次の2つの条件が付 され る.
len]t<O の とき[es]t‑[ds]t‑0 lles]tI>FLlen]tの とき
‡ l e s ] l
‑FLlen]tXSING([es]t)[4S]t‑0
(8)
(9.a)
(9.b)
ここに,両 ま粒子間の摩擦係数,SING(Z)は変数Zの 正負 を表す もの とす る. これ らの条件 は,接触点近傍 のせ ん断変形が,主 として要素間の摩擦力 によって生 ず ることを意味 し,条件式前者,式(9.a)は,非接触状 態 を,後者式(9.b)は,摩擦力の限界 をそれぞれ表 して い る.以上 よ り,時刻 tにおける2要素間の接線方向 のせ ん断力[fs]tは次式 で計算 され る.
lfs]t‑les]t+ldslt (10) 2. 2. 3 運動方程式の差分近似
注 目す る要素iと接触す るすべての要素 再 こついて 先述 した様 な形 で,接触 力[fn]t,[fs]tが計算 され る と,要素 iに関す るそれ らのX方向の分力 xi,Y方向 の分力 yiな らびに中心 回 りのモーメ ン トMi(反 時計 回 りを正)は次式で求 め られ る(Fig.2参照).
lxi]‑冒∫(‑[fn]tcosaij+lfs]tsinaid)+mag lyg]‑冒(‑[J fn]tsinαid‑[fs]tcosai,) (ll)
lMi]ニーri･∑(J[fs]t)
ここに,∑ は要素i iに接触す るすべての要素 jに関す る総和 を表 し,またmiは要素 tの質量であ り,migの 項 は重力がX方向 に作用 す るこ とによる.加速度 を 作用力 の陽関数 とみな して変形 した運動方程式 (式省 略)を,時間増 分 △才で差分近似 す れ ば,作 用 力 は式 (ll)で与 えられ るか ら,時刻 tにおける加速度 は,吹 式の様 に求 まる.
Fig.2 Contactforces,fn,fsandc0‑0rdinatex,y
153
(12)
ここに,Ziは要素 iの慣性 モーメン トを表す.
時刻tにおける変位速度 は,式(12)を時間増分△tに 関 して積分 して,
lLii]t‑[L;i]t̲At+[;'Z･]t･△t lL;i]t‑li;i]トAt+[b'i]t･△t
[右]t‑[巌]t̲At+[さi] f･△t
(13)
式(13)をさらに△tで積分すれ ば,時間増分△t間の変 位増分 は次式 となる.
[△ui]t‑([△ui]t̲At+[△L;i]t･△t)/2
[△vi]戸([△vi]t̲At+[△L;i]t･△t)/2 (14) [△Qi]t‑([△動,]卜At+[△4Si]t･△t)/2
式(3)か ら(14)に至 る一連 の過程 で定 め られ た変位 増 分 を(時刻 tか らt十△tまで)新たな変位増分 に仮 定 して,再 び演算 を繰 り返す. この ようにして,時間増 分△t毎 の変位増分が逐次計算で きる.Fig,3に, こ の繰 り返 し演算過程 を示す.
Fig.3 ProcedureofcalculationbyDEM2)
3. 実 験5)
3. 1 実験の概要
本実験 は,通常 の一面せん断試験装置 を用いて,球 形ガラス ビーズ (粒径3.0mm,均等係数 Uc‑1.0)のみ の定圧せん断試験 を行 った.後 に述べ る解析 との比較 を考慮 し,せ ん断容器等 はTablelの様 に設定 した.
試料 は,比較 的緩 く充填 した もの と,密 に充填 した も のの2種 とした.その際,充填の度合 を間隙比e,及び 平均配位 数 で表 す とTable2の様 にな る.配位 数 と は,注 目す る要素 の接 点 数 の こ とを言 う(Fig.4参 照).
154 棚橋 由彦 ･漬崎正一 ･西村 強 ･木 山英郎 Table1 Experimentalcondition
shearboX diameter:6cm bight:3cm*
shearVelocity 1mm/min(0.00167cm/S)
*解析時の1要素への負担 を軽減 させ るため通常 より 1cm高 くした
Table2 Voidratio,C0‑Ordinationnumber
0.879 6.324▲,6.386▲ ▲
0.618 9.129▲,7.421▲▲
▲smith3),▲▲fields)
J':C0‑Ordinationnumber=4
Fig.4 Difinitionofc0‑Ordinationnumber
3. 2 実験結果 と考察
Fig.5に,間隙比 e‑0.879の場合 の応力比 qと相 対変位Dの関係 を,Fig.6に体積 ひずみEvと相対変 位,Dの関係 を示す.ここに応力比 7‑ど/6であ り,T
はせ ん断応力,Jは垂直応力である. また,Fig.7,8 は同様 に間隙比e‑0.618の場合 を示 した ものであ る.
試料 にガラス ビーズ を用いたため, データにかな りの ば らつ きが見 られ る.応力比 甲と相対 変位Dの関係 か ら,間隙比e‑0.879(Fig.5)が延性的傾 向を示 して いるのに対 し,0.618(Fig.7)の場合 は脆性的傾向 を示 していると言 える.一方,体積 ひずみ Evと相対変位D の関係 では,間隙比e‑0.879(Fig.6)初期 にお ける負 のダイレイタンシー,及び間隙比e‑0.618(Fig.8)に おける正 のダイレイタンシー といった,一般的なダイ レイタンシー特性が表れている.またFig.6とFig.8 か ら,粒状体 の体積 ひずみEvが応力比 甲によ り,一義 的 に規定 され るとい う一般 の傾 向が認 め られ る.
0864210000LA40!teJSSaL)S
U=0.28
J=0.0=0.1564
0.0 0.2 0.4 0.6 0.8 1.0
relativedisplacement,D(cm)
Fig.5 0bservedrelationshipbetweenstressratio andrelativedisplacement(e‑0.879)
個別皿皿E;.01皿000000
'2(u!qS岩の∈nIタ
o=0.56 げ=0.28
0.0 0.2 0.4 0.6 0,8 1.0
relati鳩 displacement,D(cm)
Fig.6 0bservedrelationshipbetweenvolumetric strainandrelativedisplacement(e‑0.879)
64208q)4つも0111100000LALoptilSSa4S
(unitofo:kgqcm 2)
O=0.56 0=0.28
0.0 02 0.4 0.6 0.8 1.0
relati鳩 displacement,D(cm)
Fig.7 0bservedrelationshipbetweenstressratio andrelativedisplacement(e‑0.618)
.10戊朋加皿00000'2‑u凛SO.uPuJn一夕
0.0 0.2 0.4 0.6 0.8 1.0
relati鳩 displacement.D(cm)
Fig.8 0bservedrelationship between volumetric strainandrelativedisplacement(e‑0.618)
一面せん断試験 にお ける粒状体 の進行性破壊 と ダイ レイタ ンシー特性 の個別要素法 による把捉 4.解 析5)
4. 1 解析モデルの概要
等粒径 円形要素 (半径1cm)とし,DEMを用 いて定 圧一面せん断試験 をシュ ミレー トした.一面せん断試 験 をシュ ミレー トす る際のせん断速度 は1cm/Sで,せ ん断容器 の大 きさ は底 辺40cm,高 さ20cmの単位 幅 当 りとした.これ は,実際の試験装置の20/3(約6.67倍) で あ り,先 に述 べたせ ん断速度 は この場 合0.15cm/S (粒径3.0mm)となる.その容器 の最下層 に等粒径 円形 要 素 を n個, その上 層 にn‑1個 規 則 的 に配 列((n
‑1)/n配列 と呼称)した.また,この ようにしてで き る配列 (せん断前の粒子間接触 において)は,すべて配 位数4で ある.DEMに用 いた諸定数 をTable3,4に 示す.Table4中の *で示す値 は,せ ん断箱 のダイ レ イタンシーに伴 う周面摩擦 を軽減 させ るための もので ある.解析過程 は,初期配列時 に発生す る粒子浮遊状 態か ら,静止状態 を得た後,所定の上載荷重(0.1,0.2, 0.4kgf/cm2)を載荷 し(圧密過程),せん断す る(せん断 過程).なお,圧密終了 は,次式 の収束条件 に基づ くも の とす る.
tAHl<LMZT‑g･△t2･1.0×10J3 (15)
△H :時間ステ ップ毎 の載荷盤 の変位
△t:時間増分
g :重力加速度(980cm/S2)
(△t‑1.0×104の時LMZT‑9.8×109(cm))
Table3 Materialproperties4)
radius:㍗(cm) 1.0 density:p.(gf/cm.3) 2.65 young'smodules:E(kgf/cm2) 750
Table4 Interactionproperties4)
particletoparticleparticletowall Kn/pg(cm) 3.64×104 7.28×104
qn/pg(cm/S) 1.53×10 3.06×10 Ks/pg(cm) 0.91×104 1.83×104
qs/pg(cm/S) 0.76×10 1.53×10 At(S) 1.00×10‑4
frictionangle,♂ 30○ 45○,5○* frictioncoefficient,〟 0.577 1,0.087*
155 なお,解析 に用いた粒子配列 は11/12,13/14,17/18の
3パ ター ンである. それ ら各配列 の要素数,及び静止 状態 での粒 子 間接触 角αはTable5に示 す とお りで ある. これか らも分か るとお り,配列数 nが小 さ くな Table5 The number of elements and contact
angle
a汀angement Thenumberofelements contactangle,α ll/12 219 59.7560
13‑/14 189 46.9620
るにつれ,水平方向の力の伝達が卓越 した配列 となる.
4. 2 解析結果 と考察
Fig.9に17/18配列 (上載荷重0.1kgf/cm2)の時間ス テ ップ毎 のプロッター出力結果 を示す.せん断が始 ま ると,Fig.10に示す斜線部分 の接触力が減少 し,さら には消滅す る(Fig.9(g)参照). この過程 ですで に,局 所的 な破壊が生 じている と言 える(step4000‑8000,要 素番号71‑89,156‑174:黒色要素).その局所的破壊が 内部 へ進行 しせん断面 を形成す る(step12000一一20000, 要素番号71‑89,156‑174:黒色要素).
Fig.10 Schematic diagram of movemerlt each particles
また,せん断変形 に伴 う体積膨張 (正のダイレイタン シー)は,せ ん断面形成 に関与す る粒子の乗 り越 え量 に 起因す る と言 える.ここに乗 り越 え量 とは,Fig.11に 示す ような粒子i,j間 にせ ん断力が作用 し,iがjを 乗 り越 える ときのある方向の変化量(d)の ことを差 す
Fig.ll Difinitionofdanda
156 棚橋 由彦 ･潰崎正一 ･西村 強 ･木 山英郎 もの とす る.完全 に乗 り越 えた場合 のdをdmaxとす
れ ば,dmarを静止時の粒子間接触角αと粒子半径 rで 表す と,
dma3‑2r(1‑cosα) (16) とな り,αの増加 と伴 にdmarも増加 する.しか し当然
JFJIJJ JfJ llー●一一一1
●
●(a)timestep
=
1428仲)timestep=4000
也)timestep=8000
の ことなが ら,実際 には この dmaxが供試体 (集合体 と して)の最大膨張量 とはな らない.今,仮 に11/12,13/ 14配列 を例 としてdmaxの計算 を行 うと次の よ うにな
る.
dmax(ll/12)‑2･1･(1‑cos(59.7560))‑0.993cm
(d)timestep
=
12000(e)timestep=16000
(i)timestep
=
2000020 【CZd)
] 11 (EGF)
(g)forcevectors(timestel)=1000)
Fig･9 Calculatedfeatureofprogressivefailureunderdirectshear(verticalstress,♂‑0.1kgf/cm2)
一面せん断試験 における粒状体 の進行性破壊 と ダイレイタンシー特性 の個別要素法 による把握 dmar(13/14)‑2･1･(1‑cos(46.9620))‑0.635cm
ところで,解析結果 か ら最大体積 ひず みEvを読 み取 るとTable6の様 にな る.この ように,要素数200前後 であって も,複雑 に入 り組 みあって集合体 としてのダ イレイタンシー を粒子 レベルか ら求 めることは極 めて 困難である.
Table6 Maximum'volumetricstrainandexpan‑
sionamount
arrangement ll/12 13/14 volumetricstrain,Ev(0/.) 7.84 2ー17
Fig.12に13/14配列 (上載荷重♂‑0.1kgf/cm2)の配 位数,せん断応力 rと相対変位 β の関係 を示す.相対 夏位O.0cmにおける配位数3ゐ要素 は,境界 に隣接 し ているた めである.せ ん断開始 と伴 に,配位数4の要 素数が減少 し,配位数3の要素数が増加 してい る.配 位数4の要素数 と配位数 3の要素数が逆転す る ところ でほぼ ピーク(せん断応力r)を示す.この ことは,先 に 述べた局所的破壊 を思 い起 こす と容易 に理解 で きる.
つ まり,せん断が進む につれFig.13に示す斜線部分 の接触力が減衰,消滅す る. そのた め黒 く塗 りつぶ し た要素が,しだいに配位数4か ら3へ変化す る.13/14 配列 の相対変位1.2‑1.5cm付近 の急激な減少で は,配
N■quaEalaJ0)aqLunuaLA 仰000022つi11 008(ZuPghq)1gSSa.aSJeausS000O(U0
0.0 0.5 1.0 1.5 2.0 2.5 relativedisplacement.D(col)
Fig.12 Calculatedrelationshipsamongthenum‑
berofelementsbyc0‑Ordinationnumber, relativedisplacementandstressratio
せん断力
Fig.13 Featurec0‑Ordinationnumberunderdirect shear
157 位数4,配位数3ともに減 少 (それ に付 随 して配位数
2の増加)している.相対変位1.4cm付近では,配位数 2の要素数が配位数 4よ り多 くなってい る. この こと は, この付近 において,粒子の再配列が行われてい る もの と言 える.
5.実験 と解析の比較
実験値 と解析 を比較す るにあた り, まず初 めに試験 装置等の理 由か らひずみ条件 を異 にす ることを加味 し て もらいたい.前述 とあ る意味で同様 な ことであるが, 配位数 (実験 の場合 は確率変数的分布 のた め平均配位 数)も異 にす る.したが って,定量的な評価 は困難であ り,定性的な評価 に とどまる.Fig.14,15(解析結果) に比較的密 な配列 の11/12と緩 い配列 の13/14の応力比 Uと相対変位Dの関係 を示す.実験結果 に見 られ る, 密 な場合の脆性的傾 向及び緩 い場合の延性的傾 向が, 解析結果 に もあ らわれてい る. また,相対変位初期 の 応力比増分 と相対変位増分 の比が,緩 い状態の ものが 大 きい ことも実験結果 と同様 である.
Fig.16,17は,比較 的密な配列の17/18と緩 い配列
42■lr一l 08641000LLEOptLFSSaJtS
o:0.1 (unitofo:kgqcm 2)
0.0 0.5 1.0 1.5 2.0 2.5
relati鳩 displacement,D(cm)
Fig.14 Calculated relationship between stress ratio and relative displacement (dense arr∴ ll/12)
765432.10000人U00LLJopES切巴‑S
0.0 0.5 1.0 1.5 2.0 2.5
relatiledisplacement,D(cm)
Fig.15 Caicuユated relationship between stress ratioandrelativedisplacement(loosearr∴ 13/14)
