特許 6238439 放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6238439

(24)【登録日】2017年11月10日

(45)【発行日】2017年11月29日

(54)【発明の名称】放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム

(51)【国際特許分類】

   G06F 19/00 20110101AFI20171120BHJP

【ＦＩ】

   G06F19/00 110

【請求項の数】6

【全頁数】20

(21)【出願番号】特願2013-236190(P2013-236190)

(22)【出願日】2013年11月14日

(65)【公開番号】特開2015-95239(P2015-95239A)

(43)【公開日】2015年5月18日

【審査請求日】2016年11月11日

(73)【特許権者】

【識別番号】504194878

【氏名又は名称】国立研究開発法人海洋研究開発機構

(74)【代理人】

【識別番号】100100549

【弁理士】

【氏名又は名称】川口 嘉之

(74)【代理人】

【識別番号】100113608

【弁理士】

【氏名又は名称】平川 明

(74)【代理人】

【識別番号】100105407

【弁理士】

【氏名又は名称】高田 大輔

(72)【発明者】

【氏名】松田 景吾

(72)【発明者】

【氏名】高橋 桂子

【審査官】 宮地 匡人

(56)【参考文献】

【文献】 国際公開第２００９／０５７１８９（ＷＯ，Ａ１）

【文献】 吉田伸治，環境緩和効果を総合的に組み込んだ新しい３次元樹木モデルの開発，生産研究，１９９９年 １月，Vol.51 No.1，pp.11-16

【文献】 萩島理，数値計算による街路樹の暑熱緩和効果の評価，日本建築学会計画系論文集，１９９９年１１月３０日，No.525，pp.83-90

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０６Ｆ １９／００

ＪＳＴＰｌｕｓ（ＪＤｒｅａｍＩＩＩ）

ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

複数の樹木を含むシミュレーション対象空間内での放射熱輸送現象をシミュレートするシミュレーション装置において、
前記シミュレーション対象空間を、複数の面積要素と複数の体積要素とで定義した仮想的な３次元空間であって、前記シミュレーション対象空間内の前記複数の樹木の樹冠が、それぞれ、１つ以上の透過性を有する体積要素で定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、
２つの面積要素ｉ及びｊに関する形態係数として、面積要素ｉから面積要素ｊを見た形態係数と、面積要素ｊから面積要素ｉを見た形態係数とを算出し、
面積要素ｉと体積要素ｋとに関する形態係数として、面積要素ｉから体積要素ｋを見た形態係数と、体積要素ｋから面積要素ｉを見た形態係数とを算出し、
２つの体積要素ｋ及びｌに関する形態係数として、体積要素ｋから体積要素ｌを見た形態係数と、体積要素ｌから体積要素ｋを見た形態係数とを算出する
形態係数算出手段と、
前記形態係数算出手段により算出された各形態係数を用いて前記複数の面積要素及び前記複数の体積要素の中の各要素間で授受される放射熱量を算出する放射熱量算出手段と、
を備え、
前記形態係数算出手段は、
面積要素ｉから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出手段により算出される、面積要素ｉから体積要素ｋに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから面積要素ｉを見た前記形態係数として、それを用いて前記放射熱量算出手段により算出される、体積要素ｋから面積要素ｉに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから体積要素ｌを見た前記形態係数として、それを用いて前記放射熱量算出手段により算出される、体積要素ｋから体積要素ｌに到達する放射熱量が、体積要素ｋを透過する放射熱量と体積要素ｌを透過する放射熱量とに応じた分だけ減少することになる形態係数を算出し、
体積要素ｌから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出手段により算出される、体積要素ｌから体積要素ｋに到達する放射熱量が、体積要素ｌを透過する放射熱量と体積要素ｋを透過する放射熱量とに応じた分だけ減少することになる
形態係数を算出する、
ことを特徴とするシミュレーション装置。

【請求項2】

前記形態係数算出手段は、面積要素ｉから面積要素ｊを見た前記形態係数、面積要素ｊから面積要素ｉを見た前記形態係数、面積要素ｉから体積要素ｋを見た前記形態係数、体積要素ｋから面積要素ｉを見た前記形態係数、体積要素ｋから体積要素ｌを見た前記形態係数、及び、体積要素ｌから体積要素ｋを見た前記形態係数のそれぞれとして、熱量の放出源側の要素から熱量の到達先側の要素までの間における放射熱量の透過率を考慮に入れた形態係数を算出する、
ことを特徴とする請求項１に記載のシミュレーション装置。

【請求項3】

前記放射熱量算出手段により算出された放射熱量に基づき、前記複数の面積要素及び前記複数の体積要素の中の各要素の温度を算出する温度算出手段を、さらに備える
ことを特徴とする請求項１又は２に記載のシミュレーション装置。

【請求項4】

前記温度算出手段は、各要素の温度を陰解法により算出する
ことを特徴とする請求項３に記載のシミュレーション装置。

【請求項5】

複数の樹木を含むシミュレーション対象空間内での放射熱輸送現象をシミュレートするシミュレーション方法において、
コンピュータが、
前記シミュレーション対象空間を、複数の面積要素と複数の体積要素とで定義した仮想的な３次元空間であって、前記シミュレーション対象空間内の前記複数の樹木の樹冠が、それぞれ、１つ以上の透過性を有する体積要素で定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、
２つの面積要素ｉ及びｊに関する形態係数として、面積要素ｉから面積要素ｊを見た形態係数と、面積要素ｊから面積要素ｉを見た形態係数とを算出し、
面積要素ｉと体積要素ｋとに関する形態係数として、面積要素ｉから体積要素ｋを見た形態係数と、体積要素ｋから面積要素ｉを見た形態係数とを算出し、
２つの体積要素ｋ及びｌに関する形態係数として、体積要素ｋから体積要素ｌを見た形態係数と、体積要素ｌから体積要素ｋを見た形態係数とを算出する
形態係数算出ステップと、
前記形態係数算出ステップにより算出された各形態係数を用いて前記複数の面積要素及び前記複数の体積要素の中の各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
を行い、
前記形態係数算出ステップが、
面積要素ｉから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、面積要素ｉから体積要素ｋに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから面積要素ｉを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｋから面積要素ｉに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから体積要素ｌを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｋから体積要素ｌに到達する放射熱量が、体積要素ｋを透過する放射熱量と体積要素ｌを透過する放射熱量とに応じた分だけ減少することになる形態係数を算出し、
体積要素ｌから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｌから体積要素ｋに到達する放射熱量が、体積要素ｌを透過する放射熱量と体積要素ｋを透過する放射熱量とに応じた分だけ減少することに
なる形態係数を算出する、ステップである、
ことを特徴とするシミュレーション方法。

【請求項6】

複数の樹木を含むシミュレーション対象空間内での放射熱輸送現象をシミュレートするシミュレーションプログラムにおいて、
コンピュータに、
前記シミュレーション対象空間を、複数の面積要素と複数の体積要素とで定義した仮想的な３次元空間であって、前記シミュレーション対象空間内の前記複数の樹木の樹冠が、それぞれ、１つ以上の透過性を有する体積要素で定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、
２つの面積要素ｉ及びｊに関する形態係数として、面積要素ｉから面積要素ｊを見た形態係数と、面積要素ｊから面積要素ｉを見た形態係数とを算出し、
面積要素ｉと体積要素ｋとに関する形態係数として、面積要素ｉから体積要素ｋを見た形態係数と、体積要素ｋから面積要素ｉを見た形態係数とを算出し、
２つの体積要素ｋ及びｌに関する形態係数として、体積要素ｋから体積要素ｌを見た形態係数と、体積要素ｌから体積要素ｋを見た形態係数とを算出する
形態係数算出ステップと、
前記形態係数算出ステップにより算出された各形態係数を用いて前記複数の面積要素及び前記複数の体積要素の中の各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
を実行させ、
前記形態係数算出ステップが、
面積要素ｉから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、面積要素ｉから体積要素ｋに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから面積要素ｉを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｋから面積要素ｉに到達する放射熱量が、体積要素ｋを透過する放射熱量に応じた分だけ減少することになる形態係数を算出し、
体積要素ｋから体積要素ｌを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｋから体積要素ｌに到達する放射熱量が、体積要素ｋを透過する放射熱量と体積要素ｌを透過する放射熱量とに応じた分だけ減少することになる形態係数を算出し、
体積要素ｌから体積要素ｋを見た前記形態係数として、それを用いて前記放射熱量算出ステップにより算出される、体積要素ｌから体積要素ｋに到達する放射熱量が、体積要素ｌを透過する放射熱量と体積要素ｋを透過する放射熱量とに応じた分だけ減少することになる形態係数を算出する、ステップである、
ことを特徴とするシミュレーションプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、放射熱輸送現象をシミュレートするためのシミュレーション装置とシミュレーション方法とシミュレーションプログラムとに関する。

【背景技術】

【0002】

都市部の気温がその周辺部の気温より高くなる現象を、ヒートアイランド現象という。
ヒートアイランド現象を緩和するための対策として、街路樹の配置や緑地の整備が着目されている。そのため、街路樹の配置等を考慮して都市部の放射熱輸送をシミュレートする様々なシミュレーション技術（特許文献１及び非特許文献１〜９参照）が開発されている。

【0003】

ただし、既存のシミュレーション技術は、計算コスト（計算量）を抑えるために、『樹冠による放射熱の散乱を無視し、樹冠による放射熱の減衰・吸収のみを考慮する』といったように、樹冠に関する放射エネルギーの伝達過程を単純化したものとなっている。そのため、既存のシミュレーション技術では、樹冠付近や樹冠内部での放射エネルギーの吸収量等を適切にシミュレートできないことが多いのが現状である。

【先行技術文献】

【特許文献】

【0004】

【特許文献1】特開２００３−０９９６９７号公報

【特許文献2】特許５１３７０３９号公報

【特許文献3】特開２０１２−０２１６８４号公報

【非特許文献】

【0005】

【非特許文献1】吉田 伸治, 村上 周三, 持田 灯, 大岡 龍三, 富永 禎秀, 「環境緩和効果を総合的に組み込んだ新しい3次元樹木モデルの開発」, 生産研究, 51巻1号, 1999年1月

【非特許文献2】吉田 伸治, 大岡 龍三, 持田 灯, 富永 禎秀, 村上 周三, 「樹木モデルを組み込んだ対流・放射・湿気輸送連成解析による樹木の屋外温熱環境緩和効果の検討」, 日本建築学会計画系論文集, No.536, pp.87-94, 2000

【非特許文献3】坂本 雄三, 小島 悦史, 足永 靖信, 今野 雅, 「CFDを利用した樹木のクールスポット効果の数値解析 その１ ―樹木における放射と蒸散に関する計算モデル」, 日本建築学会大会学術講演梗概集, D-1, pp.689-690, 2005

【非特許文献4】小島 悦史, 坂本 雄三,足永 靖信, 今野 雅, 「CFDを利用した樹木のクールスポット効果の数値解析 その２ ―クールスポット効果のケーススタディ」, 日本建築学会大会学術講演梗概集, D-1, pp.691-692, 2005

【非特許文献5】大黒 雅之, 森川 泰成, 「街区及び敷地レベルを対象としたヒートアイランド解析評価システムの開発」, 大成建設技術センター報, No.38, pp.14-1-14-8, 2005

【非特許文献6】大黒 雅之, 森川 泰成, 本橋 比奈子, 「街区スケールを対象としたヒートアイランド解析評価プログラムの開発と適用」, 大成建設技術センター報, No.42, pp.49-1-49-8, 2009

【非特許文献7】片岡 浩人, 大塚 清敏, 赤川 宏幸, 「屋外熱環境評価のための数値予測モデルの開発 −樹木による冷却効果のモデル化−」, 日本建築学会学術講演梗概集, D-1, pp.927-928, 2008

【非特許文献8】佐々木 澄, 「数値解析に基づく街路樹がストリートキャニオン内の熱空気環境に及ぼす影響の検討」, 清水建設研究報告, No.85, pp.41-50, 2007

【非特許文献9】H.B. リジャル, 大岡 龍三 他, 「数値解析による大規模緑地のヒートアイランド緩和効果の検討」, 東京大学生産技術研究所, 生産研究, 62巻1号, 2010

【発明の概要】

【発明が解決しようとする課題】

【0006】

本発明の課題は、樹冠を含む３次元空間における放射熱輸送現象を低計算コストで良好にシミュレートできる技術を提供することにある。

【課題を解決するための手段】

【0007】

上記課題を解決するために、本発明の、放射熱輸送現象をシミュレートするシミュレーション装置は、
複数の面積要素と複数の体積要素とで定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出手段と、
前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出手段と、
を備え、
前記３次元空間内に存在する複数の樹木の樹冠が前記複数の体積要素として取り扱われるように、前記３次元空間が前記複数の面積要素と前記複数の体積要素とにより定義されている
装置として構成される。

【0008】

すなわち、本発明のシミュレーション装置は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置は、２つの体積要素に関する形態係数として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、本発明のシミュレーション装置によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。

【0009】

本発明は、上記シミュレーション装置と同様の特徴を有するシミュレーション方法や、情報処理装置（コンピュータ）を上記シミュレーション装置として機能させるシミュレーションプログラムとしても実現することが出来る。尚、本発明を、上記シミュレーションプログラムを記録したコンピュータ可読媒体としても実現することも出来る。

【発明の効果】

【0010】

本発明によれば、樹冠を含む３次元空間における放射熱輸送現象を低計算コストで良好にシミュレートすることができる。

【図面の簡単な説明】

【0011】

【図1】本発明の一実施形態に係るシミュレーションシステムの構成図

【図2】実施形態に係るシミュレーション装置の機能の説明図

【図3】２面積要素に関する形態係数についての算出式中のパラメータの説明図

【図4】面積要素から体積要素をみた形態係数についての算出式中のパラメータの説明図

【図5】体積要素の熱収支の説明図

【発明を実施するための形態】

【0012】

以下、図面を参照して、本発明の一実施形態について説明する。尚、以下の実施形態の説明は、例示であり、本発明は、実施形態の構成に限定されない。

【0013】

図１に、実施形態に係るシミュレーションシステムの構成を示す。図示してあるように、本実施形態に係るシミュレーションシステムは、シミュレーション装置１０、入力装置２０及び出力装置３０を含む。また、シミュレーション装置１０は、演算部１２、記憶部１４、インタフェース回路（Ｉ／Ｆ）１６を備える。

【0014】

シミュレーション装置１０のインタフェース回路（Ｉ／Ｆ）１６は、演算部１２が他装置と通信を行うために使用する回路である。記憶部１４は、シミュレーションプログラム１８を記憶した不揮発性の記憶装置である。この記憶部１４は、演算部１２が利用するデータや演算部１２による処理結果の記憶にも使用される。

【0015】

演算部１２は、ＣＰＵ(Central Processing Unit)、ＲＡＭ(Random Access Memory)等
を組み合わせたユニットである。演算部１２は、シミュレーションプログラム１８を記憶部１４から読み出して実行することにより、各種処理（詳細は後述）を行う。また、演算部１２は、シミュレーションプログラム１８を実行することにより、本発明に係る形態係数算出手段、放射熱量算出手段及び温度算出手段として機能する。

【0016】

入力装置２０は、シミュレーション装置１０に情報を入力するための装置である。入力装置２０は、キーボード、マウス等のポインティングデバイス等を含む。また、出力装置３０は、シミュレーション装置１０からの情報を出力するための、ＬＣＤ(Liquid Crystal Display)やＣＲＴ(Cathode-Ray Tube)等のディスプレイ、プリンタ等である。

【0017】

尚、シミュレーション装置１０は、通常、行列演算が高速に行えるコンピュータ（ベクトル型並列計算機等）にシミュレーションプログラム１８を実行させることによって実現される。従って、通常、シミュレーション装置１０（ベクトル型並列計算機等）に接続された入出力用のコンピュータの入力装置及び出力装置が、入力装置２０及び出力装置３０として機能する。

【0018】

以下、シミュレーション装置１０の機能を説明する。
シミュレーション装置１０は、複数の樹木を含む都市空間（以下、シミュレーション対象空間と表記する）内での放射熱輸送現象をシミュレートするための装置である。

【0019】

図２に模式的に示してあるように、シミュレーション装置１０の使用時には、シミュレーション対象時間範囲内の太陽の高度や方位に関する情報等が設定された計算条件設定ファイルが記憶部１４に記憶される。また、シミュレーション装置１０に対して、地形データ、建物形状データ、樹木分布データ、地表面温度データ、建物表面温度データ及び葉表面温度データが入力される。

【0020】

地形データは、シミュレーション対象空間の地形（地面の形状）を示すデータである。建物形状データは、シミュレーション対象空間内に存在する各建物の位置及び形状を示すデータである。樹木分布データは、シミュレーション対象空間内に存在する各樹木の位置及び形状と各樹木の葉面積密度分布とを示すデータである。

【0021】

シミュレーション装置１０に入力されるこれらのデータは、シミュレーション対象空間の構造（３次元都市形状及び３次元樹木分布）が分かるデータであれば良い。従って、地形データや建物形状データとして、例えば、地面／建物の高さの２次元平面データ（地面
／建物の高さと座標との対応関係を表すデータ）を使用することができる。また、樹木分布データとして、例えば、樹木インデックス（樹木の識別情報）の２次元平面データと、各樹木インデックスで識別される樹木の葉面積密度の鉛直分布を示すデータとを含むデータを使用することができる。

【0022】

シミュレーション装置１０に入力される地表面温度データ、建物表面温度データ、葉表面温度データは、それぞれ、地表面の各部の温度の初期値を示すデータ、各建物表面の各部の温度の初期値を示すデータ、各所の葉表面の温度の初期値を示すデータである。

【0023】

以下、シミュレーション装置１０（演算部１２）が行う処理の内容を説明する。尚、以下の説明では、シミュレーション装置１０に入力される地形データ、建物形状データ及び樹木分布データのことを、構造指定データと表記する。また、シミュレーション装置１０に入力される地表面温度データ、建物表面温度データ及び葉表面温度データのことを、初期値データと表記する。

【0024】

シミュレーション装置１０は、基本的には、入力された構造指定データが示しているシミュレーション対象空間内の各部の温度を、計算条件設定ファイル内の各種情報及び入力された初期値データを用いて、Δｔ毎にシミュレートする装置である。

【0025】

図２に示してあるように、シミュレーション装置１０の演算部１２が行う処理は、前処理と主処理とに分類できる。

【0026】

前処理は、シミュレーション用データ生成処理（ステップＳ１０１）とパラメータ算出処理（ステップＳ１０２）とがこの順に行われる処理である。

【0027】

ステップＳ１０１にて行われるシミュレーション用データ生成処理は、入力された構造指定データに基づき、『シミュレーション対象空間内の各建物の壁面や地面を複数の面積要素として取り扱い、シミュレーション対象空間内の各樹木の樹冠を１つ以上の透過性を有する体積要素として取り扱うためのシミュレーション用データ』を生成する処理である。

【0028】

このシミュレーション用データ生成処理により生成されるシミュレーション用データは、後述する形態係数の算出に必要な情報（各面積／体積要素の形状やシミュレーション対象空間内での位置等）が分かるものであれば良い。従って、シミュレーション用データを、例えば、面積／体積要素毎に、『通し番号、対応するシミュレーション対象空間内の計算格子の座標番号、面積要素の向いている方向／体積要素であることを示すフラグ、樹冠であるか否かを示すフラグ等からなるデータ』を含むデータとしておくことができる。

【0029】

ステップＳ１０２にて行われるパラメータ算出処理は、シミュレーション用データ生成処理により生成されたシミュレーション用データに基づき、主処理で使用される各種パラメータを算出する処理である。

【0030】

このパラメータ算出処理時に算出されるパラメータには、各２要素（面積／体積要素）に関する形態係数Ｆ、各体積要素ｋの有効表面積＜A^eff＞_k、各要素ｉの天空率ω_i、各面積要素ｉの日陰率Ｄ_i、各体積要素ｋの有効日陰率Ｄ^eff_kがある。

【0031】

まず、パラメータ算出処理時に算出される各２要素に関する形態係数Ｆについて説明する。

【0032】

パラメータ算出処理時には、２つの面積要素ｉ，ｊの組み合わせ毎に、『面積要素ｉか
ら面積要素ｊを見た形態係数Ｆ_ij』と『面積要素ｊから面積要素ｉを見た形態係数Ｆ_ji』とが、算出される。また、面積要素ｉと体積要素ｋの組み合わせ毎に、『面積要素ｉから体積要素ｋを見た形態係数Ｆ_ik』と『体積要素ｋから面積要素ｉを見た形態係数Ｆ_ki』とが、算出される。さらに、２つの体積要素ｋ及びｌの組み合わせ毎に、『体積要素ｋから体積要素ｌを見た形態係数Ｆ_kl』及び『体積要素ｌから体積要素ｋを見た形態係数Ｆ_lk』が算出される。

【0033】

パラメータ算出処理時に算出される『面積要素ｉから面積要素ｊを見た形態係数Ｆ_ij』、『面積要素ｊから面積要素ｉを見た形態係数Ｆ_ji』は、それぞれ、以下の（１）、（２）式で定義される値である。

【0034】

【数1】

【0035】

これらの式において、Ａ_i、Ａ_jは、それぞれ、面積要素ｉの面積、面積要素ｊの面積である。β_i、β_j は、図３に模式的に示してあるように、それぞれ、面積要素ｉ内の微小
面ｄＡ_iと面積要素ｊ内の微小面ｄＡ_jとを結ぶ直線と微小面ｄＡ_iの法線のなす角度、当
該直線と微小面ｄＡ_jの法線のなす角度である。また、ｒは、微小面ｄＡ_i・微小面ｄＡ_j
間の距離である。

【0036】

Ｔ_ijは、面積要素ｉと面積要素ｊの間の透過率である。Ｔ_ijは、２微小面ｄＡ_i及びｄ
Ａ_j 間の光学的厚さτ_ijを用いて次式により算出される。

【0037】

【数2】

【0038】

また、面積要素ｉと面積要素ｊとの間（微小面ｄＡ_iと微小面ｄＡ_jとの間）に樹冠が分布している場合、光学的厚さτ_ijは、樹冠の消散係数ｋ_ext及び葉面積密度ａを用いて次
式により算出される。

【0039】

【数3】

【0040】

尚、形態係数の具体的な算出手順については後述するが、上記した（１）、（２）式で定義される形態係数は、相反則を満たしている。すなわち、面積Ａ_iと面積Ａ_jと形態係数Ｆ_ijと形態係数Ｆ_jiとの間には、以下の関係がある。

【0041】

【数4】

【0042】

従って、Ｆ_jiを、（１）式から算出したＦ_ijとＡ_iとＡ_jとから算出することも、Ｆ_ijを、（２）式から算出したＦ_jiとＡ_iとＡ_jとから算出することもできる。

【0043】

パラメータ算出処理時に算出される『面積要素ｉから体積要素ｋを見た形態係数Ｆ_ik』は、以下の（３）式により定義される値である。

【0044】

【数5】

【0045】

この（３）式におけるβ_iは、図４に模式的に示してあるように、面積要素ｉの微小面
ｄＡ_iと体積要素ｊ内の微小投影面ｄＡ_k←_iとを結ぶ直線と微小面ｄＡ_iの法線のなす角度である。また、ｒは、微小面ｄＡ_i・微小投影面ｄＡ_k←_i間の距離である。

【0046】

Ａ^eff_k←_iは、体積要素ｋ自身の遮蔽率を考慮した，面積要素ｉの方向から見た体積要
素ｋの有効面積である。このＡ^eff_k←_iは、次式により算出される。

【0047】

【数6】

【0048】

ここで、Δτ_k←_iは、微小投影面ｄＡ_k←_i(図４参照)に垂直な方向の体積要素ｋの光学的厚さである。体積要素ｋが樹木の場合、Δτ_k←_iは、消散係数ｋ_ext、葉面積密度ａお
よびｄＡ_k←_iに垂直な方向の体積要素ｋの幾何学的厚さΔｓ_k←_iを用いて、次式により算出される。

【0049】

【数7】

【0050】

要するに、上記した（３）式は、（４）式を用いれば、以下のように表記できる式となっている。パラメータ処理時には、以下の（５）式に従って、面積要素ｉから体積要素ｋを見た形態係数Ｆ_ikが算出される。

【0051】

【数8】

【0052】

パラメータ算出処理時に算出される『体積要素ｋから面積要素ｉを見た形態係数Ｆ_ki』
は、次式により定義される値である。

【0053】

【数9】

【0054】

換言すれば、体積要素ｋから面積要素ｉを見た形態係数Ｆ_kiとしては、以下の式で表される相反則を満たす値が算出される。

【0055】

【数10】

【0056】

パラメータ算出処理時に算出される『体積要素ｋから体積要素ｌを見た形態係数Ｆ_kl』は、次式により定義される値である。

【0057】

【数11】

【0058】

すなわち、この形態係数Ｆ_klは、体積要素ｋの遮蔽率（“1-exp(-Δτ_k←_l)”）と体積要素ｌの遮蔽率（“1-exp(-Δτ_l←_k)”）と体積要素ｋ、ｌ間の透過率Ｔ_klとを考慮に入れた形態係数となっている。

【0059】

尚、（７）式が表している体積要素に関する形態係数は、面積要素の形態係数と同様に次式の相反則を満たす。

【0060】

【数12】

【0061】

従って、他の形態係数と同様に、２体積要素間の形態係数の双方を、（７）式に従って算出することも、一方の形態係数を（７）式に従って算出し、他方の形態係数を、当該一方の形態係数の算出結果から算出することも出来る。

【0062】

パラメータ算出処理時には、上記した各形態係数が、モンテカルロ法により算出される。

【0063】

すなわち、形態係数の算出時には、０から１までの範囲の一様乱数Ｒθ、Ｒφから、以下の式によりμ、φを算出し、算出結果に基づき、以下に示した内容の単位ベクトルｎを生成する処理が、必要とされる形態係数の精度に応じた回数、繰り返される。

【0064】

【数13】

【0065】

換言すれば、形態係数の算出時には、Lambertの余弦則にしたがって面積要素（又は体
積要素）から射出される多数の光子の進行方向単位ベクトルｎを生成する処理が行われる。

【0066】

そして、生成した各単位ベクトルｎの方向に要素ｉから等しいエネルギーＷ₀を持った
光子が射出された場合に要素ｊに入射されることになる各光子ｐのエネルギーＷ_pを積算
することにより、形態係数Ｆ_ijが求められる。

【0067】

具体的には、面積要素ｉから面積要素ｊをみた形態係数Ｆ_ijは、面積要素ｉから射出させた光子の個数がＮ、面積要素ｊに入射された光子の個数がｎであった場合には、以下の式により算出される。

【0068】

【数14】

【0069】

なお、建物壁面などの完全遮蔽性の面積要素のみの場合には、面積要素ｊに入射する光子のエネルギーＷ_pはＷ₀と等しいため、形態係数は、Ｆ_ij＝ｎ／Ｎにより算出できる。樹木のような放射透過性の要素が空間中に分布する場合には、面積要素ｊに入射する光子のエネルギーＷ_pは、面積要素ｉから面積要素ｊに到達するまでの間に減衰する。そのため
、その減衰の形態係数に対する影響を光子の持つエネルギーを減衰させることで考慮する。すなわち、次式により、Ｗ_pが計算される。

【0070】

【数15】

【0071】

ここで、Ｔ_ij,pは、光子ｐの経路に沿った透過率である。なお、光子が透過性の体積要素内で確率的に遮蔽されると考えて、面積要素ｊに入射する光子の数を減少させても、減衰の影響を考慮した形態係数を求めることが出来る。

【0072】

また、面積要素ｉから体積要素ｋをみた形態係数Ｆ_ikは、（７）式により定義される値である。従って、面積要素ｉから体積要素ｋをみた形態係数Ｆ_ikは、次式により算出される。

【0073】

【数16】

【0074】

ここで、Δτ_k←_i,pは、光子ｐの経路に沿った、体積要素ｋの内部での光学的
厚さである。Ｎは、面積要素ｉから射出された光子数であり、ｎは、体積要素ｋに入射された光子数である。

【0075】

体積要素ｋから面積要素ｉをみた形態係数Ｆ_kiの算出時には、体積要素の表面上の各点からの光子の等方的な射出が仮定される。すなわち、体積要素の表面上の各点からLambertの余弦則に従って多数の光子が仮想的に射出される。

【0076】

そして、体積要素ｋから面積要素ｉをみた形態係数Ｆ_kiが、次式により算出される。

【0077】

【数17】

【0078】

ここで、Δτ_k←_i,pは、光子ｐの射出点から光子ｐの進行方向とは逆の方向へ
バックトレースすることによって計算される，体積要素ｋの内部での光学的厚さである。また、Ｓ_kは、体積要素ｋの表面積であり、Ｍは、体積要素ｋから射出された光子数であ
り、ｍは、面積要素ｉに入射された光子数である。

【0079】

同様に、体積要素ｋから体積要素ｌをみた形態係数Ｆ_klは、次式により計算される。

【0080】

【数18】

【0081】

尚、この式におけるｍは、Ｍ個の光子を体積要素ｋから射出した結果として体積要素ｌに入射された光子数である。

【0082】

以下、パラメータ算出処理で算出される、各体積要素ｋの有効表面積＜A^eff＞_k、各要
素ｉの天空率ω_i、各面積要素ｉの日陰率Ｄ_i、各体積要素ｋの有効日陰率Ｄ^eff_kについて説明する。

【0083】

パラメータ算出処理時に算出される各体積要素ｋの有効表面積＜A^eff＞_kは、以下の式
により定義される値である。

【0084】

【数19】

この式において、ｍは、体積要素ｋの周囲に存在する（体積要素ｋから見える）要素（面積要素又は体積要素）の総数であり、ｉ（＝１〜ｍ）は、体積要素ｋの周囲に存在する面積要素又は体積要素の要素番号である。

【0085】

パラメータ算出処理時には、この有効表面積＜Ａ^eff＞_kが、次式により算出される。

【0086】

【数20】

【0087】

すなわち、有効表面積＜Ａ^eff＞_kもモンテカルロ法により算出される。

【0088】

要素（面積要素／体積要素）ｉの天空率ω_iは、要素ｉから空をみた形態係数に相当す
る値である。天空率ω_iは、要素ｉから面積要素をみた形態係数と同様の手順で算出され
る。

【0089】

面積要素ｉの日陰率Ｄ_iは、面積要素iから太陽側へ射出された光子ｐが、他の要素に入射した際に失うエネルギーΔＷｐを積算することにより求められる。より具体的には、日陰率Ｄ_iは、モンテカルロ法を用いて以下の式により算出される。

【0090】

【数21】

【0091】

尚、上記式において、Ｎは、面積要素iから射出された光子数である。

【0092】

同様に、体積要素ｋの有効日陰率Ｄ^eff_kは、モンテカルロ法を用いて以下の式により算出される。

【0093】

【数22】

【0094】

上記した各面積要素ｉの日陰率Ｄ_i、各体積要素ｋの有効日陰率Ｄ^eff_kは、太陽の位置
により値が変わるパラメータである。従って、シミュレーションを行う時間範囲が比較的に長い場合には、パラメータ算出処理時に、当該時間範囲内の各時刻における日陰率Ｄ_i
及び有効日陰率Ｄ^eff_kが算出される。

【0095】

以下、演算部１２が行う主処理（図２）の内容を説明する。
演算部１２が行う主処理は、放射熱フラックス等算出処理（ステップＳ２０１）と温度算出処理（ステップＳ２０２）とが、全時間ステップ数N_tの回数、繰り返される処理である。尚、全時間ステップ数N_tは、シミュレーション対象時刻範囲及び時間刻み幅Δtに基
づいて決定することができる。全時間ステップ数N_tの指定は、計算条件設定ファイル内に全時間ステップ数N_t又はシミュレーション対象時刻範囲及び時間刻み幅Δtの情報を設定
しておくことや、入力装置２０を用いて入力することによって、行われる。

【0096】

ステップＳ２０１にて行われる放射熱フラックス等算出処理は、パラメータ算出処理で算出されたパラメータ（形態係数等）を用いて、各要素ｉから射出される短波放射（可視光線）の放射フラックスＧ_S,i[W/m²]及び長波放射（赤外線）の放射フラックスＧ_L,i [W/m²] を算出し、算出した放射フラックスを用いて、各要素ｉに吸収される、短波放射に関する正味の放射熱Ｒ_S,i[W]及び長波放射に関する正味の放射熱Ｒ_L,i[W]を算出する処理である。

【0097】

具体的には、各要素ｉから射出される放射フラックスＧ_S,i、Ｇ_L,i [W/m²] については、それぞれ、以下の式が成立する。

【0098】

【数23】

【0099】

尚、（８）、（９）式におけるｎは、面積要素及び体積要素の総数である。<Ａ^eff>ｉ
は、要素ｉが体積要素である場合には、体積要素ｉの有効表面積であり、要素ｉが面積要素である場合には、面積要素ｉの面積である。

【0100】

α_S,i、α_L,iは、それぞれ、要素iの、短波放射、長波放射に関する反射率であり、ε_iは、要素ｉの射出率である。Ｓ_direct,iは、要素ｉに入射する太陽からの直達短波放射フラックスであり、Ｓ_diffuse,iは、要素iに入射する大気散乱短波放射の放射フラックスである。Ｌ_iは、要素ｉに入射する大気長波放射の放射フラックスであり、Ａ^eff_i←_Solar、Ａ^eff_i←_skyは、それぞれ、要素iの、太陽の方向、空の方向の有効面積である。

【0101】

Ｂ(Ｔ_i)は、要素ｉから熱放射により射出される放射フラックスである。長波放射の放
射フラックスとして、上記したＧ_L,jのみを算出する場合（長波放射に関する放射フラッ
クスを波長範囲別に算出しない場合）、Ｂ(Ｔ_i)は、ステファン・ボルツマン定数σを用
いて次式により算出される。

【0102】

【数24】

【0103】

Ｓ_direct,i、Ｓ_diffuse,iは、パラメータ算出処理で算出された天空率ω_i、日陰率Ｄ_i
を用いて以下の式により算出される。

【0104】

【数25】

【0105】

ここで，Ｓ↓は水平面に下向きに入射する太陽放射フラックスであり、Ｓ（＝(S_x,S_y,S_z)）は太陽方位ベクトルである。ｎ_iは、面積要素ｉの単位法線ベクトルであり、c_direct、c_diffuseは、直散分離のための係数である。

【0106】

尚、上記式で算出されるＳ_direct,iは、面積要素ｉのＳ_direct,iである。体積要素ｉのＳ_direct,iは，パラメータ算出処理で算出された有効日陰率Ｄ^eff_iを用いて次式により算出される。

【0107】

【数26】

【0108】

また、射出率ε_iは、要素ｉの吸収率と等しいため、射出率ε_iとして、“１−α_L,_ｉ”を使用することが出来る。

【0109】

【数27】

【0110】

上記した（８）式及び（９）式は、ｉ＝１〜ｎのそれぞれについて成立する。すなわち、以下の２つの線形行列方程式が成立する。

【0111】

【数28】

【0112】

放射フラックス等算出処理時には、これらの線形行列方程式を解くことで、各要素iか
ら射出される放射フラックスＧ_S,i，Ｇ_L,iが、算出される。

【0113】

そして、各要素ｉに吸収される、短波放射に関する正味の放射熱Ｒ_S,i[W]、及び長波放射に関する正味の放射熱Ｒ_L,i[W]が、それぞれ、以下の（１０）、（１１）式により算出される。

【0114】

【数29】

【0115】

温度算出処理（図２；ステップＳ２０２）は、放射熱フラックス等算出処理にて算出された放射熱からシミュレーション対象空間内の各部の温度を算出する処理である。この温度算出処理時における各面積要素の温度の算出手順は、放射熱が、樹冠を透過性を有する体積要素として取り扱った形態係数を用いて算出されたものであることを除けば、一般的な算出手順と同様のものである。そのため、以下では、樹冠（体積要素）の温度の算出手順のみを説明することにする。

【0116】

図５に模式的に示してあるように、樹木の樹冠の一部である体積要素には、短波放射の放射熱フラックスＲ_Sと長波放射の放射熱フラックスＲ_Lとが流入し、当該体積要素からは、顕熱フラックスＨと潜熱フラックスＬＥとが流出する。
従って、樹冠である体積要素ｉに関する熱収支は、次式により表されることになる。

【0117】

【数30】

【0118】

ここで，Ｔ_leaf,iは、要素ｉ内の葉の表面温度[K]であり、a_iは、要素ｉ内の葉の面積
密度[m²/m³]である。Ｖ_iは、要素ｉの体積[m³]であり、Ｃは、単位葉面積あたりの葉の熱容量[J/K/m²]である。Ｒ_S,i、Ｒ_L,iは、それぞれ、葉に吸収される短波放射の正味の放射熱（放射熱フラックスの強度）[W]、長波放射の正味の放射熱[W]であり、Ｌは、蒸発潜熱[J/kg]である。

【0119】

Ｈ_iは、葉から大気に放出される顕熱輸送量（顕熱フラックスの強度）[W]であり、Ｅ_i
は葉から大気に蒸散される水蒸気量[kg/s]である。

【0120】

葉から大気に放出される顕熱輸送量Ｈ_iおよび葉から大気に蒸散される水蒸気量Ｅ_iは、以下の式により算出される。

【0121】

【数31】

【0122】

ここで、Ｔ_air,iは、要素ｉ内の大気温度[K]，f_air,iは、体積要素ｉ内の大気中の水蒸気分圧[Pa]、ｆ_leaf,iは、体積要素ｉ内の葉表面の飽和水蒸気分圧[Pa]，Ｋ_hは、対流熱
伝達係数[W/m²/K]，Ｋ_wは、対流水蒸気輸送係数[kg/s/m²/Pa]、βは蒸発効率である．

【0123】

温度算出処理時には、時間ステップｎにおける葉表面温度及び熱フラックスを用いて、時間刻みΔt後の時間ステップｎ＋１における葉表面温度Ｔ_leaf,iが、計算され
る。

【0124】

具体的には、時間ステップｎからｎ＋１までの間の葉表面温度の変化量ΔＴ_leaf,i
は，時間Δtの経過に起因する葉表面温度の変化による正味長波放射，顕熱輸送量，
蒸散量の変化を考慮すると、次式により表される。

【0125】

【数32】

【0126】

従って，葉表面温度の変化量ΔＴ_leaf,iを、次式により求めることができる。

【0127】

【数33】

【0128】

温度算出処理時には、この式により、葉表面温度の変化量ΔＴ_leaf,iが求められて
から、時間ステップｎ+１における葉表面温度Ｔ_leaf,i⁽ⁿ⁺¹⁾が次式により、算出される。

【0129】

【数34】

【0130】

温度算出処理は、各部の温度の算出と、算出した各部の温度の出力（本実施形態では、記憶部１４への記憶）とが完了したときに終了する。そして、指定された回数の処理が完了していなかった場合には、再び、放射熱フラックス等算出処理が開始され、指定された回数の処理が完了したときに、主処理が終了する。

【0131】

以上、説明したように、本実施形態に係るシミュレーション装置１０は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数（（５）、（６）式参照）として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置１０は、２つの体積要素に関する形態係数（（７）式参照）として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置１０は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、シミュレーション装置１０によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。

【0132】

《変形形態》
上記した実施形態に係るシミュレーション装置１０は、各種の変形を行えるものである。
例えば、シミュレーション装置１０を、２要素間の透過率Ｔが考慮されていない形態係数を算出し、放射フラックスの算出時等に２要素間の透過率を考慮する装置に変形することが出来る。ただし、形態係数の算出時に２要素間の透過率Ｔを考慮しておいた方が正確な結果が得られるし、計算コストも低くなる。従って、上記した形態係数を採用しておくことが好ましい。

【0133】

また、樹幹の体積熱容量Caは、通常、かなり小さいので、シミュレーション装置１０を、葉表面温度を、Ca=0として算出する装置に変形することも出来る。また、シミュレーション装置１０を、一部又は全ての形態係数として、その定義式の解析解を用いる装置に変形することも出来る。

【0134】

シミュレーション装置１０は、オイラー陰解法により葉温度を算出する装置であるが、シミュレーション装置１０を、陽解法により葉温度を算出する装置や、２次精度のクランク・ニコルソン法により葉温度を算出する装置に変形することも出来る。ただし、陰解法の方が正確な値を得やすいので、葉温度の算出法には、陰解法を採用しておくことが好ましい。

【符号の説明】

【0135】

１０ シミュレーション装置
１２ 演算部
１４ 記憶部
１６ インタフェ−ス（Ｉ／Ｆ）回路

【図1】

【図2】

【図3】

【図4】

【図5】