(57)【要約】【課題】複数の建物及び複数の樹木が存在するシミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を高精度に算出できるシミュレーション装置を提供する。【解決手段】シミュレーション装置を、複数の建物及び複数の樹木が存在するシミュレーション対象空間内の各要素が吸収する放射熱の効果、及び、樹木−大気間の熱・水蒸気交換の効果を考慮に入れて当該シミュレーション対象空間における非定常な風速・気温・水蒸気の空間分布を算出する装置として構成する。【選択図】図２
2015094751.tif 000002
【特許請求の範囲】
【請求項１】
複数の建物及び複数の樹木が存在するシミュレーション対象空間を表すデータから、当該シミュレーション対象空間を、各２要素間で放射熱が授受される複数の要素として取り扱うためのシミュレーション用データを生成する生成手段と、
前記生成手段により生成されたシミュレーション用データに基づき、前記複数の要素の中の、順番を考慮した各２要素に関する形態係数を算出する形態係数算出手段と、
前記形態係数算出手段により算出された形態係数を用いて各要素が吸収する正味の放射熱を算出する第１算出手段と、
前記シミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を算出する第２算出手段と、
前記シミュレーション対象空間内の各樹木の葉の温度を、葉が吸収する正味の放射熱と葉が存在する空間の大気の温度及び当該大気中の水蒸気分圧と葉表面の飽和水蒸気分圧とを用いて算出する第３算出手段と、
前記第１〜第３算出手段を繰り返し機能させることにより、前記シミュレーション対象空間における風速、温度及び水蒸気量についての、各樹木の樹冠と大気との間の熱及び水蒸気の交換量を考慮した非定常な空間分布を算出する第４算出手段と、
を備えることを特徴とするシミュレーション装置。
【請求項２】
前記第３算出手段は、各樹木の葉の温度を陰解法により算出する
ことを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
前記生成手段は、前記複数の樹木の樹冠を体積要素として取り扱うためのシミュレーション用データを生成し、
前記形態係数算出手段は、前記シミュレーション用データに基づき、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する
ことを特徴とする請求項１又は２に記載のシミュレーション装置。
【請求項４】
コンピュータが、
複数の建物及び複数の樹木が存在するシミュレーション対象空間を表すデータから、当該シミュレーション対象空間を、各２要素間で放射熱が授受される複数の要素として取り扱うためのシミュレーション用データを生成する生成処理と、
前記生成処理により生成されたシミュレーション用データに基づき、前記複数の要素の中の、順番を考慮した各２要素に関する形態係数を算出する形態係数算出処理と、
前記形態係数算出処理により算出された形態係数を用いて各要素が吸収する正味の放射熱を算出する第１算出処理と、前記シミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を算出する第２算出処理と、前記シミュレーション対象空間内の各樹木の葉の温度を、葉が吸収する正味放射熱と葉が存在する空間の大気の温度及び当該大気中の水蒸気分圧と葉表面の飽和水蒸気分圧とを用いて算出する第３算出処理とを含む主処理であって、前記シミュレーション対象空間における風速、温度及び水蒸気量についての、各樹木の樹冠と大気との間の熱及び水蒸気の交換量を考慮した非定常な空間分布が算出されるように前記第１〜第３算出処理を繰り返す主処理と、
を行うことを特徴とするシミュレーション方法。
【請求項５】
コンピュータに、
複数の建物及び複数の樹木が存在するシミュレーション対象空間を表すデータから、当該シミュレーション対象空間を、各２要素間で放射熱が授受される複数の要素として取り扱うためのシミュレーション用データを生成する生成処理と、
前記生成処理により生成されたシミュレーション用データに基づき、前記複数の要素の
中の、順番を考慮した各２要素に関する形態係数を算出する形態係数算出処理と、
前記形態係数算出処理により算出された形態係数を用いて各要素が吸収する正味の放射熱を算出する第１算出処理と、前記シミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を算出する第２算出処理と、前記シミュレーション対象空間内の各樹木の葉の温度を、葉が吸収する正味放射熱と葉が存在する空間の大気の温度及び当該大気中の水蒸気分圧と葉表面の飽和水蒸気分圧とを用いて算出する第３算出処理とを含む主処理であって、前記シミュレーション対象空間における風速、温度及び水蒸気量についての、各樹木の樹冠と大気との間の熱及び水蒸気の交換量を考慮した非定常な空間分布が算出されるように前記第１〜第３算出処理を繰り返す主処理と、
を行わせることを特徴とするシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション装置とシミュレーション方法とシミュレーションプログラムとに関する。
【背景技術】
【０００２】
都市部の気温等をシミュレートするためのシミュレーション技術として、以下のようなもの（特許文献１、２及び非特許文献１〜１５参照）が知られている。
・建物や樹木の３次元的な分布を考慮しながら、放射熱フラックスと表面温度のみを計算する。
・建物や樹木の３次元的な分布を考慮しながら、非定常な風の流れのみを計算する。
・建物の３次元分布のみを考慮しながら、非定常な風の流れと気温・水蒸気の分布を計算する。
・建物や樹木と大気との間の熱輸送を考慮し、定常な風の流れと気温・水蒸気分布を計算する。
【０００３】
しかしながら、樹木が存在する都市部における風速、気温、水蒸気の３次元空間分布は、放射熱フラックス、樹木大気間の熱・水蒸気交換などの効果によって建物壁面温度や葉表面温度などの影響を受け、敏感に変化する。そのため、それらの効果の一部しか考慮していない既存のシミュレーション技術では、都市部における非定常な流速・気温・水蒸気の空間分布を、適切に、シミュレートすることができない。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００３−０９９６９７号公報
【特許文献２】特開２００７−０７２９９１号公報
【非特許文献】
【０００５】
【非特許文献１】吉田 伸治, 村上 周三, 持田 灯, 大岡 龍三, 富永 禎秀, 「環境緩和効果を総合的に組み込んだ新しい3次元樹木モデルの開発」, 生産研究, 51巻1号, 1999年1月
【非特許文献２】吉田 伸治, 大岡 龍三, 持田 灯, 富永 禎秀, 村上 周三, 「樹木モデルを組み込んだ対流・放射・湿気輸送連成解析による樹木の屋外温熱環境緩和効果の検討」, 日本建築学会計画系論文集, No.536, pp.87-94, 2000
【非特許文献３】阪田 升, 鈴木弘孝, 謝 俊民, 「緑化を考慮した都市ヒートアイランドのCFD」, 計算力学講演会講演論文集, No.17, pp.637-638, 2004
【非特許文献４】坂本 雄三, 小島 悦史, 足永 靖信, 今野 雅, 「CFDを利用した樹木のクールスポット効果の数値解析 その１ ―樹木における放射と蒸散に関する計算モデル」, 日本建築学会大会学術講演梗概集, D-1, pp.689-690, 2005
【非特許文献５】小島 悦史, 坂本 雄三,足永 靖信, 今野 雅, 「CFDを利用した樹木のクールスポット効果の数値解析 その２ ―クールスポット効果のケーススタディ」, 日本建築学会大会学術講演梗概集, D-1, pp.691-692, 2005
【非特許文献６】大黒 雅之, 森川 泰成, 「街区及び敷地レベルを対象としたヒートアイランド解析評価システムの開発」, 大成建設技術センター報, No.38, pp.14-1- 14-8, 2005
【非特許文献７】大黒 雅之, 森川 泰成, 本橋 比奈子, 「街区スケールを対象としたヒートアイランド解析評価プログラムの開発と適用」, 大成建設技術センター報, No.42, pp.49-1- 49-8, 2009
【非特許文献８】片岡 浩人, 大塚 清敏, 赤川 宏幸, 「屋外熱環境評価のための数値予測モデルの開発 −樹木による冷却効果のモデル化−」, 日本建築学会学術講演梗概集, D-1, pp.927-928, 2008
【非特許文献９】佐々木 澄, 「数値解析に基づく街路樹がストリートキャニオン内の熱空気環境に及ぼす影響の検討」, 清水建設研究報告, No.85, pp.41-50, 2007
【非特許文献１０】H.B.リジャル, 大岡 龍三 他, 「数値解析による大規模緑地のヒートアイランド緩和効果の検討」, 東京大学生産技術研究所, 生産研究, 62巻1号, 2010
【非特許文献１１】神田 学, 稲垣 聡, 日野 幹雄, “植生-大気境界面における大規模渦構造と運動量交換に関するLESモデルによる検討”, 土木学会論文集 No.461/II-22, pp.39-48, 1993.2
【非特許文献１２】R.H. Shaw and U. Schumann, “Large-Eddy Simulation of Turbulent Flow Above and Within a Forest”, Boundary-Layer Meteorology, Vol.61, pp.47-64, 1992
【非特許文献１３】G. Bohrer et al., “Exploring the Effects of Microscale Structural Heterogeneity of Forest Canopies Using Large-Eddy Simulations”, Boundary-Layer Meteorology, Vol.132, pp.351-382, 2009
【非特許文献１４】F. Schlegel et al., “Large-Eddy Simulation of Inhomogeneous Canopy Flows Using High Resolution Terrestrial Laser Scanning Data”, Boundary-Layer Meteorology, Vol.142, pp.223-243, 2012
【非特許文献１５】J. Dauzat et al., “Simulation of leaf transpiration and sap flow in virtual plants: model description and application to a coffee planteation in Costa Rica”, Agricultural and Forest Meteorology, Vol.109, pp.143-160, 2001
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明の課題は、複数の建物及び複数の樹木が存在するシミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を高精度に算出できる技術を提供することにある。
【課題を解決するための手段】
【０００７】
上記課題を解決するために、本発明のシミュレーション装置は、
複数の建物及び複数の樹木が存在するシミュレーション対象空間を表すデータから、当該シミュレーション対象空間を、各２要素間で放射熱が授受される複数の要素として取り扱うためのシミュレーション用データを生成する生成手段と、
前記生成手段により生成されたシミュレーション用データに基づき、前記複数の要素の中の、順番を考慮した各２要素に関する形態係数を算出する形態係数算出手段と、
前記形態係数算出手段により算出された形態係数を用いて各要素が吸収する正味の放射熱を算出する第１算出手段と、
前記シミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を算出する第２算出手段と、
前記シミュレーション対象空間内の各樹木の葉の温度を、葉が吸収する正味の放射熱と葉が存在する空間の大気の温度及び当該大気中の水蒸気分圧と葉表面の飽和水蒸気分圧とを用いて算出する第３算出手段と、
第１〜第３算出手段を繰り返し機能させることにより、前記シミュレーション対象空間における風速、温度及び水蒸気量についての、各樹木の樹冠と大気との間の熱及び水蒸気の交換量を考慮した非定常な空間分布を算出する第４算出手段と、
を備える。
【０００８】
すなわち、本発明のシミュレーション装置は、複数の建物及び複数の樹木が存在するシミュレーション対象空間内の各要素が吸収する放射熱の効果、及び、樹木−大気間の熱・水蒸気交換の効果を考慮に入れて当該シミュレーション対象空間における非定常な風速・気温・水蒸気の空間分布を算出する構成を有する。従って、本発明のシミュレーション装置によれば、複数の建物及び複数の樹木が存在するシミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を、高精度に（既存のシミュレーション技術よりも適切に）、算出できることになる。
【０００９】
本発明は、上記シミュレーション装置と同様の特徴を有するシミュレーション方法や、情報処理装置（コンピュータ）を上記シミュレーション装置として機能させるシミュレーションプログラムとしても実現することが出来る。尚、本発明を、上記シミュレーションプログラムを記録したコンピュータ可読媒体としても実現することも出来る。
【発明の効果】
【００１０】
本発明によれば、複数の建物及び複数の樹木が存在するシミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を高精度に算出できる。
【図面の簡単な説明】
【００１１】
【図１】本発明の一実施形態に係るシミュレーションシステムの構成図
【図２】実施形態に係るシミュレーション装置の機能の説明図
【図３】２面積要素に関する形態係数についての算出式中のパラメータの説明図
【図４】面積要素から体積要素をみた形態係数についての算出式中のパラメータの説明図
【図５】体積要素の熱収支の説明図
【発明を実施するための形態】
【００１２】
以下、図面を参照して、本発明の一実施形態について説明する。尚、以下の実施形態の説明は、例示であり、本発明は、実施形態の構成に限定されない。
【００１３】
図１に、実施形態に係るシミュレーションシステムの構成を示す。図示してあるように、本実施形態に係るシミュレーションシステムは、シミュレーション装置１０、入力装置２０及び出力装置３０を含む。また、シミュレーション装置１０は、演算部１２、記憶部１４、インタフェース回路（Ｉ／Ｆ）１６を備える。
【００１４】
シミュレーション装置１０のインタフェース回路（Ｉ／Ｆ）１６は、演算部１２が他装置と通信を行うために使用する回路である。記憶部１４は、シミュレーションプログラム１８を記憶した不揮発性の記憶装置である。この記憶部１４は、演算部１２が利用するデータや演算部１２による処理結果の記憶にも使用される。
【００１５】
演算部１２は、ＣＰＵ(Central Processing Unit)、ＲＡＭ(Random Access Memory)等
を組み合わせたユニットである。この演算部１２は、シミュレーションプログラム１８を記憶部１４から読み出して実行することにより、本発明に係る生成手段及び第１〜第４算出手段として機能して各種処理（詳細は後述）を行う。
【００１６】
入力装置２０は、シミュレーション装置１０に情報を入力するための装置である。入力装置２０は、キーボード、マウス等のポインティングデバイス等を含む。また、出力装置３０は、シミュレーション装置１０からの情報を出力するための、ＬＣＤ(Liquid Crystal Display)やＣＲＴ(Cathode-Ray Tube)等のディスプレイ、プリンタ等である。
【００１７】
尚、シミュレーション装置１０は、通常、行列演算が高速に行えるコンピュータ（ベク
トル型並列計算機等）にシミュレーションプログラム１８を実行させることによって実現される。従って、通常、シミュレーション装置１０（ベクトル型並列計算機等）に接続された入出力用のコンピュータの入力装置及び出力装置が、入力装置２０及び出力装置３０として機能する。
【００１８】
以下、シミュレーション装置１０の機能を説明する。
シミュレーション装置１０は、複数の樹木を含む都市空間（以下、シミュレーション対象空間と表記する）内の風速、温度（大気温度、葉表面温度、建物表面温度等）及び水蒸気量についての、各樹木の樹冠と大気との間の熱及び水蒸気の交換量を考慮した非定常な空間分布を算出するための装置である。
【００１９】
図２に模式的に示してあるように、シミュレーション装置１０の使用時には、シミュレーション装置１０に、構造指定データ、初期値データ及び計算条件データが入力される。尚、シミュレーション装置１０に入力された各データは、一旦、記憶部１４（図１）に記憶されてから演算部１２により読み出されて利用される。
【００２０】
シミュレーション装置１０に入力される構造指定データは、シミュレーション対象空間の構造（３次元都市形状及び３次元樹木分布）を指定するデータである。この構造指定データとしては、地面の高さの２次元平面データ（地面の高さと座標との対応関係を表すデータ）と、建物の高さの２次元平面データ（建物の高さと座標との対応関係を表すデータ）と、樹木インデックス（樹木の識別情報）の２次元平面データと、各樹木インデックスで識別される樹木の葉面積密度の鉛直分布を示すデータとを含むデータを使用することができる。また、構造指定データとして、各樹木の葉面積密度に関する情報を含む３次元ＣＡＤ（Computer Aided Design）データや、各樹木の葉面積密度に関する情報を含まない
３次元ＣＡＤデータと各樹木の葉面積密度に関するデータとを含むデータを使用することもできる。
【００２１】
初期値データは、シミュレーション対象空間の各部の温度や水蒸気量の初期値を指定するためにシミュレーション装置１０に入力されるデータである。この初期値データとしては、各部の温度や水蒸気量の初期値を構造指定データと同じ座標系で指定するデータや、計算対象となる計算格子／平面要素／体積要素（詳細は後述）単位で温度や水蒸気量の初期値を指定するデータを使用することができる。
【００２２】
シミュレーション装置１０に入力される計算条件データは、シミュレーション対象時間範囲内の太陽の高度及び方位を表す太陽方位ベクトル情報等の、シミュレーションに必要とされる各種情報を含むデータである。尚、計算条件データの要素としてシミュレーション装置１０に入力される太陽方位ベクトル情報は、通常、シミュレーション対象時間範囲内の各時刻の（時間Δｔ経過毎の）太陽方位ベクトルを含む情報とされる。ただし、構造指定データとして各部の緯度及び経度がわかるデータを入力する場合、そのような太陽方位ベクトル情報を計算条件データに含める代わりに、シミュレーション対象時間範囲の開始年月日時刻、緯度及び経度から太陽方位ベクトルを求めるサブルーチンをシミュレーションプログラム１８に含めておくこともできる。
【００２３】
以下、上記データに基づき、シミュレーション装置１０（演算部１２）が行う処理の内容を説明する。
【００２４】
図２に示してあるように、シミュレーション装置１０の演算部１２が行う処理は、前処理と主処理とに分類できる。
【００２５】
《前処理》
前処理は、シミュレーション用データ生成処理（ステップＳ１０１）とパラメータ算出処理（ステップＳ１０２）とがこの順に行われる処理である。
【００２６】
・シミュレーション用データ生成処理
ステップＳ１０１にて行われるシミュレーション用データ生成処理は、入力された構造指定データからシミュレーション用データを生成する処理である。
【００２７】
シミュレーション用データ生成処理により生成されるシミュレーション用データは、シミュレーション対象空間（入力された構造指定データが示している３次元都市空間）を主処理中の風速分布等算出処理に用いる複数の計算格子（３次元構造格子）の集合として取り扱えるようにするためのデータである。また、シミュレーション用データは、シミュレーション対象空間内の各建物の壁面や地面を複数の面積要素として取り扱えるように、且つ、シミュレーション対象空間内の各樹木の樹冠を１つ以上の透過性を有する体積要素として取り扱えるようにするためのデータともなっている。尚、本実施形態に係るシミュレーション用データ生成処理により生成されるシミュレーション用データは、各体積要素がいずれかの計算格子と一致し、各面積要素がいずれかの計算格子の一面と一致するデータである。
【００２８】
このシミュレーション用データの具体的な内容（データ構造）は、計算格子と面積／体積要素の対応関係等、後述する各種演算に必要な情報が分かるものであれば良い。従って、シミュレーション用データを、例えば、面積／体積要素毎に、『通し番号、対応するシミュレーション対象空間内の計算格子の座標番号、面積要素の向いている方向／体積要素であることを示すフラグ、樹冠であるか否かを示すフラグ等からなるデータ』を含むデータとしておくことができる。
【００２９】
・パラメータ算出処理
ステップＳ１０２にて行われるパラメータ算出処理は、シミュレーション用データ生成処理により生成されたシミュレーション用データに基づき、主処理で使用される各種パラメータを算出する処理である。
【００３０】
このパラメータ算出処理時に算出されるパラメータには、各２要素（面積／体積要素）に関する形態係数Ｆ、各体積要素ｋの有効表面積＜Aeff＞ｋ、各要素ｉの天空率ωｉ、各面積要素ｉの日陰率Ｄｉ、各体積要素ｋの有効日陰率Ｄeffkがある。
【００３１】
まず、パラメータ算出処理時に算出される各２要素に関する形態係数Ｆについて説明する。
【００３２】
パラメータ算出処理時には、２つの面積要素ｉ，ｊの組み合わせ毎に、『面積要素ｉから面積要素ｊを見た形態係数Ｆｉｊ』と『面積要素ｊから面積要素ｉを見た形態係数Ｆｊｉ』とが、算出される。また、面積要素ｉと体積要素ｋの組み合わせ毎に、『面積要素ｉから体積要素ｋを見た形態係数Ｆｉｋ』と『体積要素ｋから面積要素ｉを見た形態係数Ｆｋｉ』とが、算出される。さらに、２つの体積要素ｋ及びｌの組み合わせ毎に、『体積要素ｋから体積要素ｌを見た形態係数Ｆｋｌ』及び『体積要素ｌから体積要素ｋを見た形態係数Ｆｌｋ』が算出される。
【００３３】
パラメータ算出処理時に算出される『面積要素ｉから面積要素ｊを見た形態係数Ｆｉｊ』、『面積要素ｊから面積要素ｉを見た形態係数Ｆｊｉ』は、それぞれ、以下の（１）、（２）式で定義される値である。
【００３４】
【数１】
2015094751.tif 000003
【００３５】
これらの式において、Ａｉ、Ａｊは、それぞれ、面積要素ｉの面積、面積要素ｊの面積である。βｉ、βｊ は、図３に模式的に示してあるように、それぞれ、面積要素ｉ内の
微小面ｄＡｉと面積要素ｊ内の微小面ｄＡｊとを結ぶ直線と微小面ｄＡｉの法線のなす角度、当該直線と微小面ｄＡｊの法線のなす角度である。また、ｒは、微小面ｄＡｉ・微小面ｄＡｊ間の距離である。
【００３６】
Ｔｉｊは、面積要素ｉと面積要素ｊの間の透過率である。Ｔｉｊは、２微小面ｄＡｉ及びｄＡｊ 間の光学的厚さτｉｊを用いて次式により算出される。
【００３７】
【数２】
2015094751.tif 000004
【００３８】
また、面積要素ｉと面積要素ｊとの間（微小面ｄＡｉと微小面ｄＡｊとの間）に樹冠が分布している場合、光学的厚さτｉｊは、樹冠の消散係数ｋｅｘｔ及び葉面積密度ａを用いて次式により算出される。
【００３９】
【数３】
2015094751.tif 000005
【００４０】
尚、形態係数の具体的な算出手順については後述するが、上記した（１）、（２）式で定義される形態係数は、相反則を満たしている。すなわち、面積Ａｉと面積Ａｊと形態係数Ｆｉｊと形態係数Ｆｊｉとの間には、以下の関係がある。
【００４１】
【数４】
2015094751.tif 000006
【００４２】
従って、Ｆｊｉを、（１）式から算出したＦｉｊとＡｉとＡｊとから算出することも、Ｆｉｊを、（２）式から算出したＦｊｉとＡｉとＡｊとから算出することもできる。
【００４３】
パラメータ算出処理時に算出される『面積要素ｉから体積要素ｋを見た形態係数Ｆｉｋ』は、以下の（３）式により定義される値である。
【００４４】
【数５】
2015094751.tif 000007
【００４５】
この（３）式におけるβｉは、図４に模式的に示してあるように、面積要素ｉの微小面ｄＡｉと体積要素ｊ内の微小投影面ｄＡｋ←ｉとを結ぶ直線と微小面ｄＡｉの法線のなす角度である。また、ｒは、微小面ｄＡｉ・微小投影面ｄＡｋ←ｉ間の距離である。
【００４６】
Ａｅｆｆｋ←ｉは、体積要素ｋ自身の遮蔽率を考慮した，面積要素ｉの方向から見た体積要素ｋの有効面積である。このＡｅｆｆｋ←ｉは、次式により算出される。
【００４７】
【数６】
2015094751.tif 000008
【００４８】
ここで、Δτｋ←ｉは、微小投影面ｄＡｋ←ｉ(図４参照)に垂直な方向の体積要素ｋの光学的厚さである。体積要素ｋが樹木の場合、Δτｋ←ｉは、消散係数ｋｅｘｔ、葉面積密度ａおよびｄＡｋ←ｉに垂直な方向の体積要素ｋの幾何学的厚さΔｓｋ←ｉを用いて、次式により算出される。
【００４９】
【数７】
2015094751.tif 000009
【００５０】
要するに、上記した（３）式は、（４）式を用いれば、以下のように表記できる式となっている。パラメータ処理時には、以下の（５）式に従って、面積要素ｉから体積要素ｋを見た形態係数Ｆｉｋが算出される。
【００５１】
【数８】
2015094751.tif 000010
【００５２】
パラメータ算出処理時に算出される『体積要素ｋから面積要素ｉを見た形態係数Ｆｋｉ』は、次式により定義される値である。
【００５３】
【数９】
2015094751.tif 000011
【００５４】
換言すれば、体積要素ｋから面積要素ｉを見た形態係数Ｆｋｉとしては、以下の式で表される相反則を満たす値が算出される。
【００５５】
【数１０】
2015094751.tif 000012
【００５６】
パラメータ算出処理時に算出される『体積要素ｋから体積要素ｌを見た形態係数Ｆｋｌ』は、次式により定義される値である。
【００５７】
【数１１】
2015094751.tif 000013
【００５８】
すなわち、この形態係数Ｆｋｌは、体積要素ｋの遮蔽率（“1-exp(-Δτｋ←l)”）と
体積要素ｌの遮蔽率（“1-exp(-Δτl←ｋ)”）と体積要素ｋ、ｌ間の透過率Ｔｋｌとを
考慮に入れた形態係数となっている。
【００５９】
尚、（７）式が表している体積要素に関する形態係数は、面積要素の形態係数と同様に次式の相反則を満たす。
【００６０】
【数１２】
2015094751.tif 000014
【００６１】
従って、他の形態係数と同様に、２体積要素間の形態係数の双方を、（７）式に従って算出することも、一方の形態係数を（７）式に従って算出し、他方の形態係数を、当該一方の形態係数の算出結果から算出することも出来る。
【００６２】
パラメータ算出処理時には、上記した各形態係数が、モンテカルロ法により算出される。
【００６３】
すなわち、形態係数の算出時には、０から１までの範囲の一様乱数Ｒθ、Ｒφから、以下の式によりμ、φを算出し、算出結果に基づき、以下に示した内容の単位ベクトルｎを生成する処理が、必要とされる形態係数の精度に応じた回数、繰り返される。
【００６４】
【数１３】
2015094751.tif 000015
【００６５】
換言すれば、形態係数の算出時には、Lambertの余弦則にしたがって面積要素（又は体
積要素）から射出される多数の光子の進行方向単位ベクトルｎを生成する処理が行われる。
【００６６】
そして、生成した各単位ベクトルｎの方向に要素ｉから等しいエネルギーＷ０を持った
光子が射出された場合に要素ｊに入射されることになる各光子ｐのエネルギーＷｐを積算することにより、形態係数Ｆｉｊが求められる。
【００６７】
具体的には、面積要素ｉから面積要素ｊをみた形態係数Ｆｉｊは、面積要素ｉから射出させた光子の個数がＮ、面積要素ｊに入射された光子の個数がｎであった場合には、以下の式により算出される。
【００６８】
【数１４】
2015094751.tif 000016
【００６９】
なお、建物壁面などの完全遮蔽性の面積要素のみの場合には、面積要素ｊに入射する光子のエネルギーＷｐはＷ０と等しいため、形態係数は、Ｆｉｊ＝ｎ／Ｎにより算出できる。樹木のような放射透過性の要素が空間中に分布する場合には、面積要素ｊに入射する光子のエネルギーＷｐは、面積要素ｉから面積要素ｊに到達するまでの間に減衰する。そのため、その減衰の形態係数に対する影響を光子の持つエネルギーを減衰させることで考慮する。すなわち、次式により、Ｗｐが計算される。
【００７０】
【数１５】
2015094751.tif 000017
【００７１】
ここで、Ｔｉｊ，ｐは、光子ｐの経路に沿った透過率である。なお、光子が透過性の体積要素内で確率的に遮蔽されると考えて、面積要素ｊに入射する光子の数を減少させても、減衰の影響を考慮した形態係数を求めることが出来る。
【００７２】
また、面積要素ｉから体積要素ｋをみた形態係数Ｆｉｋは、（７）式により定義される値である。従って、面積要素ｉから体積要素ｋをみた形態係数Ｆｉｋは、次式により算出される。
【００７３】
【数１６】
2015094751.tif 000018
【００７４】
ここで、Δτｋ←ｉ,ｐは、光子ｐの経路に沿った、体積要素ｋの内部での光学的厚さ
である。Ｎは、面積要素ｉから射出された光子数であり、ｎは、体積要素ｋに入射された光子数である。
【００７５】
体積要素ｋから面積要素ｉをみた形態係数Ｆｋｉの算出時には、体積要素の表面上の各点からの光子の等方的な射出が仮定される。すなわち、体積要素の表面上の各点からLambertの余弦則に従って多数の光子が仮想的に射出される。
【００７６】
そして、体積要素ｋから面積要素ｉをみた形態係数Ｆｋｉが、次式により算出される。
【００７７】
【数１７】
2015094751.tif 000019
【００７８】
ここで、Δτｋ←ｉ,ｐは、光子ｐの射出点から光子ｐの進行方向とは逆の方向へバッ
クトレースすることによって計算される，体積要素ｋの内部での光学的厚さである。また、Ｓｋは、体積要素ｋの表面積であり、Ｍは、体積要素ｋから射出された光子数であり、ｍは、面積要素ｉに入射された光子数である。
【００７９】
同様に、体積要素ｋから体積要素ｌをみた形態係数Ｆｋｌは、次式により計算される。
【００８０】
【数１８】
2015094751.tif 000020
【００８１】
尚、この式におけるｍは、Ｍ個の光子を体積要素ｋから射出した結果として体積要素ｌに入射された光子数である。
【００８２】
以下、パラメータ算出処理で算出される、各体積要素ｋの有効表面積＜Aeff＞ｋ、各要素ｉの天空率ωｉ、各面積要素ｉの日陰率Ｄｉ、各体積要素ｋの有効日陰率Ｄeffkについて説明する。
【００８３】
パラメータ算出処理時に算出される各体積要素ｋの有効表面積＜Aeff＞ｋは、以下の式により定義される値である。
【００８４】
【数１９】
2015094751.tif 000021
【００８５】
この式において、ｍは、体積要素ｋの周囲に存在する（体積要素ｋから見える）要素（面積要素又は体積要素）の総数であり、ｉ（＝１〜ｍ）は、体積要素ｋの周囲に存在する面積要素又は体積要素の要素番号である。
【００８６】
パラメータ算出処理時には、この有効表面積＜Ａｅｆｆ＞ｋが、次式により算出される。
【００８７】
【数２０】
2015094751.tif 000022
【００８８】
すなわち、有効表面積＜Ａｅｆｆ＞ｋもモンテカルロ法により算出される。
【００８９】
要素（面積要素／体積要素）ｉの天空率ωｉは、要素ｉから空をみた形態係数に相当する値である。天空率ωｉは、要素ｉから面積要素をみた形態係数と同様の手順で算出される。
【００９０】
面積要素ｉの日陰率Ｄｉは、面積要素iから太陽側へ射出された光子ｐが、他の要素に
入射した際に失うエネルギーΔＷｐを積算することにより求められる。より具体的には、日陰率Ｄｉは、モンテカルロ法を用いて次式により算出される。尚、次式において、Ｎは、面積要素iから射出された光子数である。
【００９１】
【数２１】
2015094751.tif 000023
【００９２】
同様に、体積要素ｋの有効日陰率Ｄeffkは、モンテカルロ法を用いて以下の式により算出される。
【００９３】
【数２２】
2015094751.tif 000024
【００９４】
この有効日陰率Ｄeffk、上記した日陰率Ｄiは、太陽の位置により値が変わるパラメー
タである。そのため、パラメータ算出処理時には、シミュレーション対象時間範囲内の各シミュレーション時刻における日陰率Ｄｉ及び有効日陰率Ｄeffkが算出される。
【００９５】
《主処理》
以下、演算部１２が行う主処理の内容を説明する。
【００９６】
図２に模式的に示してあるように、主処理は、放射熱算出処理（ステップＳ２０１）と温度算出処理（ステップＳ２０２）と風速分布等算出処理（ステップＳ２０３）とが繰り返される処理である。尚、主処理時における放射熱算出処理等の繰り返し回数は、計算条件データの要素として入力されるシミュレーション対象時間範囲とミュレーション時間間隔Δｔとに基づき決定される。
【００９７】
・放射熱算出処理
ステップＳ２０１にて行われる放射熱算出処理は、
（１）パラメータ算出処理で算出されたパラメータ（形態係数等）と、各部の表面温度の
ｎ回目の算出結果とを用いて、各要素ｉから射出される短波放射（可視光線）の放射フラックスＧS,i[W/m2]及び長波放射（赤外線）の放射フラックスＧL,i [W/m2] を算出し、
（２）算出した放射フラックスと、各部の表面温度のｎ回目の算出結果とを用いて、各要素ｉに吸収される、短波放射に関する正味の放射熱ＲS,i[W]及び長波放射に関する正味の放射熱ＲL,i[W]を算出する
処理である。
【００９８】
尚、各部の表面温度のｎ回目の算出結果とは、ｎ＝０の場合は、各平面要素の表面温度の初期値のことであり、ｎ≠０の場合は、温度算出処理によるｎ回目の処理結果（ｎ回目の温度算出処理により得られた各平面要素の表面温度）のことである。
【００９９】
具体的には、各要素ｉから射出される放射フラックスＧS,i、ＧL,i[W/m2] については
、それぞれ、以下の式が成立する。
【０１００】
【数２３】
2015094751.tif 000025
【０１０１】
尚、（８）、（９）式におけるｎは、面積要素及び体積要素の総数である。<Ａeff>ｉ
は、要素ｉが体積要素である場合には、体積要素ｉの有効表面積であり、要素ｉが面積要素である場合には、面積要素ｉの面積である。
【０１０２】
αS,i、αL,iは、それぞれ、要素iの、短波放射、長波放射に関する反射率であり、εiは、体積要素ｉの射出率である。Ｓdirect,iは、要素ｉに入射する太陽からの直達短波放射フラックスであり、Ｓdiffuse,iは、要素iに入射する大気散乱短波放射の放射フラックスである。Ｌｉは、体積要素ｉに入射する大気長波放射の放射フラックスであり、Ａeffi←Solar、Ａeffi←skyは、それぞれ、要素iの、太陽の方向、空の方向の有効面積である
。
【０１０３】
Ｂ(Ｔｉ)は、要素ｉから熱放射により射出される放射フラックスである。長波放射の放射フラックスとして、上記したＧＬ，ｊのみを算出する場合（長波放射に関する放射フラックスを波長範囲別に算出しない場合）、Ｂ(Ｔｉ)は、ステファン・ボルツマン定数σを用いて次式により算出される。
【０１０４】
【数２４】
2015094751.tif 000026
【０１０５】
Ｓdirect,i、Ｓdiffuse,iは、パラメータ算出処理で算出されている天空率ωｉ、今回
、シミュレーション対象となっている時刻用の日陰率Ｄｉを用いて以下の式により算出される。
【０１０６】
【数２５】
2015094751.tif 000027
【０１０７】
ここで，Ｓ↓は水平面に下向きに入射する太陽放射フラックスであり、Ｓ（＝(Sx,Sy,Sz)）は太陽方位ベクトルである。ｎiは、面積要素ｉの単位法線ベクトルであり、cdirect、cdiffuseは、直散分離のための係数である。
【０１０８】
上記式で算出されるＳdirect,iは、面積要素ｉのＳdirect,iである。体積要素ｉのＳdirect,iは，パラメータ算出処理で算出されている，今回、シミュレーション対象となっている時刻用の有効日陰率Ｄeffiを用いて次式により算出される。
【０１０９】
【数２６】
2015094751.tif 000028
【０１１０】
尚、射出率εiは、要素ｉの吸収率と等しいため、射出率εiとして、“１−αＬ,ｉ”
を使用することが出来る。
【数２７】
2015094751.tif 000029
【０１１１】
上記した（８）式及び（９）式は、ｉ＝１〜ｎのそれぞれについて成立する。すなわち、以下の２つの線形行列方程式が成立する。
【０１１２】
【数２８】
2015094751.tif 000030
【０１１３】
放射熱算出処理時には、まず、これらの線形行列方程式を解くことで、各要素iから射
出される放射フラックスＧS,i，ＧL,iが、算出される。
【０１１４】
そして、放射熱算出処理時には、各要素ｉに吸収される、短波放射に関する正味の放射熱ＲS,i[W]、及び長波放射に関する正味の放射熱ＲL,i[W]が、それぞれ、以下の（１０）、（１１）式により算出される。
【数２９】
2015094751.tif 000031
【０１１５】
・温度算出処理
温度算出処理（図２；ステップＳ２０２）は、放射熱算出処理にて算出された放射熱に基づく熱収支計算によりシミュレーション対象空間内の各面積要素（建物、地面）の表面温度や体積要素（樹冠内の葉）の表面温度を算出する処理である。図２に模式的に示してあるように、この温度算出処理は、風速分布等算出処理のｎ回目の算出結果を利用するものとなっている。
【０１１６】
温度算出処理時における各面積要素の温度の算出手順は、熱収支計算に使用する放射熱が、樹冠を透過性を有する体積要素として取り扱った形態係数を用いて算出されたものであることを除けば、熱収支計算による温度の一般的な算出手順と同様のものである。そのため、以下では、樹冠（体積要素）の温度の算出手順のみを説明することにする。
【０１１７】
図５に模式的に示してあるように、樹木の樹冠の一部である体積要素には、短波放射の放射熱フラックスＲSと長波放射の放射熱フラックスＲLとが流入し、当該体積要素からは、顕熱フラックスＨと潜熱フラックスＬＥとが流出する。
従って、樹冠である体積要素ｉに関する熱収支は、次式により表されることになる。
【０１１８】
【数３０】
2015094751.tif 000032
【０１１９】
ここで，Ｔleaf,iは、要素ｉ内の葉の表面温度[K]であり、aiは、要素ｉ内の葉の面積
密度[m2/m3]である。Ｖiは、要素ｉの体積[m3]であり、Ｃは、単位葉面積あたりの葉の熱容量[J/K/m2]である。ＲS,i、ＲL,iは、それぞれ、葉に吸収される短波放射の正味の放射熱（放射熱フラックスの強度）[W]、長波放射の正味の放射熱[W]であり、Ｌは、蒸発潜熱[J/kg]である。
【０１２０】
Ｈiは、葉から大気に放出される顕熱輸送量（顕熱フラックスの強度）[W]であり、Ｅi
は葉から大気に蒸散される水蒸気量[kg/s]である。
【０１２１】
葉から大気に放出される顕熱輸送量Ｈiおよび葉から大気に蒸散される水蒸気量Ｅiは、以下の式により算出される。
【０１２２】
【数３１】
2015094751.tif 000033
【０１２３】
ここで、Ｔair,iは、要素ｉ内の大気温度[K]，fair,iは、体積要素ｉ内の大気中の水蒸気分圧[Pa]、ｆleaf,iは、体積要素ｉ内の葉表面の飽和水蒸気分圧[Pa]，Ｋhは、対流熱
伝達係数[W/m2/K]，Ｋwは、対流水蒸気輸送係数[kg/s/m2/Pa]、βは蒸発効率である．
【０１２４】
温度算出処理時には、或る時刻における葉表面温度及び熱フラックスを用いて、時間Δｔ後における葉表面温度Ｔleaf,iが、計算される。
【０１２５】
具体的には、時間Δｔの経過による葉表面温度の変化量ΔＴleaf,iは、Δtの経過に起
因する葉表面温度の変化による正味長波放射，顕熱輸送量，蒸散量の変化を考慮すると、次式により表される。
【０１２６】
【数３２】
2015094751.tif 000034
【０１２７】
従って、葉表面温度の変化量ΔＴleaf,iを、次式により求めることができる。
【０１２８】
【数３３】
2015094751.tif 000035
【０１２９】
温度算出処理時には、この式により、葉表面温度の変化量ΔＴleaf,iが求められてから、時間Δｔ経過後の葉表面温度Ｔleaf,i(n+1)が次式により、算出される。
【０１３０】
【数３４】
2015094751.tif 000036
【０１３１】
・風速分布等算出処理
主処理の要素として行われる風速分布等算出処理（図２；ステップＳ２０３）は、基本的には、密度や運動量の輸送方程式を解くことにより、シミュレーション対象空間内の各時刻における各部の流速、温度、水蒸気量を算出する処理である。また、風速分布等算出処理は、樹冠内の葉の抵抗、葉によるwake productionとdissipationの効果を考慮に入れた形で、各部の流速、温度、水蒸気量を算出する処理となっている。
【０１３２】
ただし、風速分布等算出処理は、樹冠内の葉の抵抗、葉によるwake productionとdissipationの効果を考慮するために、神田らにより提案されている樹木モデル（非特許文献１１参照）を用いたものとなっている。
【０１３３】
そのため、詳細説明は省略するが、風速分布等算出処理時には、以下に記す葉の抵抗を考慮した運動方程式と、葉によるwake productionとdissipationの効果を含む格子スケール以下の乱流エネルギー収支式とが連立して解かれる。
【０１３４】
【数３５】
2015094751.tif 000037
【０１３５】
ここで、Ｆ：樹冠による抗力，Ｃｄ：葉の抵抗係数，ａ：葉面積密度[m2/m3]，Ｕ：ス
カラー風速[m/s]，ρ：密度[g/m3]，p：圧力[Pa]，Ｒij：サブグリッドスケールのレイノルズ応力である。
【０１３６】
【数３６】
2015094751.tif 000038
【０１３７】
ここで，Ｅ：ＳＧＳ（サブグリッドスケール）乱流エネルギー[m2/s2]，Ｓ：樹冠によ
るwake productionおよびエネルギー消散の項， Ｋｍ：渦粘性係数[m2/s]，Δ：乱れの代表長さ[m]，Ｃｓ, Ｃｅ：ＬＥＤ定数，Ｃｗ：葉の抵抗により格子以上の流れ場が失うエ
ネルギーのうち格子スケール以下のwake productionに変換される比率（０〜１）である
。
【０１３８】
主処理を行っている演算部１２は、温度算出処理及び風速分布等算出処理が完了する度に、各処理による算出結果を出力（本実施形態では、記憶部１４への記憶）する。そして、演算部１２は、温度算出処理等の繰り返し回数が、シミュレーション対象時間範囲等から決定した繰り返し回数未満であった場合には、再び、放射熱算出処理、温度算出処理及び風速分布等算出処理を行い、温度算出処理等の繰り返し回数が、シミュレーション対象時間範囲等から決定した繰り返し回数未満となっていた場合に、主処理を終了する。
【０１３９】
以上、説明したように、本実施形態に係るシミュレーション装置１０は、複数の建物及び複数の樹木が存在するシミュレーション対象空間内の各要素が吸収する放射熱の効果、及び、樹木−大気間の熱・水蒸気交換の効果を考慮に入れて当該シミュレーション対象空間における非定常な風速・気温・水蒸気の空間分布を算出する構成を有する。従って、本発明のシミュレーション装置によれば、複数の建物及び複数の樹木が存在するシミュレーション対象空間における風速、温度及び水蒸気量の非定常な空間分布を、高精度に（既存のシミュレーション技術よりも適切に）、算出できることになる。
【０１４０】
また、シミュレーション装置１０は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数（（５）、（６）式参照）として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置１０は、２つの体積要素に関する形態係数（（７）式参照）として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置１０は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、シミュレーション装置１０によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。
【０１４１】
《変形形態》
上記した実施形態に係るシミュレーション装置１０は、各種の変形を行えるものである。例えば、シミュレーション装置１０を、樹冠を体積要素として取り扱わない装置（樹冠表面を面積要素として取り扱って放射熱を算出する装置）に変形することができる。また、シミュレーション装置１０を、２要素間の透過率Ｔが考慮されていない形態係数を算出し、放射フラックスの算出時等に２要素間の透過率を考慮する装置に変形することも出来る。ただし、形態係数の算出時に２要素間の透過率Ｔを考慮しておいた方が正確な結果が得られるし、計算コストも低くなる。従って、上記した形態係数を採用しておくことが好ましい。
【０１４２】
また、樹冠の体積熱容量Ｃａは、通常、かなり小さいので、シミュレーション装置１０を、葉表面温度を、Ｃａ＝０として算出する装置に変形することも出来る。また、シミュレーション装置１０を、一部又は全ての形態係数として、その定義式の解析解を用いる装置に変形することも出来る。
【０１４３】
シミュレーション装置１０は、オイラー陰解法により葉温度を算出する装置であるが、シミュレーション装置１０を、陽解法により葉温度を算出する装置や、２次精度のクランク・ニコルソン法により葉温度を算出する装置に変形することも出来る。ただし、陰解法の方が正確な値を得やすいので、葉温度の算出法には、陰解法を採用しておくことが好ましい。
【符号の説明】
【０１４４】
１０ シミュレーション装置
１２ 演算部
１４ 記憶部
１６ インタフェ−ス（Ｉ／Ｆ）回路
【図１】
2015094751.tif 000039
【図２】
2015094751.tif 000040
【図３】
2015094751.tif 000041
【図４】
2015094751.tif 000042
【図５】
2015094751.tif 000043