(57)【要約】【課題】樹冠を含む３次元空間における放射熱輸送現象を低計算コストで良好にシミュレートできるシミュレーション装置を提供する。【解決手段】シミュレーション装置を、複数の面積要素と複数の体積要素とで定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出手段と、前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出手段とを備え、前記３次元空間内に存在する複数の樹木の樹冠が前記複数の体積要素として取り扱われるように、前記３次元空間が前記複数の面積要素と前記複数の体積要素とにより定義されている装置として構成する。【選択図】図２
2015095239.tif 000002
【特許請求の範囲】
【請求項１】
放射熱輸送現象をシミュレートするシミュレーション装置において、
複数の面積要素と複数の体積要素とで定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出手段と、
前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出手段と、
を備え、
前記３次元空間内に存在する複数の樹木の樹冠が前記複数の体積要素として取り扱われるように、前記３次元空間が前記複数の面積要素と前記複数の体積要素とにより定義されている
ことを特徴とするシミュレーション装置。
【請求項２】
前記形態係数算出手段は、各２要素に関する形態係数として、各２要素間におけう放射熱量の透過率を考慮に入れた形態係数を算出する
ことを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
前記放射熱量算出手段により算出された放射熱量に基づき、各要素の温度を算出する温度算出手段を、さらに備える
ことを特徴とする請求項１又は２に記載のシミュレーション装置。
【請求項４】
前記温度算出手段は、各要素の温度を陰解法により算出する
ことを特徴とする請求項３に記載のシミュレーション装置。
【請求項５】
放射熱輸送現象をシミュレートするシミュレーション方法において、
コンピュータが、
複数の樹木の樹冠が複数の体積要素として取り扱われるように、複数の面積要素と複数の体積要素とが定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出ステップと、
前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
を行うことを特徴とするシミュレーション方法。
【請求項６】
放射熱輸送現象をシミュレートするためのシミュレーションプログラムにおいて、
コンピュータに、
複数の樹木の樹冠が複数の体積要素として取り扱われるように、複数の面積要素と複数の体積要素とが定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出ステップと、
前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
を行わせることを特徴とするシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、放射熱輸送現象をシミュレートするためのシミュレーション装置とシミュレーション方法とシミュレーションプログラムとに関する。
【背景技術】
【０００２】
都市部の気温がその周辺部の気温より高くなる現象を、ヒートアイランド現象という。
ヒートアイランド現象を緩和するための対策として、街路樹の配置や緑地の整備が着目されている。そのため、街路樹の配置等を考慮して都市部の放射熱輸送をシミュレートする様々なシミュレーション技術（特許文献１及び非特許文献１〜９参照）が開発されている。
【０００３】
ただし、既存のシミュレーション技術は、計算コスト（計算量）を抑えるために、『樹冠による放射熱の散乱を無視し、樹冠による放射熱の減衰・吸収のみを考慮する』といったように、樹冠に関する放射エネルギーの伝達過程を単純化したものとなっている。そのため、既存のシミュレーション技術では、樹冠付近や樹冠内部での放射エネルギーの吸収量等を適切にシミュレートできないことが多いのが現状である。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００３−０９９６９７号公報
【特許文献２】特許５１３７０３９号公報
【特許文献３】特開２０１２−０２１６８４号公報
【非特許文献】
【０００５】
【非特許文献１】吉田 伸治, 村上 周三, 持田 灯, 大岡 龍三, 富永 禎秀, 「環境緩和効果を総合的に組み込んだ新しい3次元樹木モデルの開発」, 生産研究, 51巻1号, 1999年1月
【非特許文献２】吉田 伸治, 大岡 龍三, 持田 灯, 富永 禎秀, 村上 周三, 「樹木モデルを組み込んだ対流・放射・湿気輸送連成解析による樹木の屋外温熱環境緩和効果の検討」, 日本建築学会計画系論文集, No.536, pp.87-94, 2000
【非特許文献３】坂本 雄三, 小島 悦史, 足永 靖信, 今野 雅, 「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
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明の課題は、樹冠を含む３次元空間における放射熱輸送現象を低計算コストで良好にシミュレートできる技術を提供することにある。
【課題を解決するための手段】
【０００７】
上記課題を解決するために、本発明の、放射熱輸送現象をシミュレートするシミュレーション装置は、
複数の面積要素と複数の体積要素とで定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出手段と、
前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出手段と、
を備え、
前記３次元空間内に存在する複数の樹木の樹冠が前記複数の体積要素として取り扱われるように、前記３次元空間が前記複数の面積要素と前記複数の体積要素とにより定義されている
装置として構成される。
【０００８】
すなわち、本発明のシミュレーション装置は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置は、２つの体積要素に関する形態係数として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、本発明のシミュレーション装置によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。
【０００９】
本発明は、上記シミュレーション装置と同様の特徴を有するシミュレーション方法や、情報処理装置（コンピュータ）を上記シミュレーション装置として機能させるシミュレーションプログラムとしても実現することが出来る。尚、本発明を、上記シミュレーションプログラムを記録したコンピュータ可読媒体としても実現することも出来る。
【発明の効果】
【００１０】
本発明によれば、樹冠を含む３次元空間における放射熱輸送現象を低計算コストで良好にシミュレートすることができる。
【図面の簡単な説明】
【００１１】
【図１】本発明の一実施形態に係るシミュレーションシステムの構成図
【図２】実施形態に係るシミュレーション装置の機能の説明図
【図３】２面積要素に関する形態係数についての算出式中のパラメータの説明図
【図４】面積要素から体積要素をみた形態係数についての算出式中のパラメータの説明図
【図５】体積要素の熱収支の説明図
【発明を実施するための形態】
【００１２】
以下、図面を参照して、本発明の一実施形態について説明する。尚、以下の実施形態の説明は、例示であり、本発明は、実施形態の構成に限定されない。
【００１３】
図１に、実施形態に係るシミュレーションシステムの構成を示す。図示してあるように、本実施形態に係るシミュレーションシステムは、シミュレーション装置１０、入力装置２０及び出力装置３０を含む。また、シミュレーション装置１０は、演算部１２、記憶部１４、インタフェース回路（Ｉ／Ｆ）１６を備える。
【００１４】
シミュレーション装置１０のインタフェース回路（Ｉ／Ｆ）１６は、演算部１２が他装置と通信を行うために使用する回路である。記憶部１４は、シミュレーションプログラム１８を記憶した不揮発性の記憶装置である。この記憶部１４は、演算部１２が利用するデータや演算部１２による処理結果の記憶にも使用される。
【００１５】
演算部１２は、ＣＰＵ(Central Processing Unit)、ＲＡＭ(Random Access Memory)等
を組み合わせたユニットである。演算部１２は、シミュレーションプログラム１８を記憶部１４から読み出して実行することにより、各種処理（詳細は後述）を行う。また、演算部１２は、シミュレーションプログラム１８を実行することにより、本発明に係る形態係数算出手段、放射熱量算出手段及び温度算出手段として機能する。
【００１６】
入力装置２０は、シミュレーション装置１０に情報を入力するための装置である。入力装置２０は、キーボード、マウス等のポインティングデバイス等を含む。また、出力装置３０は、シミュレーション装置１０からの情報を出力するための、ＬＣＤ(Liquid Crystal Display)やＣＲＴ(Cathode-Ray Tube)等のディスプレイ、プリンタ等である。
【００１７】
尚、シミュレーション装置１０は、通常、行列演算が高速に行えるコンピュータ（ベクトル型並列計算機等）にシミュレーションプログラム１８を実行させることによって実現される。従って、通常、シミュレーション装置１０（ベクトル型並列計算機等）に接続された入出力用のコンピュータの入力装置及び出力装置が、入力装置２０及び出力装置３０として機能する。
【００１８】
以下、シミュレーション装置１０の機能を説明する。
シミュレーション装置１０は、複数の樹木を含む都市空間（以下、シミュレーション対象空間と表記する）内での放射熱輸送現象をシミュレートするための装置である。
【００１９】
図２に模式的に示してあるように、シミュレーション装置１０の使用時には、シミュレーション対象時間範囲内の太陽の高度や方位に関する情報等が設定された計算条件設定ファイルが記憶部１４に記憶される。また、シミュレーション装置１０に対して、地形データ、建物形状データ、樹木分布データ、地表面温度データ、建物表面温度データ及び葉表面温度データが入力される。
【００２０】
地形データは、シミュレーション対象空間の地形（地面の形状）を示すデータである。建物形状データは、シミュレーション対象空間内に存在する各建物の位置及び形状を示すデータである。樹木分布データは、シミュレーション対象空間内に存在する各樹木の位置及び形状と各樹木の葉面積密度分布とを示すデータである。
【００２１】
シミュレーション装置１０に入力されるこれらのデータは、シミュレーション対象空間の構造（３次元都市形状及び３次元樹木分布）が分かるデータであれば良い。従って、地形データや建物形状データとして、例えば、地面／建物の高さの２次元平面データ（地面
／建物の高さと座標との対応関係を表すデータ）を使用することができる。また、樹木分布データとして、例えば、樹木インデックス（樹木の識別情報）の２次元平面データと、各樹木インデックスで識別される樹木の葉面積密度の鉛直分布を示すデータとを含むデータを使用することができる。
【００２２】
シミュレーション装置１０に入力される地表面温度データ、建物表面温度データ、葉表面温度データは、それぞれ、地表面の各部の温度の初期値を示すデータ、各建物表面の各部の温度の初期値を示すデータ、各所の葉表面の温度の初期値を示すデータである。
【００２３】
以下、シミュレーション装置１０（演算部１２）が行う処理の内容を説明する。尚、以下の説明では、シミュレーション装置１０に入力される地形データ、建物形状データ及び樹木分布データのことを、構造指定データと表記する。また、シミュレーション装置１０に入力される地表面温度データ、建物表面温度データ及び葉表面温度データのことを、初期値データと表記する。
【００２４】
シミュレーション装置１０は、基本的には、入力された構造指定データが示しているシミュレーション対象空間内の各部の温度を、計算条件設定ファイル内の各種情報及び入力された初期値データを用いて、Δｔ毎にシミュレートする装置である。
【００２５】
図２に示してあるように、シミュレーション装置１０の演算部１２が行う処理は、前処理と主処理とに分類できる。
【００２６】
前処理は、シミュレーション用データ生成処理（ステップＳ１０１）とパラメータ算出処理（ステップＳ１０２）とがこの順に行われる処理である。
【００２７】
ステップＳ１０１にて行われるシミュレーション用データ生成処理は、入力された構造指定データに基づき、『シミュレーション対象空間内の各建物の壁面や地面を複数の面積要素として取り扱い、シミュレーション対象空間内の各樹木の樹冠を１つ以上の透過性を有する体積要素として取り扱うためのシミュレーション用データ』を生成する処理である。
【００２８】
このシミュレーション用データ生成処理により生成されるシミュレーション用データは、後述する形態係数の算出に必要な情報（各面積／体積要素の形状やシミュレーション対象空間内での位置等）が分かるものであれば良い。従って、シミュレーション用データを、例えば、面積／体積要素毎に、『通し番号、対応するシミュレーション対象空間内の計算格子の座標番号、面積要素の向いている方向／体積要素であることを示すフラグ、樹冠であるか否かを示すフラグ等からなるデータ』を含むデータとしておくことができる。
【００２９】
ステップＳ１０２にて行われるパラメータ算出処理は、シミュレーション用データ生成処理により生成されたシミュレーション用データに基づき、主処理で使用される各種パラメータを算出する処理である。
【００３０】
このパラメータ算出処理時に算出されるパラメータには、各２要素（面積／体積要素）に関する形態係数Ｆ、各体積要素ｋの有効表面積＜Aeff＞k、各要素ｉの天空率ωi、各面積要素ｉの日陰率Ｄi、各体積要素ｋの有効日陰率Ｄeffkがある。
【００３１】
まず、パラメータ算出処理時に算出される各２要素に関する形態係数Ｆについて説明する。
【００３２】
パラメータ算出処理時には、２つの面積要素ｉ，ｊの組み合わせ毎に、『面積要素ｉか
ら面積要素ｊを見た形態係数Ｆij』と『面積要素ｊから面積要素ｉを見た形態係数Ｆji』とが、算出される。また、面積要素ｉと体積要素ｋの組み合わせ毎に、『面積要素ｉから体積要素ｋを見た形態係数Ｆik』と『体積要素ｋから面積要素ｉを見た形態係数Ｆki』とが、算出される。さらに、２つの体積要素ｋ及びｌの組み合わせ毎に、『体積要素ｋから体積要素ｌを見た形態係数Ｆkl』及び『体積要素ｌから体積要素ｋを見た形態係数Ｆlk』が算出される。
【００３３】
パラメータ算出処理時に算出される『面積要素ｉから面積要素ｊを見た形態係数Ｆij』、『面積要素ｊから面積要素ｉを見た形態係数Ｆji』は、それぞれ、以下の（１）、（２）式で定義される値である。
【００３４】
【数１】
2015095239.tif 000003
【００３５】
これらの式において、Ａi、Ａjは、それぞれ、面積要素ｉの面積、面積要素ｊの面積である。βi、βj は、図３に模式的に示してあるように、それぞれ、面積要素ｉ内の微小
面ｄＡiと面積要素ｊ内の微小面ｄＡjとを結ぶ直線と微小面ｄＡiの法線のなす角度、当
該直線と微小面ｄＡjの法線のなす角度である。また、ｒは、微小面ｄＡi・微小面ｄＡj
間の距離である。
【００３６】
Ｔijは、面積要素ｉと面積要素ｊの間の透過率である。Ｔijは、２微小面ｄＡi及びｄ
Ａj 間の光学的厚さτijを用いて次式により算出される。
【００３７】
【数２】
2015095239.tif 000004
【００３８】
また、面積要素ｉと面積要素ｊとの間（微小面ｄＡiと微小面ｄＡjとの間）に樹冠が分布している場合、光学的厚さτijは、樹冠の消散係数ｋext及び葉面積密度ａを用いて次
式により算出される。
【００３９】
【数３】
2015095239.tif 000005
【００４０】
尚、形態係数の具体的な算出手順については後述するが、上記した（１）、（２）式で定義される形態係数は、相反則を満たしている。すなわち、面積Ａiと面積Ａjと形態係数Ｆijと形態係数Ｆjiとの間には、以下の関係がある。
【００４１】
【数４】
2015095239.tif 000006
【００４２】
従って、Ｆjiを、（１）式から算出したＦijとＡiとＡjとから算出することも、Ｆijを、（２）式から算出したＦjiとＡiとＡjとから算出することもできる。
【００４３】
パラメータ算出処理時に算出される『面積要素ｉから体積要素ｋを見た形態係数Ｆik』は、以下の（３）式により定義される値である。
【００４４】
【数５】
2015095239.tif 000007
【００４５】
この（３）式におけるβiは、図４に模式的に示してあるように、面積要素ｉの微小面
ｄＡiと体積要素ｊ内の微小投影面ｄＡk←iとを結ぶ直線と微小面ｄＡiの法線のなす角度である。また、ｒは、微小面ｄＡi・微小投影面ｄＡk←i間の距離である。
【００４６】
Ａeffk←iは、体積要素ｋ自身の遮蔽率を考慮した，面積要素ｉの方向から見た体積要
素ｋの有効面積である。このＡeffk←iは、次式により算出される。
【００４７】
【数６】
2015095239.tif 000008
【００４８】
ここで、Δτk←iは、微小投影面ｄＡk←i(図４参照)に垂直な方向の体積要素ｋの光学的厚さである。体積要素ｋが樹木の場合、Δτk←iは、消散係数ｋext、葉面積密度ａお
よびｄＡk←iに垂直な方向の体積要素ｋの幾何学的厚さΔｓk←iを用いて、次式により算出される。
【００４９】
【数７】
2015095239.tif 000009
【００５０】
要するに、上記した（３）式は、（４）式を用いれば、以下のように表記できる式となっている。パラメータ処理時には、以下の（５）式に従って、面積要素ｉから体積要素ｋを見た形態係数Ｆikが算出される。
【００５１】
【数８】
2015095239.tif 000010
【００５２】
パラメータ算出処理時に算出される『体積要素ｋから面積要素ｉを見た形態係数Ｆki』
は、次式により定義される値である。
【００５３】
【数９】
2015095239.tif 000011
【００５４】
換言すれば、体積要素ｋから面積要素ｉを見た形態係数Ｆkiとしては、以下の式で表される相反則を満たす値が算出される。
【００５５】
【数１０】
2015095239.tif 000012
【００５６】
パラメータ算出処理時に算出される『体積要素ｋから体積要素ｌを見た形態係数Ｆkl』は、次式により定義される値である。
【００５７】
【数１１】
2015095239.tif 000013
【００５８】
すなわち、この形態係数Ｆklは、体積要素ｋの遮蔽率（“1-exp(-Δτk←l)”）と体積要素ｌの遮蔽率（“1-exp(-Δτl←k)”）と体積要素ｋ、ｌ間の透過率Ｔklとを考慮に入れた形態係数となっている。
【００５９】
尚、（７）式が表している体積要素に関する形態係数は、面積要素の形態係数と同様に次式の相反則を満たす。
【００６０】
【数１２】
2015095239.tif 000014
【００６１】
従って、他の形態係数と同様に、２体積要素間の形態係数の双方を、（７）式に従って算出することも、一方の形態係数を（７）式に従って算出し、他方の形態係数を、当該一方の形態係数の算出結果から算出することも出来る。
【００６２】
パラメータ算出処理時には、上記した各形態係数が、モンテカルロ法により算出される。
【００６３】
すなわち、形態係数の算出時には、０から１までの範囲の一様乱数Ｒθ、Ｒφから、以下の式によりμ、φを算出し、算出結果に基づき、以下に示した内容の単位ベクトルｎを生成する処理が、必要とされる形態係数の精度に応じた回数、繰り返される。
【００６４】
【数１３】
2015095239.tif 000015
【００６５】
換言すれば、形態係数の算出時には、Lambertの余弦則にしたがって面積要素（又は体
積要素）から射出される多数の光子の進行方向単位ベクトルｎを生成する処理が行われる。
【００６６】
そして、生成した各単位ベクトルｎの方向に要素ｉから等しいエネルギーＷ0を持った
光子が射出された場合に要素ｊに入射されることになる各光子ｐのエネルギーＷpを積算
することにより、形態係数Ｆijが求められる。
【００６７】
具体的には、面積要素ｉから面積要素ｊをみた形態係数Ｆijは、面積要素ｉから射出させた光子の個数がＮ、面積要素ｊに入射された光子の個数がｎであった場合には、以下の式により算出される。
【００６８】
【数１４】
2015095239.tif 000016
【００６９】
なお、建物壁面などの完全遮蔽性の面積要素のみの場合には、面積要素ｊに入射する光子のエネルギーＷpはＷ0と等しいため、形態係数は、Ｆij＝ｎ／Ｎにより算出できる。樹木のような放射透過性の要素が空間中に分布する場合には、面積要素ｊに入射する光子のエネルギーＷpは、面積要素ｉから面積要素ｊに到達するまでの間に減衰する。そのため
、その減衰の形態係数に対する影響を光子の持つエネルギーを減衰させることで考慮する。すなわち、次式により、Ｗpが計算される。
【００７０】
【数１５】
2015095239.tif 000017
【００７１】
ここで、Ｔij,pは、光子ｐの経路に沿った透過率である。なお、光子が透過性の体積要素内で確率的に遮蔽されると考えて、面積要素ｊに入射する光子の数を減少させても、減衰の影響を考慮した形態係数を求めることが出来る。
【００７２】
また、面積要素ｉから体積要素ｋをみた形態係数Ｆikは、（７）式により定義される値である。従って、面積要素ｉから体積要素ｋをみた形態係数Ｆikは、次式により算出される。
【００７３】
【数１６】
2015095239.tif 000018
【００７４】
ここで、Δτk←i,pは、光子ｐの経路に沿った、体積要素ｋの内部での光学的
厚さである。Ｎは、面積要素ｉから射出された光子数であり、ｎは、体積要素ｋに入射された光子数である。
【００７５】
体積要素ｋから面積要素ｉをみた形態係数Ｆkiの算出時には、体積要素の表面上の各点からの光子の等方的な射出が仮定される。すなわち、体積要素の表面上の各点からLambertの余弦則に従って多数の光子が仮想的に射出される。
【００７６】
そして、体積要素ｋから面積要素ｉをみた形態係数Ｆkiが、次式により算出される。
【００７７】
【数１７】
2015095239.tif 000019
【００７８】
ここで、Δτk←i,pは、光子ｐの射出点から光子ｐの進行方向とは逆の方向へ
バックトレースすることによって計算される，体積要素ｋの内部での光学的厚さである。また、Ｓkは、体積要素ｋの表面積であり、Ｍは、体積要素ｋから射出された光子数であ
り、ｍは、面積要素ｉに入射された光子数である。
【００７９】
同様に、体積要素ｋから体積要素ｌをみた形態係数Ｆklは、次式により計算される。
【００８０】
【数１８】
2015095239.tif 000020
【００８１】
尚、この式におけるｍは、Ｍ個の光子を体積要素ｋから射出した結果として体積要素ｌに入射された光子数である。
【００８２】
以下、パラメータ算出処理で算出される、各体積要素ｋの有効表面積＜Aeff＞k、各要
素ｉの天空率ωi、各面積要素ｉの日陰率Ｄi、各体積要素ｋの有効日陰率Ｄeffkについて説明する。
【００８３】
パラメータ算出処理時に算出される各体積要素ｋの有効表面積＜Aeff＞kは、以下の式
により定義される値である。
【００８４】
【数１９】
2015095239.tif 000021
この式において、ｍは、体積要素ｋの周囲に存在する（体積要素ｋから見える）要素（面積要素又は体積要素）の総数であり、ｉ（＝１〜ｍ）は、体積要素ｋの周囲に存在する面積要素又は体積要素の要素番号である。
【００８５】
パラメータ算出処理時には、この有効表面積＜Ａeff＞kが、次式により算出される。
【００８６】
【数２０】
2015095239.tif 000022
【００８７】
すなわち、有効表面積＜Ａeff＞kもモンテカルロ法により算出される。
【００８８】
要素（面積要素／体積要素）ｉの天空率ωiは、要素ｉから空をみた形態係数に相当す
る値である。天空率ωiは、要素ｉから面積要素をみた形態係数と同様の手順で算出され
る。
【００８９】
面積要素ｉの日陰率Ｄiは、面積要素iから太陽側へ射出された光子ｐが、他の要素に入射した際に失うエネルギーΔＷｐを積算することにより求められる。より具体的には、日陰率Ｄiは、モンテカルロ法を用いて以下の式により算出される。
【００９０】
【数２１】
2015095239.tif 000023
【００９１】
尚、上記式において、Ｎは、面積要素iから射出された光子数である。
【００９２】
同様に、体積要素ｋの有効日陰率Ｄeffkは、モンテカルロ法を用いて以下の式により算出される。
【００９３】
【数２２】
2015095239.tif 000024
【００９４】
上記した各面積要素ｉの日陰率Ｄi、各体積要素ｋの有効日陰率Ｄeffkは、太陽の位置
により値が変わるパラメータである。従って、シミュレーションを行う時間範囲が比較的に長い場合には、パラメータ算出処理時に、当該時間範囲内の各時刻における日陰率Ｄi
及び有効日陰率Ｄeffkが算出される。
【００９５】
以下、演算部１２が行う主処理（図２）の内容を説明する。
演算部１２が行う主処理は、放射熱フラックス等算出処理（ステップＳ２０１）と温度算出処理（ステップＳ２０２）とが、全時間ステップ数Ntの回数、繰り返される処理である。尚、全時間ステップ数Ntは、シミュレーション対象時刻範囲及び時間刻み幅Δtに基
づいて決定することができる。全時間ステップ数Ntの指定は、計算条件設定ファイル内に全時間ステップ数Nt又はシミュレーション対象時刻範囲及び時間刻み幅Δtの情報を設定
しておくことや、入力装置２０を用いて入力することによって、行われる。
【００９６】
ステップＳ２０１にて行われる放射熱フラックス等算出処理は、パラメータ算出処理で算出されたパラメータ（形態係数等）を用いて、各要素ｉから射出される短波放射（可視光線）の放射フラックスＧS,i[W/m2]及び長波放射（赤外線）の放射フラックスＧL,i [W/m2] を算出し、算出した放射フラックスを用いて、各要素ｉに吸収される、短波放射に関する正味の放射熱ＲS,i[W]及び長波放射に関する正味の放射熱ＲL,i[W]を算出する処理である。
【００９７】
具体的には、各要素ｉから射出される放射フラックスＧS,i、ＧL,i [W/m2] については、それぞれ、以下の式が成立する。
【００９８】
【数２３】
2015095239.tif 000025
【００９９】
尚、（８）、（９）式におけるｎは、面積要素及び体積要素の総数である。<Ａeff>ｉ
は、要素ｉが体積要素である場合には、体積要素ｉの有効表面積であり、要素ｉが面積要素である場合には、面積要素ｉの面積である。
【０１００】
αS,i、αL,iは、それぞれ、要素iの、短波放射、長波放射に関する反射率であり、εiは、要素ｉの射出率である。Ｓdirect,iは、要素ｉに入射する太陽からの直達短波放射フラックスであり、Ｓdiffuse,iは、要素iに入射する大気散乱短波放射の放射フラックスである。Ｌiは、要素ｉに入射する大気長波放射の放射フラックスであり、Ａeffi←Solar、Ａeffi←skyは、それぞれ、要素iの、太陽の方向、空の方向の有効面積である。
【０１０１】
Ｂ(Ｔi)は、要素ｉから熱放射により射出される放射フラックスである。長波放射の放
射フラックスとして、上記したＧL,jのみを算出する場合（長波放射に関する放射フラッ
クスを波長範囲別に算出しない場合）、Ｂ(Ｔi)は、ステファン・ボルツマン定数σを用
いて次式により算出される。
【０１０２】
【数２４】
2015095239.tif 000026
【０１０３】
Ｓdirect,i、Ｓdiffuse,iは、パラメータ算出処理で算出された天空率ωi、日陰率Ｄi
を用いて以下の式により算出される。
【０１０４】
【数２５】
2015095239.tif 000027
【０１０５】
ここで，Ｓ↓は水平面に下向きに入射する太陽放射フラックスであり、Ｓ（＝(Sx,Sy,Sz)）は太陽方位ベクトルである。ｎiは、面積要素ｉの単位法線ベクトルであり、cdirect、cdiffuseは、直散分離のための係数である。
【０１０６】
尚、上記式で算出されるＳdirect,iは、面積要素ｉのＳdirect,iである。体積要素ｉのＳdirect,iは，パラメータ算出処理で算出された有効日陰率Ｄeffiを用いて次式により算出される。
【０１０７】
【数２６】
2015095239.tif 000028
【０１０８】
また、射出率εiは、要素ｉの吸収率と等しいため、射出率εiとして、“１−αL,ｉ”を使用することが出来る。
【０１０９】
【数２７】
2015095239.tif 000029
【０１１０】
上記した（８）式及び（９）式は、ｉ＝１〜ｎのそれぞれについて成立する。すなわち、以下の２つの線形行列方程式が成立する。
【０１１１】
【数２８】
2015095239.tif 000030
【０１１２】
放射フラックス等算出処理時には、これらの線形行列方程式を解くことで、各要素iか
ら射出される放射フラックスＧS,i，ＧL,iが、算出される。
【０１１３】
そして、各要素ｉに吸収される、短波放射に関する正味の放射熱ＲS,i[W]、及び長波放射に関する正味の放射熱ＲL,i[W]が、それぞれ、以下の（１０）、（１１）式により算出される。
【０１１４】
【数２９】
2015095239.tif 000031
【０１１５】
温度算出処理（図２；ステップＳ２０２）は、放射熱フラックス等算出処理にて算出された放射熱からシミュレーション対象空間内の各部の温度を算出する処理である。この温度算出処理時における各面積要素の温度の算出手順は、放射熱が、樹冠を透過性を有する体積要素として取り扱った形態係数を用いて算出されたものであることを除けば、一般的な算出手順と同様のものである。そのため、以下では、樹冠（体積要素）の温度の算出手順のみを説明することにする。
【０１１６】
図５に模式的に示してあるように、樹木の樹冠の一部である体積要素には、短波放射の放射熱フラックスＲSと長波放射の放射熱フラックスＲLとが流入し、当該体積要素からは、顕熱フラックスＨと潜熱フラックスＬＥとが流出する。
従って、樹冠である体積要素ｉに関する熱収支は、次式により表されることになる。
【０１１７】
【数３０】
2015095239.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は、以下の式により算出される。
【０１２１】
【数３１】
2015095239.tif 000033
【０１２２】
ここで、Ｔair,iは、要素ｉ内の大気温度[K]，fair,iは、体積要素ｉ内の大気中の水蒸気分圧[Pa]、ｆleaf,iは、体積要素ｉ内の葉表面の飽和水蒸気分圧[Pa]，Ｋhは、対流熱
伝達係数[W/m2/K]，Ｋwは、対流水蒸気輸送係数[kg/s/m2/Pa]、βは蒸発効率である．
【０１２３】
温度算出処理時には、時間ステップｎにおける葉表面温度及び熱フラックスを用いて、時間刻みΔt後の時間ステップｎ＋１における葉表面温度Ｔleaf,iが、計算され
る。
【０１２４】
具体的には、時間ステップｎからｎ＋１までの間の葉表面温度の変化量ΔＴleaf,i
は，時間Δtの経過に起因する葉表面温度の変化による正味長波放射，顕熱輸送量，
蒸散量の変化を考慮すると、次式により表される。
【０１２５】
【数３２】
2015095239.tif 000034
【０１２６】
従って，葉表面温度の変化量ΔＴleaf,iを、次式により求めることができる。
【０１２７】
【数３３】
2015095239.tif 000035
【０１２８】
温度算出処理時には、この式により、葉表面温度の変化量ΔＴleaf,iが求められて
から、時間ステップｎ+１における葉表面温度Ｔleaf,i(n+1)が次式により、算出される。
【０１２９】
【数３４】
2015095239.tif 000036
【０１３０】
温度算出処理は、各部の温度の算出と、算出した各部の温度の出力（本実施形態では、記憶部１４への記憶）とが完了したときに終了する。そして、指定された回数の処理が完了していなかった場合には、再び、放射熱フラックス等算出処理が開始され、指定された回数の処理が完了したときに、主処理が終了する。
【０１３１】
以上、説明したように、本実施形態に係るシミュレーション装置１０は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数（（５）、（６）式参照）として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置１０は、２つの体積要素に関する形態係数（（７）式参照）として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置１０は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、シミュレーション装置１０によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。
【０１３２】
《変形形態》
上記した実施形態に係るシミュレーション装置１０は、各種の変形を行えるものである。
例えば、シミュレーション装置１０を、２要素間の透過率Ｔが考慮されていない形態係数を算出し、放射フラックスの算出時等に２要素間の透過率を考慮する装置に変形することが出来る。ただし、形態係数の算出時に２要素間の透過率Ｔを考慮しておいた方が正確な結果が得られるし、計算コストも低くなる。従って、上記した形態係数を採用しておくことが好ましい。
【０１３３】
また、樹幹の体積熱容量Caは、通常、かなり小さいので、シミュレーション装置１０を、葉表面温度を、Ca=0として算出する装置に変形することも出来る。また、シミュレーション装置１０を、一部又は全ての形態係数として、その定義式の解析解を用いる装置に変形することも出来る。
【０１３４】
シミュレーション装置１０は、オイラー陰解法により葉温度を算出する装置であるが、シミュレーション装置１０を、陽解法により葉温度を算出する装置や、２次精度のクランク・ニコルソン法により葉温度を算出する装置に変形することも出来る。ただし、陰解法の方が正確な値を得やすいので、葉温度の算出法には、陰解法を採用しておくことが好ましい。
【符号の説明】
【０１３５】
１０ シミュレーション装置
１２ 演算部
１４ 記憶部
１６ インタフェ−ス（Ｉ／Ｆ）回路
【図１】
2015095239.tif 000037
【図２】
2015095239.tif 000038
【図３】
2015095239.tif 000039
【図４】
2015095239.tif 000040
【図５】
2015095239.tif 000041