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

VCA 法による室内汚染物質の放出源推定
Estimation of indoor contamination source with VCA method○ 松尾 智仁, 阪大院工, 〒565-0871 大阪府吹田市山田丘 2-1, matsuo@ea.see.eng.osaka-u.ac.jp
近藤 明 , 阪大院工, 〒565-0871 大阪府吹田市山田丘 2-1, kondo@see.eng.osaka-u.ac.jp
久野 貴大, 阪大院工, 〒565-0871 大阪府吹田市山田丘 2-1, kyuno@ea.see.eng.osaka-u.ac.jp
井上 義雄, 阪大院工, 〒565-0871 大阪府吹田市山田丘 2-1, inoue@see.eng.osaka-u.ac.jp
Tomohio Matsuo, Graduate School of Engineering, Osaka University, Yamadaoka 2-1 Suita, Osaka,
565-0871
Akira Kondo, Graduate School of Engineering, Osaka University, Yamadaoka 2-1 Suita, Osaka,
565-0871
Takahiro Kyuno, Graduate School of Engineering, Osaka University, Yamadaoka 2-1 Suita, Osaka,
565-0871
Yoshio Inoue, Graduate School of Engineering, Osaka University, Yamadaoka 2-1 Suita, Osaka,
565-0871
For assessing and managing the indoor thermal and hygienic environment, it is essential to understand flow, thermaland scalar fields. And for understanding these fields, measurements and computational fluid dynamics (CFD) areusually used. However, measurements may have some missing data and CFD may have errors caused by discretizationof governing equations or incorrect initial conditions. Thus, to obtain accurate data of entire field, Variational
Continuous Assimilation method (VCA method) was applied to indoor environment. The validity of the VCA method
for estimation of location of indoor unknown contamination was verified by simulations with quasi-observed data
created by CFD.

１．はじめに
室内の温熱環境や衛生環境を把握・改善するためには、室内の流れ場や、温度や物質濃度などのスカラー場を把握することが重要となる。従来、これらの場の把握手法としては、実測と数値解析（CFD）が主に用いられている。
しかし、これらの２つの方法にはそれぞれ長所と短所がある。
実測は観測誤差を含むために得られたデータは必ずしも支配方程式を満たさず、また対象領域全体の分布を計測することは困難であるために観測値の存在しない欠測点を生じる。
他方 CFD は数値モデル化に際して近似による誤差や離散化に伴う誤差を含む。そのため、これら実測と CFD を補うあうため、以下に示すように様々なデータ同化法が開発されてきた。
もっともシンプルなデータ同化手法として、観測値に整合するように計算値にナッジング項を加えることで場を修正するナッジング法がある。ナッジング法は適用が容易である一方、同化が観測値の存在する点のみで行われるため、観測値を重視すると観測点近傍で支配方程式を満足しなくなり、計算値を重視すると修正が十分に行われないという欠点がある(1)。そのため、同化のバランスを決定するナッジング係数の検討が重要である。観測値との差分に加え、ナビエ・ストークス方程式や連続の式といった支配方程式の残差を修正に用いる手法として費用関数法がある(2),(3)。
費用関数法は対象領域全体における計算―観測誤差と支配方程式残差を合計する費用関数を導出し、費用関数を最小化するように場を修正する。
場を修正する手法の他に、初期条件、境界条件を修正する手法として、adjoint 法がある。これは誤差二乗などの評価関数を用いて観測値と計算値のズレを評価し、評価関数を最小化するような初期条件・境界条件を推定するものである。adjoint 法は誤差の原因を初期・境界条件に求めるためモデル誤差や計算誤差を扱うことができない(4)。しかしモデルに誤差が含まれる場合、正しい初期条件を用いても計算結果には誤差が含まれてしまう(5)。

放出源推定に特化した手法として、以下に示すように物質の移流拡散方程式を逆向きに解く手法がある。quasi-reversibility 法は支配方程式を離散化した後、時間変化率を負にしてナビエ・ストークス方程式を解くことで、観測された濃度から放出源を逆推定する。一般に、逆問題においては拡散項を適切に扱うことが難しいため、quasi-reversibility 法は拡散項を４次の安定化項で置き換えるなどの工夫が施される(6),(7)。また、pseudo-reversibility 法は時間変化率ではなく流れを逆転させてナビエ・ストークス方程式を解くことで放出源位置を推定している(8)。特に 10μm 以下の粒子径の浮遊粒子における放出源推定において、pseudo-reversibility 法はquasi-reversibility 法より優れた推定を与える(9)。いずれの手法も汚染が放出された正確な時刻を必要とするため、適用範囲は限られる(9)。
本研究が採用する Variational Continuous Assimilation（VCA）法
は評価関数を用いて場を修正する点で費用関数法と類似している。
しかし、以下の章で説明するように、VCA 法はある時空間における観測値と計算値の間の差分を、その時空間だけでなくその時間的、空間的な風上側の修正にも用いる。そのため、観測点近傍だけでなく、その風上側の場も修正されるという特徴を持つ。
また、修正は場に対して適用されるため、モデル誤差に対する修正に適用することも可能であると考えられる。
２.VCA 法の概要
本研究が用いる VCA 法は、LeDimet と Talagrand が導入し(10)、Derber が修正したものである(4)。

12.1 VCA 法の基礎方程式
本来、VCA 法の操作変数はモデル内の時間変動するすべての変数である。VCA 法は他の変分問題と同様、変分問題の目的関数を定義し、その目的関数を最小化することにより解を得る。目的関Copyright © 2012 by JSFM

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

数の定め方は様々ありうるが、本研究では計算値と観測値の残差の二乗和によって定義する。

ここで、 B p は A p の線形化を意味する。

同化区間における観測値の存在するタイムステップが P 個あ式(6)を式(5)に代入することで式(7)を得る。

り、そのタイムステップがt 1  t 2    t P であるとすると、モとし、計算値の存在する最後のタイムステップには観測値が存在 p( p1 ) Bp
 p

する( t N  t P )必要がある。このとき、目的関数 I は式(1)で与え式(7)は漸化式である。また、初期時刻における計算値ベクトル

デルは N （  P ）個の計算値が存在するタイムステップを必要られる。

I

1 P~~
(  p   p )T (  p   p )
2 p1

は に依存しないため、

(1)

 p
  p  B p ( p1

値ベクトル、観測値ベクトルであり、 ( )T はベクトルまたは行列の転置を意味する。

 B p1 ( p2   B 4 (3  B 22 ) ))

数値計算では、前の計算ステップ t n1 における計算値ベクトルに何らかの変換を行うことで次の計算ステップ t n の計算値ベク  A (n 1)

pp 1 p  p   (  B m )ll 2 ml 1(2)

ここで A n は  n1 から  n への変換を意味する行列式である。

pp 1
PI~
  ( p   (  B m )l )T ( p   p ) p1l 2 ml 1

式(2)に補正項を付け加えたものが式(3)である。

 A   n 1n 1n(3)によって決定する補正ベクトルである。  は  と同じ次元を持

pp 1
P
I~
  ( p   l (  ( B m )T ))( p   p ) (11) p1l 2ml 1
ここで、adjoint 変数 を導入し、式(12)のように定義する。Nつ。
VCA 法は、式(1)の I を最小化するような補正ベクトル を求めることでモデルの誤差を修正する。

は計算ステップ回数で、計算最終ステップから時間的に遡った漸化式となっている。

2.2 評価関数 I の最小化

~nnn1 T n1n  ( B )    n (   ) N~NN    

評価関数 I を最小化するような を導出するためには、少なくI
 0 が成立する必要がある。これは局所的最

1 if tn  t p
 0 otherwise

n  I
最小点を求めることは困難であるため、
 0 が成立する点を
NI  n  n
 n 2式(1)を補正ベクトル で微分することで式(4)を得る。

A (p 1) p(14)以上により、評価関数 I の補正ベクトル に対する勾配を導出できる。VCA 方の適用においては、式(12)、式(13)によって adjoint(4)変数  を求め、式(14)により評価関数 I の勾配を求める。

また、式(3)を  で微分することで式(5)を得る。
p(13)

式(12)を式(11)に代入すると、式(14)を得る。

最小点とみなし、そのような を探索する。

PI
 p T p ~ p () (   ) p1 (12)ただし、 n は式(13)のように定義する。

小点、正確には停留点を求めているに過ぎない。しかし、大域的p(10)すなわち、ここで、  は時刻 t n によって決定される正方行列、 は空間nとも関数の勾配(9)式(9)を式(4)に代入すると、VCA 法では、各計算ステップでモデルそのものに補正を行う。
n 1(8)すなわち、トルが作られる。これを式(2)で表現する。
n ( 1 ) 0 である。


したがって漸化式を繰り返し適用することで式(7)は式(8)のように変形できる。

~

ここで、  p 、  p はそれぞれタイムステップ t p における計算n(7)評価関数 I の勾配が定まれば、最急降下法を用いた繰り返し計算により、評価関数 I を最小化する補正ベクトル を求めることができる。

(5)

ここで、式(5)の右辺第１項を式(6)のようにおく。
p 1
A p ( p 1 ))p  (B

2

Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-1
VCA 法は、３次元の非定常な場の推定に用いることができる。

ないという利点がある。
VCA 計算における濃度の初期値は対象領域全体で 0 であるか

しかし、本研究では簡単のために、２次元の定常な場を用いてら、VCA 計算結果の濃度分布は補正ベクトル  によって決定さ３． VCA 法の有効性の検証

VCA 法による放出源推定の検証を行った。

れる。前述の式(3)にあるように補正ベクトルは計算値ベクトルに対する生成消滅項のように振る舞うため、補正ベクトルの分布はソース項（放出源）とみなすことができるとしている。すなわち、また、VCA 法による修正の対象は濃度場のみでなく、気流場や温度場を修正することも可能である。しかし、本研究では正しい流れ場が与えられた場合の濃度の放出源推定を目的としているた補正ベクトル  の値が高い領域の近傍に放出源が存在すると推め、移流場や温度場については既知であるとして修正を行わない。

定することができる。

VCA 法の検証は CFD シミュレーションにより行った。VCA 法
は数値計算の修正のために観測データを必要とするが、本研究でシミュレーションでは、はじめ対象領域全体におおむね等間隔に観測点を配置して VCA 法を適用する。その後、得られた補正は正しい境界条件により計算された濃度場から観測点における濃度データを抽出することで擬似観測データを作成し用いた。

ベクトル  の分布から放出源位置を絞り込み、その近傍に観測点3.1 CFD 計算の概要を追加した上で VCA 法を再び適用する。以後この繰り返しにより、必要な精度まで放出源位置を絞り込んでゆく。
Fig.4 に、今回用いた観測点パターンを示した。Fig.4 (a) がはじめの観測点パターンであり、Fig.4 (b), (c) は VCA 法の適用結果から放出源位置を絞り込んで定めた観測点パターンである。

CFD 計算には SIMPLE 法による有限差分法を用い、移流項にべき乗則、乱流モデルとしては標準 k   モデルを採用した。
3.2 対象領域
対象領域の概略を Fig.1 に示す。対象領域は吹き出し口と吸い込み口、机や仕切りといった障害物、汚染物質の放出源を持つ室内を想定した。吹き出し口からは毎秒 0.16m の風が吹き出し、吸込口は自由開放となっている。障害物として、部屋の中央に机が置かれ、その上には薄い台が置かれている。また、部屋の天井からは障害物として 2 枚の仕切りが設置されている。汚染物質はFig.1 に示した空間から放出される。汚染の放出量は時間的に一定とした。
正しい境界条件を与えた際の CFD 計算結果における流れ場、濃度場を Fig.2、Fig.3 に示す。検証では、これらの結果を真値（正解値）であるとみなしている。
Fig.2 の正解値流れ場では、机やパネルといった障害物の影響でいくつかの渦が形成されていることが見て取れる。また、Fig.3 の正解値濃度場では、Fig.2 にあらわれている渦の影響により、汚染物質が吸い込み口に向かって流れるだけでなく、一部が図の左上にむかって巻き上げられていることがわかる。

Fig.1 Objective room

3.3 VCA 計算の概要
VCA 法は、室内の濃度場と汚染の放出源位置、放出強度が未知であるとし、流れ場が既知であるとして適用した。
観測値データは、Fig.3 の濃度から観測点における濃度を抜き出すことで作成した。非定常計算においては観測データは複数の時刻において用いることができるが、本研究では定常場を仮定しているため、観測データは CFD 計算の最終ステップにおけるものだけでよい。

Fig.2 Flow field (CFD result as correct value)

3.4 観測点の配置
一般に放出源位置は未知であるから、観測点は放出源がどの位置にあってもその位置が推定できるように配置されなければならない。そのため、観測点は局所的に配置するよりも、領域全体に等間隔に配置したほうが良い。
しかし、領域全体に観測点を配置する場合、観測点の数的密度を高めると観測点数が急速に増加し、観測コストが増大する。
そのため、本研究では複数回のデータ同化により放出源位置を絞り込みつつ、観測点を配置し直す手法を採用した。この手法では、少数の観測点から VCA 法の適用を開始し、放出源の近傍に徐々に観測点を追加していくため、観測点数があまり大きくならFig.3 Scalar field (CFD result as correct value)3Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-1
４．結果と考察
VCA 法の適用結果とその考察を以下に記す。

4.1 放出源の位置の推定
3.4 節で説明したように、補正ベクトルの分布を用いて放出源位置を推定することが出来る。補正ベクトルの値は移流拡散方程式における生成項のように振る舞うため、放出源は補正ベクトルの値が大きい領域内にあると考えられる。
観測点パターン 3 における、補正ベクトルの等高線と汚染の放出源、観測点の位置を Fig.6 に示す。図中の四角の枠で囲まれた領域が汚染の放出源位置であり、曲線が補正ベクトルの等高線である。
Fig.6 より、VCA 法の繰り返し適用により、観測点の密度程度の精度まで汚染放出源の位置を推定できていることが確認できる。

4.1 補正ベクトル の分布
VCA 法における補正ベクトルの分布と、観測点の位置を重ねた図を Fig.5 に示す。
Fig.5(a)ははじめの VCA 法の適用結果である。観測点 A の風上側から、観測点 B までの領域に高い値が分布していることがわかる。そのため、その領域に観測点を追加し( Fig.5(b) )、再び VCA法を適用する( Fig.6 (b) )。
Fig.6(b)から、観測点 C の風上側から観測点 D までの領域に高い値が分布している事がわかる。そのため、その領域にさらに観測点を追加し( Fig.5(c )、再び VCA 法を適用する( Fig.6 (c) )。
計３回の VCA 法の繰り返し適用により、補正ベクトルに大きな値の入っている領域が Fig.5(c) 程度まで絞りこまれている。

BA(a) Observation pattern 1(a) Observation pattern 1CD(b) Observation pattern 2(b) Observation pattern 2(c) Observation pattern 3
Fig.5 Correction vector field (VCA result) with the observation points(c) Observation pattern 3
Fig.4 Observation points patterns

4Copyright © 2012 by JSFM

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

Fig.6 Contour lines of correction vector, contamination source, andobservation points(a) Observation pattern 1

4.2 濃度分布
VCA 法を適用し得られた濃度場を Fig.7 に示す。
Fig.7 と Fig.3 の比較より、VCA 法の繰り返し適用により、濃度分布の推定精度が向上していることが見て取れる。Fig.7(c) とFig.3 を比較すると、Fig.7(c) では濃度の分布がやや風下側にずれているが、これは前節で説明したように、補正ベクトルの分布がずれているためである。そのため、VCA 法のさらなる繰り返し適用によりさらに精度を高めることができる。
VCA 法の適用による濃度分布の修正精度を評価するため、正解値と VCA 計算結果の RMSE 値を求め、Table.1 に示した。
Table.1 より、VCA 法の繰り返し適用により、VCA 計算結果が正解値に近づいていることがわかる。

(b) Observation pattern 2

Table.1 RMSE of scalar field (VCA results)
Initial value

Observation
pattern.1

Observation
pattern.2

Observation
pattern.3

8.51

3.21

2.20

1.74

RMSE

５．結論
本研究では、Variational Continuous Assimilation (VCA) 法を用いた室内の濃度場と汚染物質の放出源位置の推定を行った。

(c) Observation pattern 3
Fig.7 Scalar field (VCA result)

VCA 法の適用に用いる CFD 計算には SIMPLE 法による有限差分法を用い、移流項にべき乗則、乱流モデルとしては標準 k   モ

参考文献(1) 加賀昭和, 塩田卓, 近藤明, 井上義雄, 山口克人, "測定値とCFD の融合のための費用関数法の非定常流れ場への適用,"
日本機械学会論文集(B 編) 682 (2003), pp. 1307-1312.
(2) 塩田卓, 加賀昭和, 近藤明, 山口克人, 井上義雄, "観測値と流体基礎方程式を融合した室内空気環境の把握," 空気調和・衛生工学会論文集, 79 (2000), pp. 29-36.
(3) 山口克人, 加賀昭和, 近藤明, 井上義雄, 山口太郎, 塩田卓, "費用関数を用いた最適化手法による PIV と CFD の融合," 日本機械学会論文集(B 編) 642 (2000), pp. 339-345.
(4) John C. Derber, "A Variational Continuous Assimilation
Technique," Monthly Weather Review 117(1989), pp. 2437-2446.
(5) Zoltan Toth, Malaquias Pena, "Data assimilation and numerical
forecasting with imperfect models: The mapping paradigm,"
Physica D 230 (2007), pp. 146-158.
(6) Tengfei Zhang, Qingyan (Yan) Chen, "Identification ofContaminant Sources in Enclosed Spaces by a Single Sensor,"Indoor Air 17(2007), pp. 439-449.
(7) Tengfei Zhang, Qingyan (Yan) Chen, "Identification of contaminantsources in enclosed environments by inverse CFD modeling,"Indoor Air 17(2007), pp.167-177.

デルを採用した。また、VCA 法の精度の検証のため、正解値となる濃度場を数値計算により作り出した。VCA 法に用いる観測値データは、正解値濃度場から観測点における濃度を抽出することで作り出した。
放出源位置の推定は、VCA 法の適用によって得られる補正ベクトルの分布を放出源位置の分布であるとみなすことで行った。シミュレーションによる精度の検証の結果、VCA 法の繰り返し適用により、観測点密度程度まで精度よく放出源位置を推定することができた。

5Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 B11-1(8) Shungo Kato, Yoshizumi Kajii, Ryunosuke Itokazu, Jun Hirokawa,Seiichiro Koda, Yoshikatsu Kinjo, "Transport of atmospheric carbon
monoxide, ozone, andhyd rocarbons from Chinese coast to
Okinawa island in the Western Pacific during winter," Atmospheric
Environment 38 (2004), pp. 2975-2981.
(9) Tengfei (Tim) Zhang, Hongzhu Li, Shugang Wang, "Inversely
tracking indoor airborne particles to locate their release sources,"
Atmospheric Environment 55 (2012), pp. 328-338.
(10) FRANÇOIS-XAVIER LE DIMET, OLIVIER TALAGRAND,“Variational algorithms for analysis and assimilation of
meteorological observations: theoretical aspects,” Tellus 38A(1986),
pp.97-110.

6Copyright © 2012 by JSFM