第 24 回数値流体力学シンポジウム
C5-3

岩石風食の数値シミュレーション
Numerical Simulation of Wind Erosion for Rock○ 桑名 杏奈, お茶大院, 東京都文京区大塚 2-1-1, E-mail: g0770610@edu.cc.ocha.ac.jp
秋保 美幸, お茶大院, 同上
河村 哲也, お茶大院, 同上
Anna KUWANA, Ochanomizu University Graduate School of Humanities and Sciences,2-1-1 Ootsuka, Bunkyo-ku, Tokyo
Miyuki AKIHO, Ochanomizu University Graduate School of Humanities and Sciences
Tetuya KAWAMURA, Ochanomizu University Graduate School of Humanities and Sciences
To reproduce strange shaped rocks appeared to Cappadocia in Turkey, we made a simple model of simulating erosion
by wind using a sand transfer equation. The numerical method is as follows: First, flow field around a rock which isgoverned by incompressible Navier-Stokes equations is solved numerically by the standard MAC method. Second, we
estimate amount of erosion that is proportional to the cube of friction velocity. Third, new shape of an eroded rock iscalculated. These three steps are repeated. The flow around a rock and its surface are visualized and the influences ofhardness of the rock on erosion by wind have been analyzed.
１．はじめにカッパドキア（Cappadocia）は 1985 年に世界遺産に登録された，トルコ中央部の岩石地帯である（Fig. 1）．海抜 1,000m のアナト
リア高原中央部に 100km2 に渡って広がるこの地帯には，Fig. 2 に示すようなキノコやラクダの形をした「奇岩」が林立している．
付近の火山（Mt. Erciyes，Mt. Hasan）から噴出された火山灰から成る凝灰岩層は，玄武岩などに比べて侵食を受けやすい．
奇岩は，岩石種の違う地層が複雑に重なりあっている地域における侵食の速度の違いにより，生成されたと考えられている．
本研究では，カッパドキア地方の奇岩形状を再現することを目標に数値シミュレーションを行った．
侵食に影響を与える要素（温度差・湿度差による風化，水や氷による侵食，水・塩類による凍結破砕，溶解・酸化などの化学的風化など）全てを考えるのは難しいため，単純なモデルとして風による侵食，特に岩石表面の粒子が飛散する現象（デフレーション）をとりあげる．
デフレーションによって生成される地形の典型が，砂丘である．
先行研究(1), (2) は，バルハン砂丘の生成など，風による砂の動きを数値シミュレーションにより解析したものであるが，今回はそこで提案された手法を岩石の風食に応用した．

２．計算方法２－１．計算領域
円筒形および円錐台形の岩石を，Fig.3 に示す円筒座標系の格子を用いて表現する．格子数は，周方向 122 x 半径方向 41 x 高さ
方向 26 = 約 130,000 点とした．遠方境界の風上側から変動風が流入し，風下側へ流出する．風は，大きさは一定であり，方向は Fig.
4 に示すように x 軸方向を主流方向として y 軸方向に±30 度，正弦関数的に 5,000step あたり 1 周期変動するものとする．計算領域の上下面は滑りなし壁とする．
岩石側面が，岩石周りの流れ場により削られると考える．岩石形状の上下境界は外挿で決定した．岩石面に沿った物理空間は侵食が進むにつれて複雑になるため，一般座標変換を施して得られる直方体領域の計算空間にて計算を行う．

(a) A rock of circular cylindrical shape.

Fig.1 Geolocation of Cappadocia.

(b) A rock of truncated-circular-cone shape.
Fig. 3 Computational grids.

(a) Mushroom shape.
(b) Camel shape.
Fig.2 Strange shaped rocks.
Fig. 4 Periodic change of wind.
1Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C5-3
２－２－３．新しい岩石面形状の決定
各格子点における岩石侵食量と，岩石表面の微小面積における砂の質量保存（Fig. 5）から，式(9) が成り立つ．式(9) は，岩石侵食量ベクトルから，岩石面の高さの時間変化が計算できることを意味している．

２－２．計算方法
２－２－１．岩石周りの流れ場を計算
下記の基礎方程式により岩石周りの流れ場を計算する．
連続の式： ∇ ⋅ u = 0(1)非圧縮性 Navier-Stokes 方程式：∂u1 2(2)+ ( u ⋅ ∇)u = −∇p +∇u∂tRedq dqdh=− 1 − 2(9)dt
dX dZ
ρ s ：砂密度，h ：岩石面に垂直な高さ，( X , Z ) ：岩石面に平行な局所座標，(q1 , q2 ) ：岩石侵食量ベクトル q（式(7)）の ( X , Z ) 成分ρs

∇ : 勾配演算子， u = (u, v, w) : 速度ベクトル，t : 時間，p : 圧力， Re : レイノルズ数 (=1,000)
基礎方程式(1), (2) は MAC 法(3) を用いて解く．
∂u1 2+ (u ⋅ ∇)u = −∇p +∇u∂tRe22(3)2∂u ∂v∂w ∂u∂v ∂w ∇ ⋅ u
 ∂u   ∂v   ∂w ∇ 2 p = −  −   −   − 2−2−2+∂y ∂x∂x ∂z∂z ∂y∆t ∂x   ∂y   ∂z (4)時間微分は一次精度の前進差分，非線形項以外の空間微分は二次精度の中心差分，非線形項は充分に細かくない格子でも安定した計算が可能(4) な三次精度の上流差分(5) で近似した．
f∂u− ui + 2 + 8(ui +1 − ui −1 ) + ui − 2x= x = f∂x i12∆x
f ui + 2 − 4ui +1 + 6ui − 4ui −1 + ui −2+12∆xFig. 5 Mass conservation of sand.
式(7) を式(9) に代入し，無次元化する．式(10) において「~」は無次元量を意味する．

(5)~dh~dt

２－２－２．表面摩擦の計算と，岩石侵食量の推定
Bagnold は，風による表面摩擦速度の大きさ u* [m/s] と砂粒の跳躍による砂輸送量 q [kg/(m･s)] の間の定量的な関係式：砂輸送
方程式(6) を導いた(5) ．
ρ
q = b1 0 u *3g(6)，b1 ：砂面形状による定数（本研究では一定とする）
ρ 0：空気密度（定数）， g ：重力加速度

dU U
dn U

(10)



ここで，式(11)によりパラメータ C を適宜する．b1, ρ 0 , Fr は定数のため，C は ρ s（砂密度）によって決定される．
ρ(11)C = −b1 0 Fr 2ρs

本研究では，砂密度を岩石密度と読み替える．岩石密度は，玄武岩で 3.0[g/cm3]，凝灰岩で 1.4～1.5[g/cm3]，軽石で 1.4[g/cm3]程度である．岩石密度の大きい（パラメータ C の値が小さい）ものを侵食されにくい岩，岩石密度の小さい（パラメータ C の値が大きい）ものを侵食されやすい岩とする．
また，岩石表面の凹凸の大きい部分は粗い格子では表現しにく
いこと，およびそういった部分は速く侵食されると考えられることから，計算においては下記の平滑化を施している．
・ 各格子点における岩石面の曲率が一定値を超えたら，その格子点での岩石面高さ hi ,k を隣接点における岩石面高さの平均値で上書きする．
 周方向

摩擦速度ベクトルu*は，式(8) で表される．
ν

~Fr ：フルード数

式(6) は，砂輸送量が摩擦速度の大きさの 3 乗に比例することを意味する．実際には，砂面上空の流れ場は場所によって様々な方向をもち，砂は各地点における風の強さに従った量がその風の向きに運ばれる．すなわち，砂輸送量は場所によって異なる方向と大きさをもつベクトル q で表される．
ρ2q = b1 0 u* u *(7)gu* =

ρ 0 2  d q1 d q2 Fr+ ~~ρs d X d Z ~= −b1(8)ifν ：空気の動粘度（定数），hi +1,k − 2hi ,k + hi −1,k> 2%
2hi +1,k + hi −1,k2ofRthenhi ,k ⇒of

R

then

hi ,k ⇒

 高さ方向U ：砂表面に平行な方向の速度d dn ：砂面に垂直な方向の微分ifhi ,k +1 − 2hi ,k + hi ,k −1> 10%
2hi ,k +1 + hi ,k −12（R は岩石上部半径の初期値）
・ ある格子点において，岩石半径が初期値 R の 40%になったら，その地点では侵食を終了する．

本研究では岩石を砂の塊と仮定し，「岩石の侵食」を「岩石表面の砂の消失」と考える．式(7) の砂輸送量ベクトルを「岩石侵食量ベクトル」と読み替えて，以下の計算に利用する．

式(10) により新しい岩石面形状が決定された後に，それに沿うように格子を作り直す．そして２－２－１．に戻って新しい岩石周りの流れ場を計算する．これを繰り返す．ただし，岩石が侵食を受ける時間スケールは流れ場が変化する時間スケールに比べて大きいので，流れ場に関する計算（２－２－１．
）100step ごとに，岩石に関する計算（２－２－２． ２－２－３．
）を 1 回行った．
2Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C5-3３．結果
３－１．円筒形の岩石はじめに，岩石の硬さを表す式(11) のパラメータ C の値の違いによる侵食の受け方の違いを比較するために，一様な硬さをもつ2 種類の円筒形の岩石について計算を行った（計算領域は Fig. 3 (a)に示したもの）．Fig. 6 は岩石面のシェーディング（カラーマップは岩石半径）と，高さ方向中央断面における流速ベクトルを示したものである．Fig. 6 (a), (b) どちらも同じく，計算が始まってから 10,000step での様子であるが，パラメータ C の値が大きい(a) の方が，より侵食を受けている様子が観察される．

３－２．奇岩の再現１
Fig. 7 (a) の写真にある岩を再現するために，パラメータ C の値が違う 3 層が重なった円筒形の岩石 Fig. 7 (b) について計算を行った．5,000step での流れ場と岩石面形状を Fig. 7 (c) に示す．パラメータ C の値が大きい LayerB が，強い侵食を受けている．

(a) A strange shaped rock.

LayerA: C = 0.5 x 10-3
LayerB: C = 0.5 x 10-2
LayerC: C = 0.5 x 10-4
Height of layers A: B: C= 2: 3: 25(b) Three layers.

Plane view

Plane view

Side view(a) C = 1.0 x 10-3 (Soft rock)

Side view

Plane View

Front view(c) Flow fields and the surface of the rock of consisting of three layers.
Fig. 7 Reproduction of a strange shaped rock 1.

Side view(b) C = 0.5 x 10-3 (Hard rock)Fig. 6 Flow fields and the surface of the rocks ofsimple circular cylindrical shape.

3Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C5-3３－３．奇岩の再現２
Fig. 8 (a) の写真にある岩（きのこ岩．カッパドキア現地では妖精の煙突と呼ばれている）を再現するために，パラメータ C の値が違う3層が重なった円錐台形の岩石Fig. 8 (b) について計算を行った（計算領域は Fig. 3 (b) に示したもの）．10,000step での流れ場と岩石面形状を Fig. 8 (c) に示す．パラメータ C の値が大きい
LayerB が，強い侵食を受けている．

(a) A strange shaped rock.

３－４．奇岩の再現３
Fig. 9 (a) の写真にあるような岩石のくぼみを再現するために，一部でパラメータ C の値が大きい円錐台形の岩石 Fig. 9 (b) について計算を行った．10,000step での流れ場と岩石面形状を Fig. 9 (c)に示す．LayerA と LayerB の間にくぼみが再現された．

LayerA: C = 0.5 x 10-3
LayerB: C = 0.5 x 10-2LayerC: C = 0.5 x 10-4
Height of layers A: B: C= 2: 1: 6(b) Three layers.

(a) Concaves of rock.

LayerA: C = 0.5 x 10-3
LayerB: C = 0.5 x 10-2
LayerC: C = 0.5 x 10-4
Height of layers A: B: C= 2: 1: 6(b) Three layers.

Plane view
Plane view

Side view
Side view

Front view(c) Flow fields and the surface of the rock of consisting of three layers.
Front view(c) Flow fields and the surface of the rock of consisting of three layers.

Fig. 9 Reproduction of a concave of rock.

Fig. 8 Reproduction of a strange shaped rock 2.

4Copyright © 2010 by JSFM

第 24 回数値流体力学シンポジウム
C5-3
４．まとめ
岩石を砂の塊だと仮定し，「岩石の侵食」を「岩石表面の砂の消失」と考えることにより，砂輸送量ベクトルを岩石侵食量ベクトルと読み替えて，単純なモデルにおける岩石風食の過程を計算した．岩石の硬さと層の構造の違いにより，カッパドキア地方にある奇岩の形状をある程度再現できた．今後は今回の単純なモデルを，より現実に則した形に修正していく予定である．

参考文献(1) M. Kan, "Application of numerical simulation to someenvironmental problems in arid land", Dissertation, Ochanomizu
University (2000).
(2) R. Zhang, M. Kan and T. Kawamura, "Numerical study of theformation of transverse dunes and linear dunes", Journal of the
Physical Society of Japan, Vol. 74 No.2 (2005) pp.599-604(3) F. H. Harlow and J. E. Welch, "Numerical calculation oftime-dependent viscous incompressible flow of fluid with freesurface", The Physics of Fluids, Vol. 8, No. 12 (1965), pp.109-116.
(4) T. Kawamura and K. Kuwahara, "Computation of high Reynolds
number flow around a circular cylinder with surface roughness",
AIAA Paper 84-0340 (1984).
(5) R.A.Bagnold, "The movement of desert sand", Proc. Roy. Soc.
A157 (1936), pp594-620.

5Copyright © 2010 by JSFM