第 26 回数値流体力学シンポジウム
講演番号 B11-2

VOIS 法計算を応用した流れの整流機能とその効果
Function of Straightening Flow by Using Calculation Procedure of VOIS Methodand Effect of the Function○ 齋藤英文, ㈱島津製作所, 京都市中京区西ノ京桑原町１番地, E-mail:h-saito@shimadzu.co.jp
Hidefumi SAITO, Shimadzu Corporation, Nishinokyo-Kuwabara-cho 1, Nakagyo-ku, Kyoto, Japan
Velocity-field Overridden by Immersed Solid (VOIS) Method is a branch of Immersed Boundary (IB) Method anduseful for calculating fluid field including complicated solid parts. The VOIS method is effective not only to graspflow fields including solids but also to create virtual flow field that is suitable for estimation. I tried to apply the VOIS
method to making straightened flow field as one of the particular fields. To verify usefulness of the straightened fieldby the VOIS method, drag coefficients of objects in the field were evaluated as an indicator of effect. Since dragcoefficient of sphere became considerably close to correct value, the straightened field was enough useful for roughestimation of drag coefficients. However, even though thickness of straightened layers was multiplied ten times, all of
the straightened flow velocities were not precisely conformed to set value.

１．はじめに流れの中に置かれた物体（たとえば流体部分に突出した部材や流体中に混入した固体物質）に対して，流体から物体に作用する力や，物体によって現れる流れの変化について，簡便な計算法で定量的な概算値を得ることができれば，設計を進める際の判断材料を得るという点で有用となる．従来から，流体内部に存在する固体を簡便な直交格子で計算し，流れの状態を知る方法として，埋め込み境界法(1)（Immersed Boundary 法，以下 IB 法と略す）が知られている．IB 法は，固体が含まれる流れ場に対し，固体の形状に捉われずに直交格子を設定し，全格子を一旦流体として計算した後，各格子内おける固体の存在割合に基づき流速を強制し，流れ場を解く方法である．IB 法の中でも，速度強制を行うための計算で，流体と固体が混在する格子に対して，運動量の加重平均した流速に上書きする VOIS 法(2)は，計算負荷を少なくできる優れた方法である．
元来 VOIS 法は，多数の固体粒子が含まれた流れ場を直接計算するため考案され，固液二相流の解明(3)に活用された．その後，VOIS 法の持つ特長を活かし，柔構造の固体と流体との相互作用の解明(4)等にも応用され，現象をうまく捉えた結果が得られている．従来の研究では，VOIS 法を利用して自然現象に含まれる複雑な状態を，適切に表現する計算モデル構築に重点が置かれた．
VOIS 法の持つ計算負荷低減効果は，設計分野で利用の可能性がある．例えば，設計対象の寸法や形状を決める際に，特定の影響因子のみを抽出し，設計対象が流体から受ける影響を調べるような場合である．
流れ場の計算で，複数の影響因子が重畳する場合，特定因子の強調や排除を行った際，流れの中に存在する対象物体が受ける力の差を比較するのが一例である．このような流れ場の設定にはいろいろな方法があるが，VOIS 法の場合，変化を与える特定の格子に対し，目標速度を上書きする操作を加えることで，異なった流れの状態を簡便に与えることができるという特長を有する．
筆者はこの点に着目し，これまでにも空間内に仮想吹出しノズ
ルを設定し，最適な吹出し形態を探索する方法において，目的とする流れ場を構築するために，仮想吹出しを移動・評価する手法を提案した(5)．今回の研究では，空間内の流れ場で，外乱を排除し整流した流れ場を形成することを試みた．その結果，VOIS 法による整流では，整流の有無によって抗力の違いを得ることが可能となった．しかし，圧力分布が残るため，厳密性を求める整流への適用は難しいことが判明した．以下，詳細に報告する．
２．評価する対象
簡便な計算を目指すために，PC レベルで計算できる事例として，Fig.1 に示す，一様な流れ場内に物体が存在する計算空間を題材にした．
xObject（Solid）

wm

dm
z

2.5 dm

6.25 dm

7.5 dm
Fig.1 Domain of Calculation

Fig.2 Drag coefficient for spheresas a function of the Reynolds number (6)1Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-2
度 Adams-Bashforth 法を用いた．また，流れ場と圧力場の連立解を求める圧力方程式の計算には，過緩和係数β=1.3 の逐次過緩和法（SOR 法）を用いた．

物体は，流れ方向を軸とした回転軸対象の固体とし，断面積最大となる直径を dm，とする．最大断面積を Am (=πdm2/4)，一様流の流速を wm，物体に作用する抗力が F のとき，抗力係数 Cd は，次式で与えられる．

Cd =F.
ρ wm 2 Am 2(2)．境界条件
計算空間の境界条件は，流れ状態が前後左右に伝達されることを基本とするため，６面に対し周期境界条件を与えた．ただし，z方向の入口面では整流の操作を行う格子を設定した．Fig.1 の中で整流操作を行う部分を赤色に，他の周期境界条件面を緑色に表す．
赤色部分の対象格子に対して，VOIS 法に従い，速度を次の値に強制した．

(1)

物体が球体の場合には，抗力係数が Fig.2 に示される値となることが，古くからよく知られている．
３．支配方程式
計算の支配方程式は，非圧縮性流体の方程式を用いる．

∇⋅ u = 0 ,(2)∂u1= −∇ ⋅ (uu) − ∇P + ∇ ⋅ (2ν D) ,
ρ∂t

(3)

w = wm ,u = v = 0,(9)
流体内に存在する物体の周辺格子では，流体：固体の割合に基づき，VOIS 法に従って，速度強制を行った．
５．流体により物体に作用する力の計算(1)．物体の表面への作用力の計算
流れ方向を軸とした回転対称な物体に作用する力は，変動分を除くと，軸に対し回転対称となる．そこで，Fig.3 に示すように，格子幅で輪切りにした部分表面を設定した．切り出した部分表面
では，変動分を時間平均することで，同一周上では，部分表面が流体から受ける応力の大きさは一定の値になり，応力のｚ方向成分も一定となる．物体への作用力を求める手順として，部分表面上に対応する各位置で作用力のｚ方向成分を抽出し，一定時間と周方向とで平均化した値を求め，部分表面の x-y 断面への投影面積を乗じ，当該の部分表面のｚ方向作用力とした．各部分表面のｚ方向作用力を，物体の先端から後端まで総和を取り，物体に作用する力とした．

ここに，u ：流速ベクトル，
D ：ひずみ速度テンソル，

P：圧力，
ρ：密度，
ν：動粘性係数．
本モデルでは，物体の後側に現れる渦の動的変化を捉えるとともに，この渦を整流した際の状態を知るために，LES による計算を採用した．格子以下のスケール（Sub-Grid Scale，以下 SGS と略す）は，van Driest 型減衰関数を採用したスマゴリンスキーモデルとする．

ν = ν 0 +ν SGS ,(4)

ν SGS = ( Cs f s Δ ) D ,(5) − y+  ,f s = 1 − exp  25 (6)2(8)

Its z-component of the stressvector on partial surface

yy+ =luτ ,z(7)ν0x

ここに，ν0：流体の分子拡散による動粘性係数，

νSGS：SGS 渦粘性，()Δ = 3 ΔxΔyΔz ：フィルタ幅，Solid object
StreamCS ( = 0.1) ：スマゴリンスキー定数,

Ring-like partial surface

D ：ひずみ速度テンソルの大きさ，Fig.3 Partial surface and its stress vector

f S ：減衰関数，
l ：固体からの最短距離．
壁からの無次元距離 y + の算定には，一様流れに平行な壁面で(2)．抗力の圧力寄与分
物体に発生する抗力は，圧力によるものと流体剪断摩擦によるものとの和となる．まず，圧力の寄与分の計算方法を示す．
圧力値は，各格子中心に対して設定されており，物体表面では定義されていない．そこで，格子中心が物体の表面に近い位置にあるほど，圧力値を大きく重み付けする計算方法を検討した．
Fig.3 に示す部分表面毎に，流体の割合が０を超え 1 未満となる格子（すなわち内部に物体表面の一部が含まれる格子）の圧力値を抽出し，格子内に存在する物体表面の面積に応じて，圧力値を加重平均した．この操作で得られた平均圧力が，当該の円環状の部分表面に均一に作用しているものと仮定して，当該の部分表面への作用力を求めた．

境界層を発達させた際の，壁面直近（1/2 格子位置）に現れる壁面近傍速度から得られる摩擦速度 uτ を全計算空間に適用する方法を用いる．
４．計算スキームと境界条件(1)．計算スキーム
本研究では，Fig.1 に示す空間を x×y×ｚの各方向について100×100×160 の Staggered 格子に区切った流れ場の計算は，空間差分に 2 次精度中心差分を，移流項・粘性項の時間進行は SMAC 法を，時間進行予測段階の計算は 2 次精
2Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-2
領域を対象としており，②の発生領域には該当せず，本研究では影響しないと考える．
以上の検討から，剪断摩擦による抗力の計算には課題が残されているものの，本研究における抗力の計算は，有効なものと判断するに至った．

(3)．流体摩擦によって生じる抗力
剪断摩擦による応力は，物体表面近傍の，表面に垂直な方向に現れる速度勾配に粘性係数を乗じて得られる．
流速は格子境界面に離散的に与えられているため，通常，物体表面での正確な速度勾配を得るには，十分な格子解像度が必要になる．しかし，まず今回は概算値を得ることを狙い，粗い格子の流速値から，物体に作用する剪断応力を求めることを試みた．速度が設定される位置が，物体表面の外側に存在し，かつ物体表面からの距離が格子幅より小さくなり，その速度は強制されていないものを抽出した．物体表面からの距離と抽出した速度から速度勾配を求め，剪断応力野のｚ方向成分の平均値を，Fig.3 に示す部分表面毎に計算した．

８．VOIS 法を活用した整流による効果
本研究での課題である VOIS 法を活用した整流による効果について評価した．一様流の流速を 0.25m/s とし，流体内の物体を球体とした場合で，整流操作を行わない場合（整流操作の対象格子が 0 層の場合）と，整流操作の対象格子を１層，5 層，10 層とした各場合で得られた球体の抗力係数を Fig.4 に示す．この図から明らかな通り，整流が行われない時には，抗力係数が 0.27 以下となっていたものが，整流が行われることで 0.36 以上の値となり，速度強制を行う層を厚くすることで，Fig.2 の Re=9.5×103 における抗力係数（＝0.40）に近づく結果が得られた．設計判断のため，例えば抗力係数に 10%程度の誤差が許容される場合には，1 層の速度強制であっても，判断に必要な情報が得られる．
一方，速度強制を 10 層に増加しても，Fig.2 に示される 0.40 に達することができなかった．VOIS 法を活用した整流の効果は，1層でも効果が得られる反面，整流対象を 10 層まで拡大しても，厳密な整流状態を形成することが難しいという結果となった．

６．設定条件Fig.2 に基づき，抗力係数が Re 数に対しほぼ一定となる Re=103〜105 を狙って評価することを試みた．物体の寸法は，dm=40mmを採用し，流体の物性値を Table.1 のとおり設定した．
一様流の流速 0.25m/s を基本とし，基本流速より遅い流れ0.10m/s と，速い流れ 0.50m/s の３種類に流速で計算を行う．表中に，dm を代表寸法として算出した Re 数を示す．
Table 1 Physical properties of fluid
1.075×10-3viscosity coefficient μ(Pa s)
1.049×10-6
kinematic viscosity ν(m2/s)
1.024×103
density ρ(kg/m3)
maximum diameter dm (m)0.040
stream velocity wm (m/s)
0.10
0.25
Re-number
3810
9526

0.50
19051

７．抗力計算の有効性
文献(6)(7)によれば，境界層厚さδは，δ≈dm ,Re(10)のレベルとなるので，本研究の Re=3.8×103〜1.9×104 では，計算格子間隔で境界層を捕えるには十分な解像度が得られない．そこ
で，５項に示す抗力計算についての有効性を検討した．
境界層が格子に対して薄く，格子解像度が不十分となった場合
でも，境界層の外の流れ場が正しく捉えられていれば，物体周りに現れる圧力分布は正しく評価される．しかし，①境界層内の流速で決まる物体直近での速度勾配は正しく評価されず，物体表面の剪断応力が正しく計算されない，②境界層が乱流遷移した場合，境界層の剥離点が移動するが，境界層が正しく捉えられていないと，境界層の乱流遷移を検出することができず，剥離点の位置や近傍の圧力分布の変化が正しく評価できない，という問題を生じる．①の問題について，今回の計算で得られた抗力の剪断摩擦による寄与分を調べた．その結果，例えば一様流速が 0.25m/s の場
合，圧力による寄与分 1.43×10-2 (N)に対し，剪断摩擦による寄与分は 1.23×10-4 (N)となり，圧力による寄与分と比較すると，剪断摩擦による寄与分は２〜３桁小さな値となった．実際には境界層が薄く，速度勾配が計算結果より数倍大きくなった場合でも，圧力による寄与分が支配的なため，抗力全体はせいぜい 2〜3%の違いとなって現れる程度である．②については，急激な抗力係数低下が現れる所が，境界層の乱流遷移に伴う剥離点移動の発生に該当する．Fig.2 から分かるとおり，急激に抗力係数が低下するのはRe=3×105 付近となっている．今回の計算は Re 数が１桁以上低い

Fig.4 Layers for velocity forcing vs. drag coefficientFig.4 に示された結果が得られる理由を明らかにするために，速度強制を行わない場合，速度強制を 1 層で行った場合，速度強制を 10 層で行った場合のそれぞれについて，物体の回転軸を通る面に現れた流れ場の様子（Fig.5〜7）を観察した．図では u，w をベクトル，v をコンタで示す．速度強制が 1 層の場合，速度強制層により，物体の後方に発生し物体前方に回り込む渦のほとんどが除去されているが，流れ方向に筋状となる速度ムラ（周囲と異なる流速の領域）が残った．速度強制の範囲を 1 層から 10 層に拡大しても，速度強制層を通過した流れには，筋状の速度ムラが確認され，完全な整流には至らなかった．
物体の上流側の流れに含まれる渦等の乱れが顕著なほど，物体の前面にて圧力の低下する箇所が発生し，抗力係数が低下し易くなる．Fig.6〜7 と比較すると，Fig.5 では明らかに物体の上流側の乱れが大きく，抗力係数の低下に至ったと判断できる．Fig.6〜7は完全な整流状態となっておらず，物体前面に圧力が低下する箇所が残れば，Fig.2 に示される抗力係数に到達できない理由となる．
以上，Fig.5〜7 で得られる知見は，Fig.4 の結果と符合する．VOIS法の速度強制による整流は，簡便な評価に供される流れ場では利用可能であるが，厳密な整流状態の実現は難しいと言える．
3Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-2
速度強制での整流層を通過した流れの中に，速度ムラが残る理由について考察した．VOIS 法に基づく速度強制は，各時間ステップにおいて，一旦計算空間全体を流体とみなし，Navier-Stokes方程式による解を得た後，強制する速度値に上書きする．上書きされた後の速度分布を基に，次の時間ステップの計算に進む．
これらの段階における，整流の対象となる格子群に現れる流れ場のイメージを Fig.8 に示す．図では，圧力分布はコンタで示される．(a)は Navier-Stokes 方程式による解が得られた状態を示す．
この時点で流速と圧力の分布は，方程式を満たす状態となっている．(b)は，(a)で得られた速度に対し，整流する流速値を上書きした後の様子を示す．この時の圧力分布は，(a)で得られた状態が残されている．(c)は，(b)の状態から次の時間ステップとして，1 ス
テップ分の時間経過に伴う移流分や拡散分を加え，Navier-Stokes方程式による解を再度得た状態となっている．(c)では，流速と圧力は再配置される．

Contour: Pressure
Arrow: Velocity
Fig.5 Flow pattern with non-layer for velocity forcing(a) Before treatment by straightening

Fig.6 Flow pattern with 1-layer for velocity forcing(b) Overriding velocities by VOIS method

Red arrow: AcceleratingBlue arrow: Deceleratingor reversed flow(c) Resolving Navier-Stokes Equation
Fig.7 Flow pattern with 10-layers for velocity forcing
Fig.8 Shift image of velocity and pressure distribution
in straightened area based on the VOIS method
4Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-2
する抗力係数の計算を試みた．球形の対象物体に対する抗力係数を求めたところ，1 層の整流でも，実際の値(0.40)に近い0.36 を得ることができた．また，整流機能を止めた場合には，物体前方から流れてくる乱れの影響により，抗力係数は 0.27に低下した．

その後，(b)と(c)が繰り返される．この繰り返しを進める過程
で，整流対象域では，流速と圧力が均一に分布する状態に近づく．
しかし，整流対象域内の流れは，周囲の計算領域に存在する圧力の分布や，上流側からの乱れを含む流れの流入によって，常に影響を受けており，乱れを完全に消去した流れ状態を形成し維持するには，現行の VOIS 法の計算操作に対し，何らかの改良を加える必要があると考える．

・９．抗力係数 Cd 値の計算結果例
VOIS 法を活用した整流を用い，異なる流速や異なる物体形状で，抗力係数を求めることを試みた．抗力係数は概算値で良いと仮定し，整流を行う層は 1 層のみとした．一様流の流速が基本流速（0.25m/s），遅い流れ（0.10m/s），速い流れ（0.50m/s）のそれぞれについて，球形に加え，z 方向の軸長が短い形状（=0.5dm）と長い形状（=1.5dm）の物体で，それぞれ抗力係数の計算を行った．
結果を Fig.9 に示す．図中には，物体形状のイメージ図も併せて示す．
Fig.9 から，① 計算した流速範囲（Re=3.8×103〜1.9×104）であれば，流速が変化しても抗力係数はほぼ一定の値となる，② z方向軸長の，最大断面直径 dm に対する軸比が大きくなるに従い，抗力係数は低下し，軸比が 0.5→1.0，1.0→1.5 のそれぞれの変化で，抗力係数は約半減する，という結果を確認した．

厳密性が求められる整流された流れ場を形成するには，VOIS法の計算手法に中に，圧力場も強制すると言ったような改善を加えて行く必要がある．この点は，今後の課題としたい．
今回の抗力係数の計算では，対象物体表面に現れる境界層内の流速勾配を正しく評価する点に課題が残った．より広い範囲の Re数の流れに対し抗力係数を評価して行くには，この課題解決は必要であり，さらに検討を加えたい．
謝辞
本研究を遂行するに際して，大阪大学大学院 工学研究科 機械工学専攻の梶島岳夫教授，竹内伸太郎准教授から貴重なご助言を賜りました．ここに深く謝意を表します．
参考文献(1) Mittal, R. and Iaccarino, R., “Immersed Boundary Methods”,
Annual Reviews Fluid Mechanics, Vol.37, (2003), pp.239-261.
(2) 梶島，竹内, 「固気液混相流の直接数値シミュレーション」,
機械学会論文集 B, Vol.76(765), (2010), pp. 712-719.
(3) Kajishima, T. and Takiguchi, S., “Interaction Between Particle
Clusters and Fluid Turbulence”, Int. J. Heat and Fluid Flow, Vol.23,
Issue 5, (2002), pp.639-646.
(4) Takeuchi, S., Yuki, Y., Ueyama, A., and Kajishima, T., “Interaction
Between Multiple Deformable Particless and Fluid Flow by a
Coupled Finite-Difference and Finite-Element Technique with aConservative Momentum-Exchange Algorithm”, Int. J. NumericalMethods in Fluids, Vol.63, Issue 5, (2010), pp.1084-1101.
(5) Saito, H., and Kajishima, T., “Searching Optimal Jet Point in
Aircraft Cabin Flow Field Coupled with Thermal Field”, Proc. 3rd
Asian Symposium on Computational Heat Transfer and Fluid Flow,(2011), #124.
(6) Schlichting, H., “Boundary-Layer Theory (Seventh Edition)”,McGraw-Hill, (1979), pp17 and Chapter VII.
(7) 吉澤, 「流体力学」, 東京大学出版会 (2001), 4.4 節．

Fig.9 Calculated drag coefficient
１０．結論
IB 法の一手法である VOIS 法を活用し，流体の計算空間内部に整流機能を付与し，整流された流れ状態について評価を行った．
整流状態の評価として，流れ方向に軸を持つ回転体の抗力係数を指標として求め，その値について考察した．本研究により，次の結論を得た．
・VOIS 法による速度強制の手法を活用した流れの整流操作は，1 層のみの整流であっても，設計に必要な係数等の概算値を得るための流れ場を形成することが可能となる．VOIS 法の速度上書き操作は，このような特異な流れ場を形成するために有効な方法である．

・

厳密な整流を目指し，速度強制を行う層を 10 層に増加し，VOIS 法による流れの整流を試みた．しかし，整流後の空間には筋状の流速ムラが残り，厳密な均一流速空間を形成することは困難であった．

・VOIS 法の速度強制手法による整流化した流れ場を用い，設計に必要な係数の概算する例として，回転対称の物体に作用対象物体の形状や，流れ場の流速を変え，抗力係数を計算した結果，Re が 3.8×103〜1.9×104 となる範囲では，同形状で流速が変化しても抗力係数はほぼ一定に保持されること，物体形状の軸比（z 方向軸長／dm）の増加に応じ，抗力係数が低下する結果とを得ることができた．

5Copyright © 2012 by JSFM