車内が臭いときにどの窓を開けたらいいのか

　車のどの窓を開けたら良いのか気になって調べてみたのですがいまいち良い答えにたどり着けませんでした。そこでオープンソースソフトウェア（OSS）を用いて車内の換気効率について計算し解析してみました。

最初に結論

　60km/hで走るヴィッツの窓は「前2つは全開、後ろは片方が全開、もう片方は半分くらい」で開けると最適に換気が行われる。

環境

• OS : Ubuntu 16.04 LTS
• モデリング : SALOME-MECA-8.3.0
• 流体計算 : OpenFOAM-v1712
• ポストプロセス : ParaView-5.4.0 (OpenFOAMv1712同梱)
• データ処理 : Jupyter lab, pandas

前提条件

　60km/hで走るTOYOTAのヴィッツを想定します。3Dモデルはarchibaseplanet.comからお借りしました。

　このモデルを使用して窓を開けたヴィッツの外と中の空気の流れを計算してみます。
空気の流れの計算にはOpenFOAMというOSSを使用しますが、そのためにまずSALOME-MECAというOSSでメッシュを作成します。（左がお借りしたモデル。右はそれを元に作成したメッシュ）

openfoam.org

www.code-aster.org

　今回の計算では窓の開け方のパターンをどう表現し整理するかという部分が重要なポイントになってきます。
そこで窓を以下のように4分割し、各部位で[空気が通過する or 壁になっている]と切り替えることで色々なパターンを作成しました。 具体的には4つの窓をそれぞれ[0%, 25%, 50%, 75%, 100%]の5段階ということで5×5×5×5-1=624通りの計算を行うこととなります。（全窓全閉は計算しない）

　なお条件のネーミングとして1111〜5554のように名付けています。
これは[右前, 右後, 左前, 左後]の開閉を表しており、[1:100% 2:75% 3:50% 4:25% 5:0%]のようになっています。
(例: 1225 →　[右前: 100%, 右後: 75%, 左前: 75%, 左後: 0%])

　OpenFOAM側の条件としては乱流モデルにRASのk-ωモデルを使用し、まず速度場の定常計算を行ったあとに、その速度場を用いて空気の新鮮さを非定常計算しました。

結果

　では624個の結果のうち一例を見てみましょう。

[case1111] 全窓を全開にした場合

　まずは全ての窓を全開にしたときの結果を見てみます。

　色々な情報を得ることができましたが、これだけでは優劣はわかりません。
我々の経験上、換気が早く、風の強いパターンであることは予測できます。
別の結果と比較をしながら考察を進めていきましょう。

図の説明

　この図のような結果が各開閉パターンで得られます。
今回OpenFOAMにより結果として得られたのはU(速度)、p(圧力)、ｋ,ω(乱流パラメータ)に加えてオリジナル物性値の新鮮さです。 左の図はまず車のCGに重ね合わせて空気の流れを線として描画しています。そして新鮮さで色付けをしています。青に近いほど新鮮な空気です。
　その上には全体平均・ドライバーの顔の横・助手席の人の顔の横の速度・新鮮さ・圧力を示しています。また右はヒストグラム（体積分布）です。平均値からはわからない値の偏りなどを示しています。

各物性値について

・U(速度)
　風速です。m/sで表しています。
　wikipediaによると1.5m/sまでは無風、3.5m/sまでは風を感じる程度、6m/sまでは木の葉を揺らす程度だそうです。

・p(圧力)
　今回の様な"速い流れと接した空間"はベンチュリー効果によって圧力が下がります。 圧力は低すぎると耳が痛くなることもありますのでこれは避けたいです。ちなみに飛行機の中は-20000Paだそうです。 ベンチュリ - Wikipedia

・ｋ,ω(乱流パラメータ)
　kとωは乱流に関するパラメータですが今回は注目しません。

・新鮮さ(空気が車内に入ってきてからの時間)
　空気が車内に入ってから何秒経っているかを新鮮さとして計算させました。 例えばこれが20秒であれば、そこにある空気は車内に入ってから20秒経過していることを示します。 この数字が大きいほど古い空気であるとしています。
※"空気齢"という言葉があるそうですが、私は専門外であり今回のとは違う概念である可能性がありますので、今回は"新鮮さ"という表現で記載したいと思います。

[case5515] 助手席のドアのみが全開の場合

　次に助手席側の窓のみが空いているときの結果を見てみましょう。

　平均の新鮮さが30秒も長くなっています。つまり空気が全窓全開よりも30秒古いということがわかります。 その割に助手席に座る人の顔に当たる風の強さは先程とあまり変わりませんのであまりいい窓の開け方とは言えません。

[case5551] 左後ろの窓のみが全開の場合

　次は後窓の片方を全開にした場合です。 　こちらは更に空気が古く、およそ5分前の空気が流れていることがわかります。

　このようにして全開閉パターンについて計算しましたが、この比較方法では「どれが最も良いのか」を判断することができません。 そこでpandasなどを用いてデータを整理して比較していきます。

• 全結果の比較

　以下に各パターンでの平均値の分布を示します。
縦軸は全て新鮮さ[s]であり、横軸には左に速度U、右に圧力pです。

　まず左の速度の結果を見ると、新鮮さと平均速度には相関があることがわかります。つまり、「換気はされるけどあまり風が吹かない」といった都合のいい条件がないことが確認できました。

　次に右の圧力の結果を見てみます。
最も低いパターンでも-75Pa程度でありこれは耳などに影響するほどのものではありません。 ということで今回は圧力の結果は無視してしまって良いと判断しました。

• 最も換気が早い条件

　まずは新鮮さでソートしてみます。

case	速度(平均)	速度(顔付近)	速度(最大)	速度(最小)	新鮮さ
1113	1.93056	0.931776	13.3781	0.030810	19.9936
2213	1.84915	0.929595	14.1396	0.034294	20.5013
2112	1.91422	0.618329	13.4901	0.018785	20.5772
2111	1.95503	0.632469	13.8089	0.030394	20.6549
1111	2.00505	1.269572	13.6723	0.024419	20.7344
1212	1.97671	0.755584	13.7637	0.041585	20.7677

　平均値で一番空気が新鮮なのは"ケース1113(1311)"であることがわかりました。
つまり前2つは全開、後ろは片方が全開、もう片方は半分くらい。という条件です。

　結論は出てしまいましたがもう少し掘り下げてみます。

　速度の結果について拡大してみましょう。

　これらは新鮮さを25秒以下に絞ったものです。ここに顔の横を通る空気の風速を追加て表示してみます。

　「平均の速度はそこそこでも顔に当たる風が強すぎる」といったケースもあるかと思い描画してみました。 中には顔付近の風だけ速くなってしまう結果も見られますが、先程の最も換気される条件である"1113"は顔の周りに強風が吹いているというわけではなさそうです。

　最後にこれらのそれぞれのヒストグラム（体積分布）を重ね合わせてみます。

　もうここまで近いデータだと差がわかりません。
というわけで先程の"1113”を今回の”最適な窓の開け方”としてしまいましょう。

再び結論

[case1113] 前2つは全開、後ろは片方が全開、もう片方は半分くらい

のパターンで最適に換気が行われる。

最後に

　使用したツールやテクニックなどの概要を少し掘り下げて紹介します。
更に詳しい方法はいずれQiitaに載せる予定です。

モデリング

　外気の流れについては上の窓のメッシュが同一になるように以下のようなメッシュを追加しております。
この3つのメッシュについてそれぞれOpenFOAM形式に変換したあと、OpenFOAMのmergeMeshにより一つのメッシュとして計算を行います。
そして窓のメッシュのうち気体の移動が行われる部分（開いている窓）についてはstitchMeshを用いて縫い合わせることで、各窓パターンを表現します。

速度場の計算(OpenFOAM::simpleFoam)

　以上のメッシュを使用してまずは速度場について計算を実行します。 今回使用するソルバーはsimpleFoamと呼ばれる"定常状態"の速度-圧力計算を行うためのソルバーです。 乱流モデルはRASのk-ωモデルを使用しました。 初期条件として時速60km/h≒17m/sの空気が前方より吹き込まれます。

空気の新鮮さの計算(OpenFOAM::ソルバーを作る)

　新鮮さを示す指標として、今回は空気が何秒間滞在しているかというオリジナルな物性値を求めてみました。
このような自由度の高い改造もOSSの魅力ですね。
入り口の空気を0秒として、simpleFoamで計算した速度場に乗って空気が流れていきます。
⊿t秒経過する毎にその新鮮さも⊿t秒ずつ加算されていき最終的に落ち着いた状態が空気の新鮮さの分布となる、といったコンセプトです。

OpenFOAMで記述すると以下のようになりました。

while (pimple.loop()) {
    fvScalarMatrix TEqn (
        fvm::ddt(T)
      + fvm::div(phi, T)
      - fvm::laplacian(DT, T)
    );

    TEqn.solve();
}

T += runTime.deltaT();

計算の実行

• 元のメッシュをコピー
• stitchMeshで窓を作成
• simpleFoamで速度場を計算
• 手作りソルバーで空気の新鮮さを計算

という動作を624パターンについて行うことになるため、シェルスクリプトを書いて一気に計算させました。

結果の後処理(paraview)

　各計算結果を表示するためにはParaViewを使います。
今回ParaViewをJupyterLab内で呼び出す方法に挑戦しました。
ParaViewはコンパイルする際にPythonのAPIをラッピングすることができるため、パスを通すことでimport paraviewをすることができます。 今回のように多くの結果を比較したい場合などは、ケースの名前を変更するだけで図を出力するようにしておくと便利です。

　また、

import paraview.simple as pv

foam = pv.OpenFOAMReader('path/to/case') # ケースを読み込む

probe = pv.ProbeLocation(Input=foam,ProbeType='Fixed Radius Point Source') # ある点の情報を読むこむ
probe.ProbeType.Center = [0.0, 0.0, 0.0]

p = probe.PointData['p'].GetRange(0)[0] # その点の値を読み込む。

# 注意：GetRangeはそのデータの[最小値, 最大値]を返します。
# 今回は一点のみの抽出ですので最大値最小値どちらも同じ値が返ってきます。

などのようにすればある点での値を読み込むことができます。

まだあまり使いこなしていないため以下のコードを見ながら勉強中です。 github.com ある程度理解してきたらどこかにまとめたいと思います。

結果の後処理（ユーティリティ）

　すでに計算したケースの値の平均値を計算させる方法を最初はParaViewとPyFoamで探していました。
ただ今回は見つからなかったためC++で簡単なコードを書いて計算させました。
pythonでcommands.getoutputで実行して返ってきた標準出力を使う、といった方式です。
時間がかかるため、この辺は今後やり方を考えなければいけないと思っています。

以下のようなものもありましたので今後試してみたいと思います。
github.com

結果の整理 (pandas)

　paraviewなどにより抽出した値はPythonのフレームワークであるpandasでデータベース化して利用しました。
今回のように大量のケースを比較する際にはとても有用でした。

最後に

　初ブログで拙い文章であるにも関わらずここまで読んでくれた人、本当にありがとうございます。
　今後は2ヶ月に一回くらい何か計算をして、それに関する記事をQiitaかここに上げていけたらと考えています。