数理解析研究所講究録 1472 巻 2006 年 71-80

71

砂山形成過程における頂点移動のダイナミクス II

京都大学・人聞・環境学研究科

千由 (Chiyori Urabe)

占部Graduate School of Human and Environmental,
Kyoto Univiversity
概要

砂山形成過程のダイナミクスについてシミュレーションを用いて研究を行なった。
離散要素法を用い, 砂山上空の一点から 1 粒子つつ時間間隔 $T$ で自由落下させ供給し続けた。 頂点位置の時系列のパワースペクトルを計算すると、パワースペクトルは $T$ によらず広い周波数領域でベキ的になるが、 巾の指数は $T$ に依存する。 $T$ が短いときにはべ $*$ の指数は一 1 に近付き, $T$ がある値を越えて大きくなると $T$ の広い範囲で指数が近似的に一 1.4 になるの領域が現れた。 2 次元の砂山では $T$ をさらに大きくした場合に指数がさらに小さくなる。 $T$ が小さくパワースペクトルが $1/f$ 的になる場合に着目し、 運動エネルギーの時聞変化を用いて雪崩を測定するした。 雪崩の起きる斜颪は時問によって入れ替わり、 その入れ替わりを表す 2 値化した時系列のパワースペクトルは頂点の時系列のパワースペクトルと等しい指数をもつ巾関数で近似されることがわかった。 3 次元の場合には、全粒子の運動量の床と平行な2 成分それぞれの平均を極座標表示した時の角度成分 を雪崩の方向を表す変数とした。 の時系列のパワースペクトルも 2 次元の場合と同様に頂点位置の時系列のパワースペクトルと近似的に等しいことがわかった。 2 値化された時系列の待ち時間分布が指数 $D$ のべ\not\simeq 関数であり, 時系列のパワースペクトルもベキ関数で表されると
B. Lowen と .
き, パワースペクトルの指数は $-3-D$ となることが 1993 年に$1/f$の場合には待ちパワースペクトルがC. Teich により解析的に求められている。
$K$の時系列の待ち時聞分布を時間分布の巾の指数が -9. になると考えられ、 実際に測定すると近似的に一 2 となることから理論により待ち時間分布とパワースペクトル$K(t)$$\theta$$\theta$$\mathrm{S}$,$\mathrm{M}$の関係が説明されることがわかった。

1

はじめに砂山形成過程は雪崩と粒子の積み上げの相互作用からなり、 砂山表面の形状は時々刻々変化し続ける。 本研究では砂山の形成過程を特徴づけるものとして頂点の位置に着目した。 砂山の概形は頂点とすそのを通る円錐形であり、 砂山に供給された粒子は砂山表面に沿って転がるため供給された粒子の堆積するおおよその方向が頂点の位置の情報から予測することができる。 また雪崩によって頂点を含む斜面が劇れれば、 頂点の位置は雪崩と逆方向に移動する。 頂点の位置は形状変化と雪崩の両方とそれぞれ常に密接な関係があるため、 頂点移動のダイナミクスを考えることは砂山形成過程の理解に貢献できるものがあるのではないかと期待される。
本研究では数値シミュレーションを用いて 2 次元の砂山と 3 次元の砂山の形成過程の頂点移動について調べた。 砂山形成過程を維持するために、 時間間隔 $T$ で一粒ずつ砂山に72粒子を供給し続けた。 $T$ が小さいときは雪崩が続いている途申で粒子を供給するため雪崩の起きている砂山斜面は流動状態を保ちやすく、 逆に $T$ が大きいときは雪崩がおわった後で粒子を供給することになるため新たな雪崩が発生しにくく雪崩の発生頻度は低い。
頂点の時系列のパワースペクトルを計算すると、 の大きさによらず巾的になるが、 巾の指数は の大きさに依存することがわかった。 が小さいときには頂点のパワースペク
$T$

$T$

$T$

トルは $1/f$ 的になり特徴的時間が存在しない。 $T$ がやや大きい場合には $T$ の広い範囲で岡じ指数が現れ、 をさらに大きくした場合には指数が更に小さくなる傾向が見られた。
頂点のパワースペクトルが $1/f$ 的になる場合に着目して考えると、 常にどこかで雪崩が起きており、 その雪崩によって頂点は頻繁に動かされる。 雪崩を測るために 2 次元で$T$は砂山を左右半分に分け、 それぞれの領域にある粒子の運動エネルギーの総和を計算し、運動エネルギーの大きさをその瞬間での雪崩の大きさとした。 左右の領域の運動エネルの片方が大きくもう一方はそれに比べて非ギー ,の時系列を見ると、 あるいはの大小は時間によって変化する。 左斜面で雪崩が起き常に小さくなる傾向があり、を右モードと呼び、 左モードのているとき伽 $>k_{T}$ ) を左モード、 右斜面の場合の時系列から 2 値化した時系
蒔を一右モードの時を +1 として運動エネルギー と列 $K$ を作り、 $K$ のパワースペクトルを計算すると、 頂点のパワースペクトルと弼様の指数一 1 の巾関数となることがわかった。 $K$ の時系列と頂点の蒔系列が相関があり、 パワースペクトルが由的であり巾の摺数が 致することより、 頂点移動がモードの switchirig で鍛き換えることができた。 3 次元の砂山の場合には雪崩を測定するために粒子の運動量の平均を調べ、 砂山底面と平行な 2 成分を極座標表示したときの角度成分を調べた結果$h_{l}^{\rho}$$k_{l}$$k_{r}$$k_{r}^{\wedge}$$k_{f)}^{\wedge}k_{r}$$(k_{f}^{\wedge}>k_{l})$$k_{f}^{\wedge}$$1_{\text{、}}$$k_{r}$$-\cdot-$

2 次元の場合と同様に頂点の時系列との問でパワースペクトルの巾の指数が一致することがわかった。
$K$のような 2 値化された時系列のパワースペクトルについては既存の研究がある。 ある

2 値化された時系列の同じ値が継続する蒔問 (待ち時間) がそれぞれ独立であり、 待ち疇
問分布とパワースペクトルが市的になる場合には、 待ち時間分布の巾の指数 $D$ とパワースペクトルの巾の指数 の聞に $\alpha+D=-3$ $(-3<D<-1)$ という関係が解析的に導き出される
実際にシミュレーションにより $K$ の待ち時間分布を計算すると車的になり巾の指数は近似的に -2 となったことから、 モードの切り替わりの待ち時間が指数-2 の巾関数で表されることが事様で時系列のパワースペクトルが $1/f$ 的になることがわ$\alpha$$[1]_{\text{。}}$

かった。

2設定

粒子の運動をシミュレーションするために離散要素法 (DEM) を屠いた [2]。粒子は粘着性がなく、 薗 (体積) 密度は 定の円盤または球であるとする。 接触中の粒子問には弾性力、 粘性力、 クーロン摩擦が、 接触していない場合は重力のみが働くとする。 弾性力を表すために粒子間に線形バネがあると仮定した。最大粒径と最大粒子質量を と $m$ とし、粒径分布は $(0.8d, d)$ で一様分布させる。 2 次元の場合はクーロン摩擦係数とし、$-\cdotarrow$$d$$\mu=\mathrm{t}1.5$

73

図 1: シミュレーションによる砂山。 床平面と平行には床の中心とする。
から粒子を供給する。

$x_{top}$は頂点位置の$x$軸、 鉛直方向に 軸をとり、 原点成分とする。 軸上の頂点からの高さ $H$ の位置
$x$

$y$

$y$

接線方向・法線方向のバネ係数はそれぞれ $k_{t}=2.0\cross 10^{3}[g/d],k_{n}=1.0\rangle(10^{4}[g/d]$ とし、接線方向・法線方向の粘性係数は $\eta_{n}=1.0><10^{2}[\sqrt{g/d}],\eta_{t}=0.0[\sqrt{g/d}]$ とし、 3 次元の場.
合は $\mu=0.2,$$h_{t}^{\wedge}=2.5\cross$$10^{2}[g’/d],k_{n}=\mathrm{L}\mathrm{O}$ $\mathrm{x}10^{4}[g/d],$$\eta_{\tau\iota}=\sqrt{2\pi k_{n}^{\wedge}/3}[\sqrt{g/d}],\eta_{t}=\eta_{n}/2.0$とする。 これらの値を適用するとどちらの場合も跳ね返り係数は約 02 となる。 シミュ
レーションを行う上での時間刻 $At=1.0\rangle\langle$ $10^{-3}[\sqrt{d}/g]$ とする。

頂点の位置は砂山を構成する粒子の中で粒子の重心がもっとも高い位置にあるものとする 図
2 次元の場合は砂山を支える床は粒径 の粒子 80 個を水平に並べたものと$($$d$$1)_{\text{。}}$

軸は床に沿ってとり 軸を鉛直方向上向きとし原点を床の申心とする。
3 次元の場合は床は直径 $30d$ の固定された円で、 床平面を $x-y$ 平面とし鉛直方向上向きに 軸をとり、 原点は 2 次元の場合と同様に床の中心とする。
砂山形成過程における頂点のゆらぎについて調べるため、 粒子を時間聞隔 で砂山に供給し続ける。 供給される粒子は床の中心の真上から自由落下させる。 2 次元の場合は図し、 床は固定し

$x$

$y$

$z$

$T$

1 のように常に頂点からの高さ $H$ で粒子を放つが、 3 次元の場合は実験でも同じ設定で行いやすいように $H$ は床からの高さとする。 また、 床を覆うようになるまで成長させた砂山を初期状態として使う。 簡単のため 3 次元の砂山では内部の粒子を固定した。 床は有限で、 床からこぼれ落ちた粒子は計算から除外するため、 砂山の大きさはほとんど変化しない。

74
初期の砂山の上に粒子を供給し続け、 頂点の時系列が定常状態に落ち着いたあとの時系列を取り出し、 そのパワースペクトル $S(f)$ を計算する。 パワースペクトルを欝算する
で $M$ 等分し、 $m$ 番目の蒔系列のパワースペクトルために、 時系列を時間間隔$T^{(s)}$$S_{\tau r\iota}(f_{j})$を(1)$\mathrm{S}_{m}^{\mathrm{v}}(f_{j})\equiv\frac{1}{N}|\sum_{n=1}^{N}x_{top,n}^{(n\iota)}e^{-\frac{2\pi\dot{x}n}{N}j}|^{2}$と定義する。 但し、$x_{t_{op},7\iota}^{(rr\mathrm{z})}=x_{top}((m+ \frac{\tau\iota}{\mathit{1}^{\mathrm{t}}\mathrm{v}^{\tau}}).T^{(s)}),f_{j}=\frac{j}{T^{(s)}}\circ S_{7\Gamma l}$$S(f_{j})= \frac{1}{\Lambda I}\sum_{m=0}^{\Lambda f-1}S_{m}(f_{j})$の平均.

(2)をパワースペクトル $S(f$ } とする。

32 次元の砂山における結果を除いて、 $x=0$ を中心として左右 2 つの部分に別け、 それぞれの部分に含まれる粒子の運動エネルギーの総和 $K_{l}(t),$ $K_{r}(t)$ を計算した。

$x\leq d$

75
シミュレーションの観測から雪崩は砂山上の 2 つの斜面で同時に起きることは希であり、雪崩が砂山の右 $(x>0)$ の斜面あるいは左の斜面 $(x<0)$ で起きている状態をそれぞれ右モード, 左モードと呼ぶことにする。 2 つのモードの切り替わりを調べるために、 以下のように 2 値化した時系列を定義する。
$K(t)=\{$+1$I\mathrm{f}_{l}(t)<K_{r}(\ell)$-1\emptyset & き(3)

それ以外の符号は時刻で主にどちらの斜面で雪崩が起きているかを示す。 .
$2\sqrt{d}/g$ ,$H=20d$ のときの $K(t)$ との時系列のパワースペクトルを図 4 に示した。
$K(t)$ のパワースペクトルは長いタイムスケーのパワースペクト
ルで $1/f$ 的であり、ルと巾の指数が 致する。また、任意の $x_{top}(t)$に対する $K(t)=-1$ の条件付確率を調べるが増加するに従い増と、 その確率は$K(t)$$t$$T$$=$$x_{t\sigma p}(t)$$x_{b\sigma}(t\}t)$$–\vee\wedge$$x_{t\sigma_{l^{J}}}(t)$す。 このことから それぞれのモードにおいは主に雪崩が起きている斜面とは反対方向に動く傾向があることがわかる。 つまのゆらぎは 2 つのモード聞の切り替わりりとほぼ 致している。て$x_{lo_{\mathit{1}}},(t)$図 3:. の $T$ に対する依存性。

$(\mathrm{a})_{\Omega}$$(\mathrm{t}))\alpha$の $H$ に対する依存性。 エラーバーの長さは標準偏差。

$x_{t_{CJ}r)}$$-\wedge$

4

3 次元の砂山における結果, ztop) とし. テーブ
頂点の座標を (ル表面と平行な成分は (y 一とする。形$x_{t\sigma p},$$y_{8\sigma p}$$x_{\ell op1}$

成過程における” 雪崩と逆方向への頂点移動”
も”’ 積み上がりによる頂点が移動’ もどちらの場合でも水平方向の動きで考えられることか$\circ$ら$(x_{t\mathrm{o}p}, y_{\mathrm{f}o\mathrm{p}})$に着目する。 テーブルの形が円形であり粒子を供給する位置はテーブルの中央上空であるため、$(x_{top)}y_{top})$を極座標表示図 4:$K(t)$と$x_{\text{ゅ}}(t)$ベクトル。 丁 $=2\sqrt{d}/g$する。
$r\equiv\sqrt{x_{t\sigma p}^{2}+y_{top}^{2}}$,時系列のパワー’$\phi\equiv \mathrm{a}\mathrm{I}^{\cdot}\mathrm{C}\tan\frac{y_{top}}{x_{\ell op}}$$(0\leq\phi<2\pi)$.

$H=20d$

76$\mathrm{T}=3\mathrm{s}\mathrm{q}\mathrm{r}\mathrm{t}\{\mathrm{d}/\mathrm{g}),$$\mathrm{T}=6\{\}\mathrm{s}\mathrm{q}\mathrm{r}\mathfrak{t}(\mathrm{d}/\mathrm{g}$},$\mathrm{H}_{-}^{\wedge}30\mathrm{d}$$\mathrm{T}=3\mathrm{s}\mathrm{q}\mathrm{r}\mathrm{t}(\mathrm{d}/\mathrm{g}\rangle,$$\mathrm{H}=30\mathrm{d}$T略$\mathrm{H}=30\mathrm{d}$,$0\mathrm{s}\mathrm{q}_{\mathrm{I}}\mathrm{e}(\mathrm{d}/\mathrm{g})$$\mathrm{H}=30\mathrm{d}$図 5: \phi (右) と T(左) の時系列。 縦軸が時間であり、 $200T$ の間のデータを表示した。 $T=$
$3\sqrt{d/g}$ の場合 (上) と $T=60\sqrt{d/g}$ の場合 (下) について示す。の時間変化を図 5 に示す。 図 5 はそれぞれ 200 粒子が供給される間の時間の範囲での時系列を表したものであり、 $T=3\sqrt{d/g}$ の場合 (上) と $T=60\sqrt{d}/g$ の場合 (下) で程は時間の長さが 20 倍違い、 $T=60\sqrt{d/g}$ の の図では時間にしておおよそ 3000度のの聞頂点が同じ位置にあることがわかる。 長時間頂点が同じ位置に留まるということはほとんど雪崩が起きていないことを示す。
次に $T=3\sqrt{d/g}$ のときと $T=60\sqrt{d/g}$ のときの の時系列からパワースペクトルは広い範囲で由的になるが、 市の指数は異なを計算すると、 どちらの場合もる 図が小さいときには 5\phi (のは $1/f$ 的になり、 $T$ が大きい場合には巾の指数は一 1より小さくなる。
$T$ とパワースペクトルの巾の指数の関係を調べる。 $f=1/T$ で粒子の供給の影響でパワースペクトルにピークや段差が生じるため、 $f=1/T$ を避け範囲 10を$f\leq 1/(2T)$ でを両対数表示したグラフにたいして最小自乗法を用いて指数フィットすると、 $T$ に対して図 7 のように変化する。 $T<10\sqrt{d}/g$ では \mbox{\boldmath $\alpha$}ゆは $T$ が小さくの領域が現れる。
なるにつれて急激に一 1 に近づき、 $T>10\sqrt{d}/g$ の範囲では$T$ を変化させてもがほぼ同じ値をとる plateau 領域は 2 次元の砂山でも現れ、 $T$ が小$\phi$と$\phi,$

$r$

$r$

$\sqrt{d/g}$$\phi$$\phi$$S_{\phi}(f)$$S_{\phi}(f)$$($$6)_{\text{。}}T$$S_{\phi}(f)$$\alpha_{\phi}$$\mathrm{x}10^{-3}\leq$$S_{\phi}(f)$$\alpha_{\phi}$$\alpha_{\phi}\cong-1.3$$\alpha_{\phi}$

77
が一 1 に近づくことも 2 次元の場合で現れた特徴と同じである。 2さいときに急激に次元の場合はパワースペクトルが近似的に $1/f$ となるとき、 左右の斜面の雪崩の切り替わりの時系列も $1/f$ 的になったことより、 3 次元の砂山でも雪崩の方向の時系列を測定する。 時刻 での運動量の 成分、 成分の平均をとし、 を定義したように雪崩をの方向$\alpha_{\phi}$

$t$

$x$

$(p_{x}, p_{y})$$y$$\phi$$\theta$$(0\leq\theta<2\pi)$$\theta=\arctan\frac{p_{\tau y}}{p_{x}}$の時系列のパワースペクトルを計算すると図 8 のように、 頂点のパワースペクトルとよく一致することがわかった。
3 次元の砂山でも 2 次元の場合と同様に頂点移動と雪崩の方向の切り替わりが密接に関係していると考えられる。また、 3 次元の場合
では砂山内部の粒子を固定しシミュレーションを行ったが、 それでも 2 次元の場合と同様の結果が得られたため、 本研究で現れたような頂点の時系列の性質が粒子の供給の時間間隔 によって変化するという現象については理論やシミュレーションでは 2 次元の砂山の表面流を考えることにより解明できると考えられ、 次元によらず現れる現象であることより実験でも観測が可能であると予想される。とする。

$\theta$

$T$

5

議論
(3) で定義された $K(t)$ のような 2 値化された時系列について、 時系列の待ち時間が巾分布に従い、 それぞれの区間が独立であれば、そのパワースペクトルも巾関数で表されることが解析的に知られている $[1,3]_{0}$ ここで待ちを二値化された時系列の隣り合う切り替わりの間の時間間隔として定義する。 のを次式のような $\tau=a,$ で確率分布密度カットオフのある巾関数とする。
時間$\tau$$\tau$$b$$p(\tau)$$p(\tau)=\{$

但し、$a,$$b$$c\tau^{D}$$a<\tau<b$ のとき

0

それ以外

図 7:

$\alpha_{\phi}$と$T(H=3\square \mathrm{J}d)_{0}$(4)は定数で $a<<b_{0}c$ は規格化定数。 $1/b<<f<<1/a$ の範囲で二値化された時系列のパワースペクト)$\triangleright S_{b}(f)$は$S_{b}(f)\sim f^{-(D+3)}$,$(-3<D<-1)$$(_{\backslash \mathrm{J}}^{r})$78となる。 $T$ が小さい場合、 シミュレーションでは $K(t)$ のパワースペクトルの巾の指数は近似的に一 1 であるため、 $p(\tau)$ は指数 $D\cong-2$ の巾関数となることが予想される。 しかしながら $K(t)$ の待ち時間は短いノイズによって寸断されるために、 $\tau\cong 100\sqrt{d/g}$ より長い待ち時問間隔はほとんど現れない。 それゆえに時問幅 $60\sqrt{d}/g$ の median fflter を利用し、 粗磁化された時系列の持ち時間分布 $p(\tau’)$ を計算した。 $p(\tau’)$ は図 9 のように近似的に巾則に従い、 その指数は一 2 であることがわかった。

頂点のパワースペクトルの巾の指数-1.4 の plateau 領域は図 3, 図 7 で の広い範$\alpha\cong$$T$

囲で現れる 9 しかしながら、 $T$ が大きくなる

につれて雪崩のモードの切り替わりは $K(t)$ によって定義しにくくなる。頂点のゆらぎは粉体流の密度ゆらぎとなんらかの関係がある可能性がある。 流体でみたされた鉛直パイプ内での粉体流の実験では密度の時間のパワースペクトルについて一 1, -4/3, -3/2 が報告されて
おり [4, 5, 6, 7, 8, 9, 10, 11]、交通流 ASEP では一 1.4, -3/2 が報告されている $[12, 13, 14]_{0}$これらの指数は -1.4 t こ近い。
は一 2 に近が十分大きくなるにつれて付くと思われる。頂点は $T$ に比べて十分長い$T$$\alpha$

図 8:$\phi$と$\theta$の時系列のパワースペクト

ル。 $(H=30d, T=3\sqrt{d/g})$

時間ほぼ同じ場所に留まり、 頂点移動は雪崩の瞬間的な力によって引き起こされる。 このとき頂点はランダムウォークのように動き、 そのパワースペクトルは Lorentzian 的での減衰が現れると予想される。

$f^{-2}$6まとめ

砂山形成過程は雪崩と供給された粒子の積み上げによって成り立ち、 雪崩の発生の仕方は粒子の供給量によって大きく変化する。 本
研究では砂山頂点の移動を通して、 形成過程を考えた。頂点の位置は雪崩によって雪崩と逆の方向に移動し、 また粒子が雪崩を経験せずに積み上がることによって供給位置に引き寄せられるという性質をもつことから、 頂点の位置のダイナミクスは形成過程と密接に結びついている。 本研究ではシミュレーション

図 9: $H=11\mathrm{O}d$ で $T=2\sqrt{d}/g$ のときの $K(t)$ の待ち時間分布。

により 2 次元と 3 次元の砂山の頂点移動について調べた。

砂山上空の 1 点から蒔間間隔$T$で 1 粒子づっ落とすことによって粒子の供給を行うこ79ととした。 $T$ が小さいときには供給量が大きいため常に砂山斜面のどこかで粒子が雪崩ている状態を実現し、 逆に $T$ が大きいときには粒子供給がほとんどないため表面は長い間固体的状態を保ち、 雪崩は間欠的に発生する。 頂点の時系列を測定しパワースペクト
ルを計算すると、 $T$ によらず広い波数領域で巾的になり、 巾の指数は $T$ によって変化することがわかった。 $T$ が非常に小さいときにはパワースペクトルは $1/f$ 的になり、 $T$ がやや大きくなると巾の指数は減少し、 $T$ がある値を越えて大きくなると $T$ の広い範囲で近似的に同じ指数をもつ巾に従い、 2 次元では $T$ をさらに大きくすると指数が更に減少する。 $T$ の変化に対する頂点の時系列の性質の変化は粒子の供給量によって砂山形成過程の様相がかわることと対応していと考えられる。 また、 $T$ が小さいとき系の運動エネルギーや運動量を用いて雪崩を計測し、 雪崩の方向の時系列を生成すると、 そのパワースペクトルは頂点の時系列のパワースペクトルとよく一致し $1/f$ 的になることがわかった。 2 次元では頂点からすそのに下る 2 つの斜面しかないため、 雪崩の方向の時系列は2 値化された時系列 $K$ で書き表すことができ、 2 値化された時系列について時系列のパワースペクトルが巾則に従い時系列の待ち時間分布も指数 $D$ の巾関数で表されるときパワースペクトルの巾の指数は一 $D-3$ と表されることが解析的に求められており、 これよリパワースペクトルの指数が一 1 のときに待ち時間分布の指数が一 2 であることが予想される。 実際に $K$ の待ち時間分布を測定すると、 近似的に指数一 2 の巾が現ることがわかった。 このことより、 $K(t)$ の時系列のパワースペクトルが $1/f$ 的になることは待ち時間分布が指数 -2 の市関数がであることから説明された。

参考文献
Teich. Fractal renewal process generate[1]$E,$$1/\mathrm{f}$noise. Phys. Rev.

Vol. 47, p. 992, 1993.

[2] P. A. Cundall and 0. D. L. Strack. A discrete numerical model for granularblies. Geotechnique, Vol. 29, pp. 47-65, 1979.

$\mathrm{a}_{\wedge}\mathrm{s}\mathrm{s}\mathrm{e}\mathrm{m}-$[3] Nobuko Fuchikami and Shunya Ishioka. Statistics of level crossing intervals. InZolt\’an Gingi, Jose’ M. Sancho, Lutz Schimansky-Geier, ar Janos Kertesz, editors,
Proceeding of Noise in Complex Systems and Stochastic Dynamics $II,$ Vol. 5471 ofProc. SPIE, pp. 29-37. SPIE, 2004.
$\iota \mathrm{d}$[4] Shin Horikawa, Akio Nakahara, Tomoo Nakayama, and Mitsugu Matsushita. Selforganized critical density waves of granular material flowing through a pipe. J. Phys.
Soc. Jpn., Vol. 64, pp. 1870-1873, 1995.
[5] Shin Horikawa, Takeshi Isoda, Tomoo Nakayama, Akio Nakahara, andsushita. Self-organized critical density wave of granular particles flowing through apipe. Physica $A,$ Vol. 233, pp. 699-708, 1996.
$\mathrm{M}\mathrm{i}\mathrm{t}_{\mathrm{S}\mathrm{u}\mathrm{g}}1\mathrm{M}\mathrm{a}\mathrm{t}\sim$80[6] Akio Nakahara and Takeshi Isoda.
density fluctuation at the slugging transition
point of granular flows through a pipe. Phys. Rev. , Vol. 55, pp. 4264-4273, 1997,$1/f^{\alpha}$$E$[7] Osam Moriyama, Naoya Kuroiwa, Mitsugu Matsushita, and Hisao Hayakawa. 4/3. Left., Vol. 80,
law of granular particles flowing through a vertical pipe. Phys.
pp. 2833-2836, 1998.
$\mathrm{u}$$f_{\iota}^{)}ev$[8] O. Moriyama, N. Kuroiwa, T. Isoda, T. Arai, S. Tateda, Y. Yariiazaki, anMatsushita. Dynamics and structure of granular flow through a vertical pipe. In
M. Fukui, Y. Sugiyama, M. Schreckenberg, and D. E. Wolf, editors, TRAFFICAND GRANULAR, pp. 437-448. Springer, 2001.
$1\mathrm{d}\mathrm{M}$$FI,OW’ \mathit{0}\mathit{1}$[9] Yoshihiro Yamazaki, Sayako Tateda, Atsushi Awazu, Taichi Arai, Osamu Moriyarna,and Mitsugu Matsushita. Dynamics of granular flow through a vertical pipe undertle control of infiow rate. J. Phys. Soc. Jpn., Vol. 71, pp. 2859-2862, 2002.
$\lfloor \mathfrak{f}10]$

Osamu Moriyama, Naoya Kuroiwa, Sayako Tateda, Taishi Arai, Atsushi Awazu,
Yoshihiro Yamazaki, and Mitsugu Matsushita. }to density waves of granular particles flowing through a vertical pipe. Prog. Theor. Phys. Supp., Vol. 150,
pp. 136-146, 2003.
$7$

[11]

.

of granular particles flowing through a vertical glass pipe.
thesis, Ch uo University, 2005.

$\mathrm{A}\_{J}\mathrm{s}\mathrm{u}\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{A}\mathrm{w}\mathrm{a}^{r}A11$$\mathrm{M}\mathrm{a}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{r}.\mathrm{s}^{\backslash }$$\mathrm{a}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$\mathrm{D}\mathrm{y}_{\mathrm{I}1}\mathrm{a}\mathrm{r}\mathrm{n}\mathrm{i}\mathrm{c}\mathrm{s}$[12] Toshimitsu Musha and Hideyo Higuchi. The
fluctuation of a traffic current on
an expressway. , $fpn$ . J. Appl. Phys.. Vol. 15, p. 1271, 1976,$1/\mathrm{f}$[13] Toshimitsu Musha and Hideyo Higuchi. Traffic current fluctuation and the burgers
equation. $Jpn$ . J. Appl. Phys., Vol. 17, pp. 811-816, 1978.
[14] Shinji Takesue, Tetsuya Mitsudo, and Hisao Hayakawa. Power-law behavior in thepower spectrum induced by brow nian motion of a domain wall. Phys. Rev. , Vol. 68,
p. 015103(R), 2003.
$E$