第 25 回数値流体力学シンポジウム
C03-1

デジタルカラー画像を用いた非圧縮性流れに対する仮想領域法
A domain embedding method for incompressible viscous ﬂow using digital color image
吉井孝子, 千葉大院, 〒 263-8522 千葉市稲毛区弥生町 1-33, E-mail：takako-yoshii@graduate.chiba-u.jp

周曼, 千葉大院, 〒 263-8522 千葉市稲毛区弥生町 1-33, E-mail：zhouman@graduate.chiba-u.jp

腰越秀之, 千葉大院, 〒 263-8522 千葉市稲毛区弥生町 1-33, E-mail：koshigoe@tu.chiba-u.ac.jp

Takako Yoshii, Department of Urban Environment Systems,Chiba University, Chiba 263-8522,JapanMan Zhou, Department of Urban Environment Systems,Chiba University, Chiba 263-8522,JapanHideyuki Koshigoe, Department of Urban Environment Systems,Chiba University, Chiba 263-8522,Japan
In this article we propose a numerical method to analize a incompressible viscous ﬂow modelled by the NavierStokes equations around circular cylinders and buildings in Lisbon. The solution method combines ﬁnite diﬀerence
method, image processing and ﬁctitious domain method. The image processing we used here is very usefull toestablish a ﬁnite diﬀerence method and to construct a regular cell for the complex shaped domain since thedigital color image consists of so many picxels and each picxel forms a small square, is addressable and has acolor intensities. Developing these techniques makes it possible to compose a new numerical algorithm of ﬁnite
diﬀerence method and to apply to various problems in enviornmental mathematics.

1.

序

本稿では，デジタルカラー画像を用いた障害物（円柱）周りの流れの数値計算法を示し，その応用として，障害物がリスボン市街地の建物の場合の風の流れの数値シミュレーションについて考察する。ここでの数値計算の特徴
はデジタルカラー画像と有限差分法の結合である。それは，デジタル画像は極微小な正方形のピクセルから構成され，各ピクセルは，赤 (R), 緑 (G), 青 (B) の輝度値，および，アドレスを有しているので，その情報を有限差分法に取り込むことである。この方法は，どのような領域形状でも数値計算ができるし，いろいろな環境数理問題にも広く応用が可能である。

Γ1 から流れを与えた時，Ω1 での Navier-Stokes
方程式：∂v
1+ (v · ∇)v = −∇p +Δv in (0, T ) × Ω1
∂tRe∇·v =0in (0, T ) × Ω1を考察する。境界条件は (a-1)，(a-2) である。

(a-1) 円柱 Ω2 , Ω3 , Ω4 の側面 γ 上の境界条件：∂p=0∂nv=0(1) 本稿で用いる画像は Fig.1 および Fig.2 である。

on (0, T ) × γ(1)

on (0, T ) × γ(2)(a-2) 外周 Γi (i = 1, 2, 3, 4) 上の境界条件:
v = (0.2, 0.0)

on (0, T ) × Γ1

∂v=0∂np = 0.0on (0, T ) × (∂Ω1 \(Γ1 ∪ γ))∂p=0∂nFig. 1: Digital color image of 3 circular cylinderson (0, T ) × Γ2on (0, T ) × (∂Ω1 \(Γ2 ∪ γ))(b) Fig.2(Google Earth) と関連資料から, 横 1000ピクセル，縦 300 ピクセルのカラー画像を作成する。

Fig. 2: Lisbon and the Tejo river (Google Earth)Fig. 3: Digital image of Lisbon city
(2) 上の画像を用いて，次の (a), (b) の数値計算を行う。
Ω1 は白い領域で道路と広場 (A,B,C) を表し，Ω2は黒い領域で建物 (障害物) を表す。また，Ω3 は青い領域でテージョ川を表す。

(a) 3 つの円柱を含む画像 Fig.1 では，白い領域を
Ω1 , 3 つの円柱をそれぞれ Ω2 , Ω3 , Ω4 ，および外周の境界 Γ1 , Γ2 , Γ3 , Γ4 とする。
1c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
C03-1

ここでは，テージョ川から吹く風を与えて，次の Navier-Stokes 方程式：に対して、次の方程式系を設定する。なお，(13)-(14) がγ 上の接合境界条件である。

−Δu1 = fin D1
in D2−Δu2 = 0on γu1 = u2∂u1∂u2= on γ∂n∂non Γ10u1 = 0∂u1= 0on Γ11∂n∂u
この時， → 0 時， 1 → 0 on γ が成り立つ。
∂n∂v
1+ (v · ∇)v = −∇p +Δv∂tRein (0, T ) × (Ω1 ∪ Ω3 )∇·v =0in (0, T ) × (Ω1 ∪ Ω3 )を考察する。境界条件は (b-1)，(b-2) および付加条件は (b-3) である。
(b-1) 建物 (障害物)Ω2 の境界面 γ 上の境界条件：∂p= 0 on (0, T ) × γ∂nv = 0 on (0, T ) × γ(3)(4)(b-2) 外周 Γ 上の境界条件:

この状況を次の１次元問題で説明する。

(i) 速度ベクトルの境界条件：
v = (0.2, 0.0)
on (0, T ) × Γ1(5)∂v= 0 on (0, T ) × (Γ2 ∪ Γ3 ∪ Γ4 ) (6)∂n(例） ディリクレ・ノイマン境界問題
d 2 u1(x) = π 2 cos(πx) (0 < x < 1)dx2
u1 (0) = 1
du1(1) = 0dxを考察する。この解は u1 (x) = cos(πx) である。
−(ii) 圧力の境界条件：
p = 0.0on(0, T ) × Γ2(7)∂p= 0 on (0, T ) × (Γ1 ∪ Γ3 ∪ Γ4 ) (8)∂n(14)(15)(16)

(17)
(18)
(19)

この境界値問題に対して，次の方程式系を準備する.

(b-3) Ω2 と Ω3 における付加条件：
v = (0.0, 0.0) in (0, T ) × Ω2
v = (0.2, 0.0) in (0, T ) × Ω3d2 u1(x) = π 2 cos(πx) (0 < x < 1)dx2
d 2 u
− 22 (x) = 0 (1 < x < 2)dxu1 (1) = u2 (1)du1du(1) =  2 (1)dxdxu1 (0) = 1u2 (2) = 0(9)(10)−

ここで, v, p および Re は流体の速度ベクトル, 圧力およびレイノルズ数を表す. n は境界面での外向き単位法線ベクトルである。

2.

(11)
(12)
(13)

ノイマン境界条件と仮想領域法

圧力に関して障害物 (３つの円柱および建物) の境界にノイマン境界条件 (1)(または (3)) が課されている。このノイマン境界条件の解析的近似法として，障害物の領域を計算領域内に埋め込み (domain embedding)，障害物領域には特異摂動の方程式とその境界 γ 上には接合境界条件を設定するという方法がある (4) 。この方法を流体の圧力ポアソン方程式に適用する。まず，この方法の概要を示す。

(20)(21)∂u1(1) → 0 に収束する．
∂x
実際，この解を構成する。未知数 α を導入し，(20) からこのとき， → 0 の時，u1 (1) = u2 (1) = α
とおく．各区間 (0,1) と (1,2) で α をディリクレ境界値として解き，接合境界条件 (21) を満たすように計算すると

γΓ11D1Γ10αD2= −11+1)x + cos(πx)
1+u1 (x)=(1 −u2 (x)=1(x − 2)1+を得る。明らかに， → 0 の時

Fig. 4: Fictitious domainu1 (x) → cos(πx)(= u1 (x))

領域 D1 におけるポアソン方程式：−Δu1 = fin D1∂u1=0on γ∂nu1 = 0on Γ10∂u1=0on Γ11∂ndu1(1) → 0dxに収束する。
下図は，上の方程式系の問題に対する数値解（赤色）を求め， → 0 の時，真の解（青色）に収束するアニメーションの１コマである。

2c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
C03-1

ノイマン境界条件 (3) が課されている。このノイマン境界条件を２節，３節を用いて近似する (1),(2),(3) 。
計算セル上の関数：
1 on (Ω̃1 ∪ Ω̃3 )a(x, y, t) =(22)εon Ω̃2および，Fig. 5: Fictitious domain method3.

f˜(x, y, t) =f (x, y)
0on (Ω̃1 ∪ Ω̃3 )on Ω̃2(23)を定義する。ここで仮想領域とデジタル画像

デジタル画像は障害物領域を画像として自動的に含む。
また，デジタル画像は極微小な正方形をしたピクセルから構成されている。このことは，計算対象領域をデジタル画像にすると，障害物もピクセルで近似できることを示している。Fig.3 の画像は横 1000 ピクセル縦 300 ピクセルから構成されている。この 3 × 105 個の各ピクセルはアドレスと (R,G,B) の輝度値の情報を持つ。我々は横ｘ縦が Nc × Nc ピクセルの集合をセルと定義する。我々が作成したプログラムは 24 bitBMP ファイルの画像からピクセルのアドレスと (R,G,B) の輝度値の情報を読み取り、セルがどの部分領域に属するかを判定する計算セル生成プログラムである。ここでは、Fig.3 に対する計算セル生成の１例を示す。
まず，セルの輝度値 (yR , yG , yB ) を次のように定める：

f (x, y, t) =∂(∇ · v)
1+ ∇ · (v · ∇)v −Δ(∇ · v)
∂tReである。このとき，次のように計算する。

[数値計算アルゴリズム](4-1) 全領域 Ω において，次の形の方程式−∇ · (a(x, y, t)∇p) = f˜(x, y) in (0, T ) × Ω
を外周の境界条件 (7)，(8) のもとで数値解 p を求める。

(4-2) この数値解 p を Navier-Stokes 方程式:∂v
1+ (v · ∇)v = −∇p +Δv∂t
Rein (0, T ) × Ω

セルに含まれるピクセルの赤 (R) の輝度値の合計
yR =,Nc × Nc
セルに含まれるピクセルの緑 (G) の輝度値の合計,
yG =Nc × Nc
セルに含まれるピクセルの青 (B) の輝度値の合計.
yB =Nc × Ncに代入して，外周の境界条件 (5)，(6) および付加条件 (9)，(10)：

このセルの輝度値を用い，計算セル Ω̃1 , Ω̃2 , Ω̃3 を次のように定める：のもとで、速度ベクトル v を計算する。この計算では，陰解法と三次精度上流差分を用いている。

Ω̃2= { セルの集合；(yR , yG , yB ) = (0.0, 0.0, 0.0)},Ω̃3= { セルの集合；(yR , yG , yB ) = (0.0, 0.0, 255.0)},Ω̃1= { セルの集合；Ω\(Ω̃2 ∪ Ω̃3 )}.

v = (0.0, 0.0)

in (0, T ) × Ω̃2v = (0.2, 0.0)

in (0, T ) × Ω̃3(4-3) 以上の計算を各時間ステップ毎に繰り返す。
5.

数値計算結果

数値計算アルゴリズム (4-1),(4-2),(4-3) を適用する。
（１）3 つの円柱周りの流れの数値シミュレーション：
画像のピクセル数：4.8 × 105
セル=8 × 8 ピクセル
レイノルズ Re = 2, 000
時間の刻み幅 t = 0.01
時間ステップ数 s = 5, 000のもとで計算し，流れの方向を矢印で，流れの速さをカラーで表すと下図になる：

これより作られた計算セルは次の通りである。

Fig. 6: Regular cell for the computation
4.

計算セルと仮想領域法による数値計算アルゴリズム

２節と３節で述べた方法によって、3 つの円柱周りやリスボン市街地の風の流れの数値計算を行う。ここでは，リスボン市街地の場合を例に計算セルと仮想領域法を用いて Navier-Stokes 方程式に対する数値計算アルゴリズムを示す。
2 次元非圧縮性 Navier-Stokes 方程式の数値計算において，圧力ポアソン方程式を用いる。このとき，γ 上に3Fig. 7: Numerical result for the folw around 3 circularcylinder.

c 2011 by JSFM
Copyright 第 25 回数値流体力学シンポジウム
C03-1（２）リスボン市街地における風の流れの数値シミュレー
ション：
画像のピクセル数：3 × 105
セル=5 × 5 ピクセル
レイノルズ Re = 2, 000
時間ステップ数 s = 1, 000, 000
時間の刻み幅 t = 0.01のもとで計算し，流れの方向を矢印で，流れの速さをカラーで表すと下図になる：

いた。
参考文献

(1) H.Fujita, H.Kawahara and H.Kawarada, Distribution theoretic approach to ﬁctitious domain method
for Neumann problems, East-West J.Number
Math.3(2)(1995),pp.111-126.
(2) 河原田秀夫, 自由境界問題 (理論と数値解法), 東京大学出版 (1989).
(3) H.Koshigoe,T.Shiraishi and M.Ehara, Distributionalgorithm in ﬁnite diﬀerence method and its application to a 2D simulation fo temperature inversion, Journal of computational and appliedmathematics,vol.231(2009)pp.102-108.
(4) J.L.Lions,Perturbations Singulieres dansles Problems aux Limites et en Controle Optimal, Lecture
Notes in Mathematics 323,1973.

Fig. 8: Wind ﬂow in Lisbon city(5) 佐伯公美，小島智子，腰越秀之，地理情報を利用した差分法の数値計算-仮想領域と OpenGL を用いた可視化について-，日本応用数理学会 2007 年度年会
講演予稿集，pp.256-257(2007).

この数値シミュレーションからテージョ川から流入した風が大きな通りだけでなく小さな道にも入り，市街地の広場 B，C の方向に流れている様子がわかる。

(6) 吉井孝子，増田潤一，腰越秀之，デジタルカラーを用いた産業廃棄物中の流れの数値シミュレーション，第２２回数値流体力学シンポジウム，Proceedings CD，
K6-5(2008).

まとめ
仮想領域とデジタルカラー画像を用いて，２次元Navier-Stokes 方程式に対する数値計算法を示した。この計算では，次のプログラムを作成し，数値実験を行った：（C1）24bitBMP ファイルの画像から計算セルを生成するプログラム，（C1）の計算セ（C2）Navier-Stokes 方程式の計算では，
ルを用いて，LU 分解，３次精度上流差分，陰解法を基にした有限差分法のプログラム，（C3）流れの可視化では，OpenGL を用い (5) ，流れの向きはベクトルで、その大きさはカラーで表示するプログラム。
このように有限差分法はデジタル画像を組み込むことによって，計算セル上の数値計算を可能にしている。この方法は，現在，瀬戸内海・豊島の産業廃棄物・地層中の地下水流れの数値シミュレーション (6),(7) に応用されている。今後，下図のような有限な高さの円柱(7) T.Yoshii and H.Koshigoe,A computational method
for groundwater ﬂow through industrial waste byuse of digital color image,Lecture Notes in Computer Science 6329,Springer(2010),pp.288-296.

Fig. 9: Fig.5 Three circular cylinders
周りの流れやリスボン市街地の建物の高さを考慮した３次元の数値シミュレーションを仮想領域とデジタルカラー画像を用いて行う。
謝辞
本研究は科研費 (22540113) の助成を受けているものである。また，増田潤一氏（荏原製作所）からご援助を頂

4c 2011 by JSFM
Copyright