円柱群を過ぎる流れにおける
undisturbed flow
の基本特性
FUNDAMENTAL PROPERTIES OF UNDISTURBED FLOWIN ARRAYS OF CIRCULAR CYLINDERS
横嶋 哲1) 平井 俊也2) 安藤 丈央3) 風早 祐輔3) 内田 龍彦4) 河原 能久5)
Satoshi YOKOJIMA1), Toshiya HIRAI2), Takefumi ANDO3), Yusuke KAZEHAYA3), Tatsuhiko UCHIDA4) and Yoshihisa KAWAHARA5)
ABSTRACT
A drag-force model is essential in practical predictions of flows around vegetation/urban canopies, since it is too prohibitive to resolve details of the canopy elements and asso-ciated fluid motions, especially in geophysical applications. The model, however, has a serious difficulty inherent in its formulation: how one can specify the representative velocity? In ideal situations such as flows past an isolated obstacle, it is the velocity of the inflow and is easily available. In flows past multiple obstacles, on the other hand, a straightforward extension of the idea behind the ideal situations leads to introduce undisturbed flow, which is the flow that would exist at an obstacle location in the ab-sence of that obstacle but with all other obstacles present. Here the undisturbed flow has been evaluated directly in fully resolved computations of a two-dimensional flow past circular cylinders and the fundamental properties of the flow are discussed. Key Words: Canopy flow, Drag-force model, Undisturbed flow, CFD
1. はじめに 都市域における風環境や氾濫流解析，あるいは樹木群を有する大気境界層や河川流の予測において，個々の 構造物や植生要素のスケールは流れの代表スケールと比較して圧倒的に小さいため，その周囲の流動を個々の 要素形状を織り込んで再現するには工学的に許容できないコストを要することが多い．したがってそれらキャ ノピー層内の流動と抗力の関係を何らかの巨視的パラメータを用いて適切に表現・モデル化することが強く望 まれる．一般に見付面積A，抗力係数CDの任意の物体に作用する力Fは，ニュートンの抵抗則に基づいて
F = (1/2)CDρA|urep|urep (1)
で表される．ここでurepは物体に対する流体の接近速度（fluid velocity seen by object）の代表値，ρは流体
密度を表す．よって，この物体が単位体積当たりの流体に及ぼす力fは f =−F/V = −(1/2)CDρλ|urep|urep (2) ————————————————————————————————————————— 1) 静岡大学 学術院工学領域 数理システム工学系列 准教授（〒432-8561浜松市中区城北3-5-1） 2) 静岡大学 大学院総合科学技術研究科 工学専攻 数理システム工学コース 大学院生（同上） 3) 静岡大学 工学部 数理システム工学科 大学生（同上） 4) 広島大学 大学院先進理工系科学研究科 社会基盤環境工学プログラム 准教授（〒739-8527東広島市鏡山 1-4-1） 5) 広島大学 学術・社会連携室 副学長（同上）
U0 U0 U0 U0 FD=(1/2)ρAurep2 FD=(1/2)ρAurep2 urep=U0 urep=? (a) (b) (c) (d) 図1 式(1)に基づいて巨視的に円柱抗力を評価する際に必要となる代表流速urep．(a)孤立円柱を過ぎる流れ， と(b) (a)の円柱に対する代表流速の考え方（当該円柱を除去した場合に当該円柱位置を過ぎる流速）． (c)円柱群を過ぎる流れ，と(d) (c)のグレーで色付けされた円柱に対する代表流速の考え方（当該円柱は 存在せず，それ以外の円柱は全て存在する場合に，当該円柱位置を過ぎる流速）． のように近似でき，いわゆる抗力モデルとして，これまでに多くの適用例が報告されてきた1)．ここでV は式 (2)を適用する検査体積を表し，λ≡ A/V は障害物の表面積密度と呼ばれる． 障害物群が幾何学的に単純であれば表面積密度λはその定義に基づいて一意に決まるのに対して，適切な 抗力係数CDを導く標準的な手法は単純な障害物群に対しても存在せず，このアプローチの大きな課題であ る．Yokojima et al.2)は円柱群（樹木群模型）を過ぎる直線開水路流れのラージ・エディ・シミュレーション （LES）において，円柱群内部の抗力係数としてある値を一様に与えた場合，円柱群が主流方向に連続するケー
ス（Case 1）でのLESの予測精度は高い反面，円柱群が主流方向にパッチ状に配置された場合（Case 2）の流
れ特性の再現性は不十分であることを報告した．彼らは抗力係数CDの値を円柱群内部で0.6, 0.8, 1.0と一様 に与えた3種類のLESをそれぞれのケースに対して実行し，対応する水理実験データとの比較によってCase 1ではCD= 0.6，Case 2では1.0がより良い結果を与えることを示した．このように，同じ障害物群でも，そ の配置パターンによって適切な抗力係数の値は変わり得ることが問題を難しくしている． 対象流れに対してどのような観測値も存在しない場合にも，適切な抗力係数CDを推定する方法の確立を目 指して，横嶋・河原3,4)は以下のようなアプローチを試みた．すなわち，上述の円柱群を過ぎる直線開水路流 れ2)の幾何学的最小ユニットを対象として，個々の円柱周りの流動および各円柱に作用する流体力を支配方程 式に基づいて近似無しに算出した数値シミュレーション（文献3,4)では予備解析と呼称）結果から円柱群内部 のCD分布を推定し，得られたCD分布を3次元LES解析（文献3,4)では主解析と呼称）に導入することで， 文献2)で得られた予測結果と比べてどれほどの精度向上が可能かを調べた．個々の円柱の抗力係数C Dの評価 に必要な代表速度スケール（式(1)中のurep）として，横嶋・河原3)は断面平均流速Ubを全ての円柱に対し て一様に用いたのに対して，横嶋・河原4)は個々の円柱に対する接近流速u aを個別に算出した（uaの定義等 は3.1節で詳述）．水理実験結果2)と比較した場合に，接近流速u aに基づく抗力係数CD,localは円柱群抵抗を 過大に，断面平均流速Ubに基づくCD,globalはやや過小に評価する傾向が認められたものの，特にCD,globalに ついては抗力係数を一様に与えた場合2)の最も良好な結果と同程度の予測が可能であることが示された． 上述のように，著者らのグループではこれまで，適切な抗力係数の推定法を求めていくつかの試みを行って きた．本研究ではその一環として，特に代表速度urepに注目する．流入流速U0を有する系に孤立円柱が存在 し，この円柱に働く抗力FDを式(1)で評価する場合を考える（図1(a)）．代表速度urep（urepの主流方向成 分）は円柱が存在しない場合にその円柱位置を過ぎる流速（当該円柱の影響を受けない場合の流速であるので，
以下ではundisturbed velocity uunと呼称）であり，この場合，明らかにurep= U0となる（図1(b)）．孤立
円柱のCD値もよく知られているので，抗力FDは式(1)を用いて容易に評価できる．次に，円柱群中のある
円柱に働く抗力を評価する場合を考える（図1(c)）．この場合，先述の孤立円柱を過ぎる流れにおける代表流
速の考え方を拡張すれば，当該円柱に対する代表流速（undisturbed velocity）は，その円柱は存在せず，それ 以外の円柱は全て存在する場合に当該円柱位置を過ぎる流速となり，一般にurep̸= U0であり，その算出は全 く容易ではない（図1(d)）．このように抗力モデル（式(1)）は，実用的な問題への適用を考えた場合に，本
x1 x3 Vegetation zone 27cm26.5cm 80cm 198cm (a) 26.5cm 99cm 99cm 198cm (b)
図2 検討対象とした円柱群（樹木群模型）を有する直線開水路流れの概要: (a) Case 1; (b) Case 2．破線の 長方形領域は2.2節に記載の数値解析の計算対象領域（図3参照．Case 2の幾何学的最小ユニット）を 示す． 質的な困難を有している． Undisturbed velocity uunが必要となる別の例を以下に挙げる．微小粒子の流体中での振る舞いを追跡する 際に，流体運動をオイラー視点で，個々の粒子運動を質点近似に基づくMaxey-Riley方程式5)を用いてラグラ ンジュ視点で表現することはよく行われる6)．微小粒子に働く流体抵抗は，慣性が非常に強い場合にはストー クス抵抗のみで精度良く近似できるため，粒子（半径a，質量mp）および流体（密度ρf，粘性係数µ）に対 する運動方程式は，例えば以下のように表される． mp du(k)p dt = 6πµa ( u(k)rep− u(k)p ) � �� � F(k)St , (3) ρf Du Dt =−∇p + µ∇ 2u− 1 Vc ∑ k F(k)Stδ(x− x(k)p ). (4) ここでx(k)p , u(k)p はk番目の微小粒子の位置と速度，uとpは流体流速と圧力を表す．式(4)の右辺第3項は 各粒子が受けるストークス抵抗の流体への反作用を表し（Vcは式(4)に関わる参照体積，δはデルタ関数を示 す．），粒子の存在が流体場に影響しない（いわゆる1wayカップリング）場合には，無視される．この場合， k番目の粒子に働くストークス抵抗F(k)St の算出に必要な代表流速には，式(4)から求まる流速の粒子位置での 値を与えればよい（すなわち，u(k)rep= u ( x(k)p ) ）．つまり，1wayカップリング解析では，代表流速の算出は空 間補間の問題に帰着する．他方で粒子の流体場への影響が無視できない場合（2wayカップリング），式(4)か ら求まる流速uは全ての粒子の影響を受けたものであり，他方でk番目の粒子に働くストークス抵抗F(k)St の 算出に必要な流速は，k番目の粒子は存在せず，それ以外の粒子は全て存在する場合に，k番目の粒子位置を 過ぎる流速であり，その評価方法は自明ではない．つまり，2wayカップリング解析では，全ての粒子の影響
を受けた流れ場（disturbed flow field）から，特定の粒子のみの影響が排除されたundisturbed velocityを推
定する必要がある．この問題の解決策が近年になっていくつか提案された7−9)ものの，その試みはまだ緒に 就いたばかりであり，先行研究の大半では，2wayカップリング解析においても1wayカップリングと同様に， 流速場uの単なる粒子位置への補間値が代表流速とされているのが現状である10)． 本研究では，先行研究2−4)で検討がなされた円柱群を過ぎる完全発達流れに関して，個々の円柱を除去した 系を対象とした新たな数値実験を実施することでundisturbed flowを直接的に求め，その基本特性を検討した． 2. 検討対象流れとその数値解析法 2.1検討対象流れ 全長24 m，幅80 cm，水路勾配1/555の直線水路上に，中心間距離sを3 cmとして正方格子状に配置された 直径D = 3 mmの円柱からなる円柱群（樹木群模型，円柱群の表面積密度λ = 1/30 cm−1）を過ぎる流れ2)を 検討対象とした．水路床および側壁は水理学的滑面である．詳細はYokojima et al.2)を参照されたい．概略図を 図2に示す．幅27 cmの円柱群を主流方向に連続的に水路中央に沿って配置した場合をCase 1，99 cm× 27 cm
の円柱群パッチを99 cm間隔で主流方向に離散的に配置した場合をCase 2とする．図2中の‘vegetation zone’
（グレーで塗られた領域）は，各円柱が3 cm× 3 cm × H（ここでHは水深を表す）の領域を受け持つと考え
た場合のキャノピー部を表し，このキャノピー部を基準とした場合の樹木群密度λ（式(2)に含まれる），plan area index，およびfrontal area indexの値はそれぞれ(1/30) cm−1，0.0079，0.1となる．いずれの場合にも
x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (a) x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (b)
図3 水平2次元モデルにおける計算領域および円柱配置: (a) Case 1; (b) Case 2．主流x1方向には周期条件 を課した．赤色の円柱位置におけるundisturbed flow uunを算出するため，それらを個別に除去した系 を対象として数値実験を行った．見易さのため，円柱は実際のサイズよりも大きく描画した．
表1 水平2次元モデルに基づく数値実験の計算条件．
ReD≡ Domain size # of grid points Grid spacing
UbD/ν L1/D L3/D n1 n3 ∆x1/D ∆x3/D ∆tUb/D All runs 700 660 800/3 7920 1274 1/12 1/12∼ 20/12 1/64 量を9000 cm3/sとしたときの平均水深HはCase 1の場合に4.78 cm，Case 2の場合に4.71 cmとなり，円柱 群は常に非水没な状態にあった．断面平均流速Ubと平均水深に基づくレイノルズ数は11250である． 2.2数値解析法 室内水理実験結果2)によれば，検討対象流れの平均流構造は水路床近傍を除けば水深方向に概ね一様であ る．よって本研究ではこの流れに対して図3に示すような水平2次元モデルを導入する．この水平2次元モデ ルの計算領域はCase 2における幾何学的な最小ユニットに当たり，どちらのCaseにおいても主流方向に周期 境界条件（周期長：198 cm）を課した． 数値実験においては，全てのケースにおいて，断面平均流速Ubと円柱直径Dに基づく円柱レイノルズ数 ReDが常に700となるように，流れの駆動力fdriveを瞬時レベルで調節した11)（後述の図4(b)参照）．この条 件は室内水理実験における水表面近傍の流況（ReD≈ 790）に概ね一致する．孤立円柱を過ぎる流れでは，レ イノルズ数ReDが200を超えたあたりから，たとえ流入流れや円柱形状が円柱軸方向に一様であっても，円 柱後流が3次元性を示すことはよく知られている12)．しかしながら，この知見が円柱が複数近接する場合に どの程度成立するのかは不明である．よってここではゼロ次近似として水平2次元モデルを採用する． 円柱の存在が流動に及ぼす影響の表現や，円柱に働く流体力の評価には，Uhlmann13)によって提案された 埋め込み境界法を用いた．この埋め込み境界法では剛体表面に分布する点の集合で剛体形状を表し，それぞれ の点上で粘着・不透過条件が課される．流体解析は変数をスタガード配置したデカルト計算格子上で行った． 計算格子間隔は主流方向には一様に∆x1 = D/12を与え，流路横断方向には円柱群内部と流路側壁近傍で最 も密な∆x3 = D/12とし，それ以外の領域では最大で∆x3= 1.67Dとなるような不等間隔な計算格子を用い た．計算条件に関する情報を表1にまとめた．本研究で実施する数値解析の手法や物理/計算条件は，横嶋・河 原3,4)において予備解析と呼ばれたものと同一である．違いは，2.3節で詳述するように，本研究では各円柱 位置におけるundisturbed flowを直接的に算出するために，円柱を個別に除去した新たな系を対象として数値 実験を繰り返した点にある． 図4はどの円柱も除去しない本来の対象流れの数値実験から得られた，運動エネルギーと流れの駆動力の時 系列分布を示す．Case 1では流速の初期値を，円柱内部でゼロ，それ以外で断面平均流速Ubを主流方向成分 として与えた．Case 2ではCase 1のt = 650D/Ubの瞬時場を初期条件として課した．図4の分布から各々の 系が統計的定常状態に達した時刻を判断し，それ以後のデータ（図4の破線部，Case 1では5000D/Ub，Case
2では8000D/Ub時間に渡る）を用いて統計量を算出した．その一例として，平均主流速分布を図5に示す．
2.3 Undisturbed flowの直接評価
tUb=D 0 5000 10000 K = U 2 b 0.5 0.7 0.9 7 K = 0:895Ub2 7 K = 0:801Ub2 tUb=D 0 5000 10000 (fd ri v e D =U 2 b) # 1 0 ! 3 0 1 b 7 fdrive= 8:15 # 10!4(Ub2=D) 図4 どの円柱も除去しない本来の検討対象流れの数値実験から得られた時系列分布: (a)単位質量当たりの運 動エネルギーKの系全体での平均値; (b)単位質量当たりの流れの駆動力fdrive．青と赤の実線はそれぞ れCase 1とCase 2の結果を，破線は統計量算出のためのデータ収集時間を，f¯はfの時間平均値を示す． 図5 どの円柱も除去しない本来の検討対象流れの数値実験から得られた平均主流速分布u¯1/Ub: (a) Case 1; (b) Case 2．等値線は0.1間隔で描画し，u¯1/Ub = 1の等値線のみラベルを付した． CDを直接評価した．具体的には，(i)埋め込み境界法13)に基づいて個々の円柱に働く流体力の時間平均値F¯D を算出，(ii)個々の円柱に対する代表流速urepを，時間平均主流速u¯1を用いて評価，(iii)個々の円柱に対する 抗力係数CDを，2 ¯FD/(ρDu2rep ) により算出した．横嶋・河原3)は，代表流速u repとして，断面平均流速Ub を全ての円柱に対して用い，得られた抗力係数をCD,globalと呼称した．他方で横嶋・河原4)は，代表流速urep として個々の円柱に対する接近流速uaを個別に評価し，求まった抗力係数をCD,localと呼んだ．uaの定義や 値等については3.1節で示す． 本研究では1.節で示した定義に従って各円柱位置でのundisturbed flow uunを直接的に算出する．図3に 示した計算領域内には，Case 1で594本（= 9行×66列）の，Case 2で297本（= 9行×33列）の円柱が含 まれる．以下では，図3において上からi番目，左からj番目の円柱を‘ricj’と呼称する．N本の円柱位置で undisturbed flowを直接評価するには，Nケースの数値実験（注目する円柱が除去され，それ以外の(N-1)本 の円柱が残った系に対する）が必要となることに注意されたい．横嶋・河原3,4)と同様に，本研究でも式(1)は 時間平均流れのレベルで成立するものとする．このとき，円柱配置の主流方向に関する周期性（Case 1）や流 路中央に関する対称性（Case 1および2）のため，統計的に独立な円柱の数はCase 1では5本（= 5行×1列） に，Case 2では165本（= 5行×33列）となる．本研究では，Case 1に関して，これら5本全てを個別に除去
した系を対象とした数値実験を行い，undisturbed flowを直接評価した．他方でCase 2に関しては，165ケー
ス全てを対象とすることは難しいため，円柱群領域内でのuunの空間分布の特徴を捉えることを目指す．すな
わち，流路横断x3方向には3断面: r1, r3, r5を，主流x1方向については12断面：c1, c2, c3, c4, c5, c9, c17,
c25, c30, c31, c32, c33に注目し，計36（= 3× 12）ケースの数値実験を行い，各円柱位置でのundisturbed
flowを直接的に算出した．ある円柱位置でのundisturbed flow uunには，主流速u1の時空間平均値を与えた．
すなわち，当該円柱内部に配置されるu1の空間平均値（重みは全て等しい）を時間方向にも平均することで，
tU
b=D
4000 8000 12000u
un=U
b -0.50 0.5 -0.50 0.5 -0.50 0.5 -0.50 0.5 -0.50 0.5r1c1
r2c1
r3c1
r4c1
r5c1
(a)
T
sampleU
b=D
2000 4000 6000 8000u
u n=
U
b 0 0.5 1.0(b)
図6 Case 1におけるundisturbed flow uunの(a)時系列分布と，(b)その時間平均値のデータ収集時間Tsample に対する依存性．
c
1 9 17 25 33u
un=
U
b;
u
a=
U
b 0.5 1.0 r = 1 r = 2 r = 3 r = 4 r = 5 Circle: Square: Asterisk: Cross: uunin Case 1 uunin Case 2 uain Case 1 uain Case 2 図7 Undisturbed velocity uunおよび接近流速uaの分布．断面平均流速Ubで規格化． られたある瞬時場の情報を与え，12000D/Ubに渡る数値シミュレーションを行い，4000 ≤ tUb/D ≤ 12000 に渡る時間平均操作によって統計量を算出した．図6(a)はそのようにして得られた，Case 1のr1c1, r2c1, r3c1, r4c1, r5c1の位置でのundisturbed flow uunの時系列分布を示す．uunは時間に関して大きく揺らぐこ とがわかる．これは，特に今回は正方格子状に円柱が配置されているために，上流側に隣接する円柱の後流 の影響に因るものと考えられる．(b)では，(a)の時系列データについて，tUb/D = 4000を平均操作開始時刻として2000D/Ub, 4000D/Ub, 6000D/Ub, 8000D/Ub時間に渡る時間平均値をそれぞれ算出し，時間平均
値のデータ収集時間Tsampleへの依存性を表示した．時間平均値は完全に収束しているとは言えないものの，
Tsample ≥ 4000D/Ubではそれぞれの円柱位置でのuunの時間平均値の大小関係は不変であり，本研究の目的 には供すると判断した．ただし，例えば，Bagchi and Balachandar14)は粒子位置でのundisturbed flowとし
て，粒子中心位置での流速に加えて，粒子直径の2倍および10倍の領域に渡る空間平均値についても検討し ている．本研究においてもそういった検討はなされるべきであり，今後の検討課題としたい． 3. 計算結果と考察 3.1 Undisturbed flowと接近流速 得られたundisturbed flow uunの分布を，図7に示す．横軸cの数字は図3において左から何番目の円柱位 置（での代表流速）かを表す．Case 1においては，uunは円柱群側端r1において約0.75Ubの値をとり，r2で 約0.5Ubに急減したのち，さらに円柱群内部に向かうと緩やかに減少する．Case 2のr1断面では，円柱群上
x1[cm] 0 40 80 120 160 200 7u1 = U 0.0 0.5 (a) Case 1. x1[cm] 0 40 80 120 160 200 7u1 = Ub 0.0 0.5 1.0 (b) Case 2. 図8 どの円柱も除去しない場合の数値実験から求まった時間平均主流速u¯1の，各円柱中心を通る主流方向の 5断面（r1, r2,...,r5）での分布. 黒，赤，青，マゼンタ，シアン色の曲線はそれぞれr1(x3 = 52 cm), r2, r3, r4, r5(x3 = 40 cm)断面での時間平均主流速u¯1を表す．黒の点線は各円柱の中心位置を示す． 流端c1から流下方向に進むにつれて，uunはまず増加し，c3位置で最大値およそ1.1Ubをとった後，0.7前後 の値（これはCase 1のr1におけるuunの値に近い）に漸近する様子が見られた．この円柱群上流端近傍での uunの増加は，円柱群パッチを迂回する加速流の存在によると解釈できる．Case 2のr3およびr5断面では， 流下するにつれてuunの値は概ね単調に減少する傾向にあるものの，c17断面以降で両者は異なる分布を示す． r5断面でのuunはc17以降も減少傾向を示し，c31付近で最小値およそ0.3Ubをとった後に下流端c33で0.4Ub まで回復する．下流端での値はr3とr5断面であまり相違せず，また，Case 1のr3およびr5におけるuunの 値にも近い． 横嶋・河原4)で導入された接近流速u aも図7に表示した．図8は，横嶋・河原4)で実施された，どの円柱 も除去しない場合の数値実験から求まった時間平均主流速u¯1の，各円柱の中心を通る主流方向の5断面（r1, r2,...,r5）での分布を示す．主流方向にj番目の円柱への接近流速uaは，図8において(j-1)番目とj番目の円 柱間におけるu¯1の最大値で定義される4)．Undisturbed flow uunは，個別に円柱を除去した系を新たに導入 しない限りは入手できない情報であるため，disturbedな流れ場（どの円柱も除去しない，本来の検討対象流 れ）からuunを推定する方法を確立することは非常に重要である．図7中のuunとuaの分布を比較すると，ua はuunの基本的な特徴は少なくとも定性的には捉えられており，最大誤差は10%程度である．しかしながら， 抗力係数CDは代表流速urepの2乗に反比例するので，urepを10%過小評価するとCDは20%以上も過大評 価される．この点については3.2節で改めて触れる． 3.2抗力係数
図9には，3.1節で得られたundisturbed flow uunを用いて，各円柱（Case 1では5本，Case 2では36 本，図3中の赤く塗られた円柱）の抗力係数CDと円柱レイノルズ数ReD関係をプロットした．2.3節の手順 (iii)に従ってCDを評価する際に用いる時間平均抗力F¯Dは，横嶋・河原3,4)と同じく，どの円柱も除去しな い本来の検討対象流れの数値シミュレーションから得られたもの（文献4)の図–6を参照）である．Finnemore and Franzini15)の教科書に掲載の，孤立円柱に対するC D–ReD関係，2次元流体解析から得られた孤立円柱 のCD–ReDデータ16,17)も併せて示す．⃝は本研究と同一手法を用いた孤立円柱周り流れの2次元数値解析 から得られた結果で，円柱中心を原点として主流x1方向に[−20D, 40D]，横断x3方向に[−20D, 20D]の領域 を435× 800の計算セルに分割（円柱近傍では∆x1= ∆x3 = D/20で等間隔）し，上流端で一様流入(Ub, 0) を，下流端で対流流出条件∂u/∂t + Ub∂u/∂x1 = 0を，横断方向には周期条件を課したものである．図9よ り，円柱群中の個々の円柱に対するCD–ReD関係は，たとえ代表流速にundisturbed flow uunを用いた場合 にも，孤立円柱に対するCD-ReD関係（本研究では，円柱群を過ぎる流れを平面2次元モデルで扱ったので，
比較対象には2次元流体解析から得られた孤立円柱周り流れを使用）から乖離することが読み取れる．すなわ
ち，円柱群中の個々の円柱に，孤立円柱に対するCD–ReD関係を当てはめることは適切でない．他方で，円柱 群パッチを迂回する加速流の影響を強く受けるr1c3や，最もuunが低下するr5c25やr5c32といった一部を除 けば，多くの円柱のCD,unは2.0± 0.5程度の値をとることも判明した．
ReD 101 102 103 104 CD 1 2 3 4
図9 Undisturbed velocity uunに基づくCD–ReD関係．実線–孤立円柱周り流れ(Finnemore and Franzini15)); 白抜きシンボル– 2次元流体解析に基づく孤立円柱周り流れ(△, Park et al.16)-境界適合格子;▽,樽川・ 平野17)- 境界適合格子;⃝,本研究 -埋め込み境界法);中塗りシンボル– 円柱群を過ぎる流れ（本研究， 各シンボルの意味は図7と一致）． 最後に，undisturbed flow uunおよび接近流速uaを代表流速として求めた抗力係数の円柱群領域内での空 間分布を図10に示す．比較のため，全ての円柱に対して断面平均流速Ubを代表流速とした場合の抗力係数3) （CD,global ≡ 2 ¯FD/ ( ρDUb2)）の分布も並示した．Case 2におけるuunは36点でのみ直接評価されたので，そ れ以外の円柱位置でのuunは周囲からの線形補間によって評価し，それを利用して全ての円柱位置でCDを算 出した．接近流速uaは，注目する円柱の表面での粘着・不透過境界条件の影響によって，当該円柱が静止状 態にあるなら，uunよりも一般に小さな値をとること，およびこれらの流速はいずれも断面平均流速Ubより は一般に小さいと予想される．その結果，CD,global < CD,un ≡ 2 ¯FD/(ρDu2un
) < CD,local ≡ 2 ¯FD/(ρDu2a ) の関 係がよく成立すると期待される．実際に，図7および図10から，上述の予想が概ね正しいことがわかる． 横嶋・河原 4)では，C D,global ≡ 2 ¯FD/(ρDUb)は円柱群抵抗を過小評価する傾向がみられること，および CD,localは円柱群抵抗を過大評価する傾向にあること，が示された．既に述べたように，接近流速u¯aは注目す る円柱の周囲の円柱の影響だけでなく，当該円柱自体の影響も受けた流速である．他方で，断面平均流速Ub はどの円柱の影響も受けない流速と捉えられる．注目する円柱以外の全ての円柱の影響を受けるundisturbed flow uunはこの両者の間に位置すると解釈でき，円柱群抵抗の予測に関する上述の問題を克服することが期待 できる． 4. おわりに 本研究では，多数の障害物を過ぎる流れを巨視的抗力モデル（式(1)）で表現する場合に必要な抗力係数に 関して，特にその代表速度として，undisturbed flow（注目する障害物が存在せず，それ以外の障害物は全て 存在する系において当該障害物位置を過ぎる流速）を用いることの是非の検討を試みた．具体的には，円柱群 （樹木群模型）を過ぎる直線開水路流れ（室内水理実験）を水平2次元モデルで近似し，個々の円柱周りの流 動および各円柱に作用する流体力を支配方程式に基づいて算出する数値実験を実施した．個々の円柱に対する undisturbed flowを，実際に当該円柱を除去した系をシミュレートすることで直接的に評価し，その基本特性 を検討した．
円柱群中の個々の円柱に対するCD–ReD関係が，代表流速にundistured flowを用いた場合に，孤立円柱
のCD–ReD関係とどの程度一致するのかを検討した．両者には明瞭な差異が認められ，たとえ代表流速に
undistured flowを用いた場合にも，孤立円柱のCD–ReD関係を円柱群中の個々の円柱に適用することは適切
でないことが確認できた．他方で円柱群中の多くの円柱のCDは，代表流速にundistured flowを用いた場合，
2.0± 0.5程度の値をとることもわかった．
著者らの先行研究3,4)では，代表流速に断面平均流速U
x1[cm] 0 40 80 120 160 x3 [c m ] 0 40
0:34 5 C
D;global5
1:06;
hC
D;globali = 0:59
0.4 0.8 1.2 1.6 2.0 x1[cm] 0 40 80 120 160 x3 [c m ] 0 400:34 5 C
D;global5
2:12;
hC
D;globali = 0:92
0.4 0.8 1.2 1.6 2.0 x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (c)2:17 5 C
D;local5
3:23;
hC
D;locali = 2:71
1.0 2.0 3.0 4.0 x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (d)1:79 5 C
D;local5
4:15;
hC
D;locali = 2:46
1.0 2.0 3.0 4.0 x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (e)1:90 5 C
D;un5
3:06;
hC
D;uni = 2:39
1.0 2.0 3.0 4.0 x1[cm] 0 40 80 120 160 x3 [c m ] 0 40 80 (f)1:36 5 C
D;un5
4:13;
hC
D;uni = 2:13
1.0 2.0 3.0 4.0図10 水平2次元モデルの数値実験から得られた，円柱群の抗力係数分布．(a) Case 1, CD,global≡ 2 ¯FD/(ρDUb2);
(b) Case 2, CD,global; (c) Case 1, CD,local ≡ 2 ¯FD/(ρDu2a
)
; (d) Case 2, CD,local; (e) Case 1, CD,un ≡
2 ¯FD/
(
ρDu2un); (f) Case 2. ⟨CD⟩は当該ケースの全ての円柱のCD値の算術平均を表す．CD,globalの値
は他の2種類のCDよりもかなり小さいため，カラーバーの範囲が異なることに留意されたい．また，
見易さのため，円柱は実際のサイズよりも大きく描画した．
を用いて得られたCD,globalを巨視的抗力モデル（式(1)）に適用した場合に，円柱群抵抗が過小評価されるこ と，および，代表流速に接近流速ua（詳細は3.1節参照，全ての円柱の影響を受けた流速）を用いて得られた CD,localを抗力モデルに与えた場合，円柱群抵抗が過大評価される傾向にあることが示された．注目する円柱 以外の円柱の影響を受けるundisturbed flowを代表流速とした場合，そこから得られる抗力係数CD,unは一般 にCD,globalとCD,localの間の値をとることが示された．よってCD,unを巨視的抗力モデルに適用した場合，上
述の過小/過大評価が解消され，より正確なキャノピー流れ解析が実現する可能性が期待できる．この点につ いては今後の検討課題としたい． 謝辞 本稿に対して査読者から有益なご指摘を頂きました．また，本研究の一部は科研費補助金（20K04703）の 支援を受けて行われました．記して謝意を表します． 参考文献 1)例えば，榎木 康太, 石原 孟: 一般化キャノピーモデルの提案と都市域における風況予測への応用, 土木学 会論文集 A1, 68(1), pp.28–47, 2012.
2) Yokojima, S., Kawahara, Y. and Yamamoto, T.: Impacts of Vegetation Configuration on Flow Structure and Resistance in a Rectangular Open Channel, J. Hydro-environ. Res., 9, pp.295–303, 2015.
3)横嶋 哲,河原 能久: キャノピー内部の抗力係数分布を考慮した樹木群を過ぎる流れのLES,土木学会論文 集B1, 71(4), pp.I 1051–I 1056, 2015.
4)横嶋 哲, 河原 能久: 個々のキャノピー要素への接近流速を考慮した抗力係数分布に基づく樹木群流れの
LES, 土木学会論文集A2, 71(2), pp.I 703–I 711, 2015.
5) Maxey, M.R. and Riley, J.J.: Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids, 26, pp.883–889, 1983.
6)例えば，横嶋 哲,島田 佳昭: 微小粒子の乱流中のクラスタ特性に対するバセット履歴力の影響, 土木学会 論文集 B1, 75(2), pp.I 595–I 600, 2019.
7) Horwitz, J.A.K. and Mani, A.: Accurate calculation of Stokes drag for point-particle tracking in two-way coupled flows, J. Comput. Phys., 318, pp.85–109, 2016.
8) Ireland, P.J. and Desjardins, O.: Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling, J. Comput. Phys., 338, pp.405–430, 2017.
9) Akiki, G., Jackson, T.L. and Balachandar, S.: Pairwise interaction extended point-particle model for a random array of monodisperse spheres, J. Fluid Mech., 813, pp.882–928, 2017.
10) Eaton, J.K.: Two-way coupled turbulence simulations of gas-particle flows using point-particle tracking, Int. J. Multiphase Flows, 35, pp.792–800, 2009.
11) You, J., Choi, H. and Yoo, J.Y.: A modified fractional step method of keeping a constant mass flow rate in fully developed channel and pipe flows, KSME Int. J., 14(5), pp.547–552, 2000.
12)例えば，Williamson, C.H.K.: Vortex dynamics in the cylinder wake, Annu. Rev. Fluid Mech., 28, pp.477– 539, 1996.
13) Uhlmann, M.: An immersed boundary method with direct forcing for the simulation of particulate flows, J. Comput. Phys., 209, pp.448–476, 2005.
14) Bagchi, P. and Balachandar, S.: Effects of turbulence on the drag and lift of a particle, Phys. Fluids, 15, pp.3496–3513, 2003.
15) Finnemore, E.J. and Franzini, J.B.: Fluid Mechanics with Engineering Applications (10th edition), McGraw-Hill, New York, 2001.
16) Park, J., Kwon, K. and Choi, H.: Numerical solutions of flow past a circular cylinder at Reynolds numbers up to 160, KSME Int. J., 12, pp.1200–1205, 1998.
17)樽川 智一,平野 廣和: 広範囲なReynolds数域における円柱まわり流れの数値流体解析, 第20回風工学シ ンポジウム論文集, pp.66–71, 2008.
