关于辐射热输送现象的模拟装置、模拟方法和模拟程序技术领域
本发明涉及用于模拟辐射热输送现象的模拟装置、模拟方法和模拟程序。
背景技术
将城市部的气温高于周边部的气温的现象称为热岛现象。作为用于缓解热岛现象
的对策，着眼于路旁树的配置、绿地的配备。因此，研发出考虑路旁树的配置等而模拟城市
部的辐射热输送的各种模拟技术(参照专利文献1和非专利文献1～9)。
但是，现有的模拟技术为了抑制计算成本(计算量)，例如“忽略树冠对于辐射热的
散射，仅考虑树冠对于辐射热的衰减、吸收”，使与树冠相关的辐射能量的传递过程简单化。
因此，现状是凭借已有的模拟技术，大多无法适当地模拟树冠附近或树冠内部的辐射能量
的吸收量等。
专利文献1：日本专利公开公报特开2003-099697号
专利文献2：日本专利公报第5137039号
专利文献3：日本专利公开公报特开2012-021684号
非专利文献1：吉田伸治，村上周三，持田灯，大冈龙三，富永祯秀，“综合地编入环
境缓解效果的新的三维树木模型的研发(環境緩和効果を総合的に組み込んだ新しい3次
元樹木モデルの開発)”，生产研究，51卷1号，1999年1月
非专利文献2：吉田伸治，大冈龙三，持田灯，富永祯秀，村上周三，“设置有树木模
型的对流、辐射、湿气输送耦合解析所产生的树木的室外温热环境缓解效果的研究(樹木モ
デルを組み込んだ対流·放射·湿気輸送連成解析による樹木の屋外温熱環境緩和効果
の検討)”，日本建筑学会计划系论文集，No.536，pp.87-94，2000
非专利文献3：坂本雄三，小岛悦史，足永靖信，今野雅，“利用CFD的树木的冷汇效
果的数值解析其1―关于树木的辐射和蒸腾的计算模型(CFDを利用した樹木のクールスポ
ット効果の数値解析その1―樹木における放射と蒸散に関する計算モデル)”，日本建筑学
会大会学术演讲概要集，D-1，pp.689-690，2005
非专利文献4：小岛悦史，坂本雄三，足永靖信，今野雅，“利用CFD的树木的冷汇效
果的数值解析其2―冷汇效果的专题(CFDを利用した樹木のクールスポット効果の数値解析
その2―クールスポット効果のケーススタディ)”，日本建筑学会大会学术演讲概要集，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.rijal，大冈龙三等，“基于数值解析的大规模绿地的热岛缓解
效果的研究(数値解析による大規模緑地のヒートアイランド緩和効果の検討)”，东京大学
生产技术研究所，生产研究，62卷1号，2010
发明内容
本发明的课题在于提供能够以低计算成本良好地模拟含有树冠的三维空间的辐
射热输送现象的技术。
为了解决所述课题，本发明的模拟装置构成为用于模拟辐射热输送现象，其包括：
形态系数计算单元，计算与由多个面积要素和多个体积要素定义的虚拟三维空间内的各2
要素相关的形态系数，并且计算减少与透过一个或者两个体积要素的辐射热量相应的大小
的值的形态系数，作为与含有该一个或者两个体积要素的2要素相关的形态系数；以及辐射
热量计算单元，使用由所述形态系数计算单元计算出的各形态系数，计算在各要素之间交
换的辐射热量，以将存在于所述三维空间内的多颗树木的树冠视为所述多个体积要素的方
式，由所述多个面积要素和所述多个体积要素定义所述三维空间。
即，本发明的模拟装置将各树冠视为一个以上的具有透过性的体积要素，且计算
如下的形态系数作为与一个体积要素和一个面积要素相关的形态系数，该形态系数减少的
值与透过该一个体积要素的辐射热量对应。另外，模拟装置计算如下的形态系数作为与两
个体积要素相关的形态系数，该形态系数减少的值减少与透过该两个体积要素的辐射热量
对应。此外，模拟装置使用计算出的各形态系数，来计算在各要素之间交换的辐射热量。因
此，根据本发明的模拟装置，能够以无需另行计算树干内部的状态的方式(换句话说，以低
计算成本)，良好地模拟含有树冠的三维空间的辐射热输送现象。
本发明也能够作为具有与所述模拟装置相同的特征的模拟方法、使信息处理装置
(计算机)作为所述模拟装置发挥功能的模拟程序来实现。此外，也可以将本发明作为记录
有所述模拟程序的计算机可读介质来实现。
根据本发明，能够以低计算成本良好地模拟含有树冠的三维空间的辐射热输送现
象。
附图说明
图1为表示本发明的一个实施方式的模拟系统的结构图。
图2为实施方式的模拟装置的功能的说明图。
图3为与2面积要素相关的形态系数的计算式中的参数的说明图。
图4为从面积要素观察体积要素的形态系数的计算式中的参数的说明图。
图5为体积要素的热平衡的说明图。
附图标记说明
10模拟装置
12运算部
14存储部
16接口(I/F)电路
具体实施方式
以下，参照附图对本发明的一个实施方式进行说明。此外，以下的实施方式的说明
为例示，本发明并不局限于实施方式的结构。
图1中示出实施方式的模拟系统的结构。如图所示，本实施方式的模拟系统包括模
拟装置10、输入装置20和输出装置30。另外，模拟装置10具备运算部12、存储部14和接口电
路(I/F)16。
模拟装置10的接口电路(I/F)16为用于运算部12与其他装置进行通信的电路。存
储部14为存储模拟程序18的非易失性的存储装置。该存储部14也被用于存储运算部12所使
用的数据和运算部12所进行的处理的处理结果。
运算部12为由CPU(CentralProcessingUnit)、RAM(RandomAccessMemory)等
组合而成的单元。运算部12通过从存储部14读取并执行模拟程序18，从而进行各种处理(详
情后述)。另外，运算部12通过执行模拟程序18，从而作为本发明的形态系数计算单元、辐射
热量计算单元和温度计算单元发挥功能。
输入装置20为用于向模拟装置10输入信息的装置。输入装置20包括键盘、鼠标等
定点设备等。另外，输出装置30为用于输出来自模拟装置10的信息的LCD(LiquidCrystal
Display)、CRT(Cathode-RayTube)等显示器、打印机等。
此外，通常通过使可高速进行矩阵运算的计算机(向量型并行计算机等)执行模拟
程序18来实现模拟装置10。因此，连接于模拟装置10(向量型并行计算机等)的输入输出用
的计算机的输入装置和输出装置通常作为输入装置20和输出装置30发挥功能。
以下对模拟装置10的功能进行说明。
模拟装置10为用于模拟含有多颗树木的城市空间(以下称为模拟对象空间)内的
辐射热输送现象的装置。
如图2中示意性所示，在使用模拟装置10时，将设定有与模拟对象时间范围内的太
阳的高度、方位相关的信息等的计算条件设定文件存储于存储部14。另外，对模拟装置10输
入地形数据、建筑物形状数据、树木分布数据、地表温度数据、建筑物表面温度数据和叶表
面温度数据。
地形数据为表示模拟对象空间的地形(地面的形状)的数据。建筑物形状数据为表
示存在于模拟对象空间内的各建筑物的位置和形状的数据。树木分布数据为表示存在于模
拟对象空间内的各树木的位置和形状以及各树木的叶面积密度分布的数据。
向模拟装置10输入的所述数据只要是可得知模拟对象空间的结构(三维城市形状
和三维树木分布)的数据即可。因此，作为地形数据、建筑物形状数据，例如可以使用地面/
建筑物的高度的二维平面数据(表示地面/建筑物的高度与坐标的对应关系的数据)。另外，
作为树木分布数据，例如可以使用含有如下数据的数据：树木索引(树木的识别信息)的二
维平面数据；以及表示由各树木索引识别出的树木的叶面积密度的铅直分布的数据。
向模拟装置10输入的地表温度数据、建筑物表面温度数据、叶表面温度数据分别
是表示地表各部分的温度的初始值的数据、表示各建筑物表面的各部分的温度的初始值的
数据、表示各个位置的叶表面的温度的初始值的数据。
以下，对模拟装置10(运算部12)进行的处理的内容进行说明。此外，在以下的说明
中，将输入至模拟装置10的地形数据、建筑物形状数据和树木分布数据称为结构指定数据。
另外，将输入至模拟装置10的地表温度数据、建筑物表面温度数据和叶表面温度数据称为
初始值数据。
模拟装置10基本上为如下装置：使用计算条件设定文件内的各种信息和输入的初
始值数据，对输入的结构指定数据所示的模拟对象空间内的各部分的温度以每个Δt进行
模拟。
如图2所示，模拟装置10的运算部12进行的处理可以分类为前处理和主处理。
前处理为依次进行模拟用数据生成处理(步骤S101)和参数计算处理(步骤S102)
的处理。
在步骤S101中进行的模拟用数据生成处理为基于输入的结构指定数据，生成“用
于将模拟对象空间内的各建筑物的壁面、地面视为多个面积要素，并将模拟对象空间内的
各树木的树冠视为一个以上的具有透过性的体积要素的模拟用数据”的处理。
利用该模拟用数据生成处理生成的模拟用数据只要是可得知后述的形态系数的
计算所必要的信息(各面积/体积要素的形状、在模拟对象空间内的位置等)即可。因此，模
拟用数据可以为如下数据：例如对于每个面积/体积要素，包含“由序列号、对应的模拟对象
空间内的计算网格的坐标编号、面积要素的朝向/表示要素为体积要素的标志、表示是否为
树冠的标志等构成的数据”。
在步骤S102中进行的参数计算处理为如下的处理：基于由模拟用数据生成处理生
成的模拟用数据，计算在主处理中使用的各种参数。
在该参数计算处理时算出的参数具有与各2要素(面积/体积要素)相关的形态系
数F、各体积要素k的有效表面积＜Aeff＞k、各要素i的天空因数(skyfactor)ωi、各面积要
素i的阴影率Di、各体积要素k的有效阴影率Deffk。
首先，针对与参数计算处理时算出的各2要素相关的形态系数F进行说明。
在参数计算处理时，对于两个面积要素i、j的每个组合，计算“从面积要素i观察面
积要素j的形态系数Fij”和“从面积要素j观察面积要素i的形态系数Fji”。另外，对于面积要
素i和体积要素k的每个组合，计算“从面积要素i观察体积要素k的形态系数Fik”和“从体积
要素k观察面积要素i的形态系数Fki”。进而，对于两个体积要素k和l的每个组合，计算“从体
积要素k观察体积要素l的形态系数Fkl”和“从体积要素l观察体积要素k的形态系数Flk”。
在参数计算处理时算出的“从面积要素i观察面积要素j的形态系数Fij”、“从面积
要素j观察面积要素i的形态系数Fji”分别是由以下的式(1)、式(2)定义的值。
[数学式1]
F i j = 1 A i ∫ A i ∫ A j T i j cosβ i cosβ j πr 2 dA i dA j ... ( 1 ) ]]>
F j i = 1 A j ∫ A i ∫ A j T i j cosβ i cosβ j πr 2 dA i dA j ... ( 2 ) ]]>
在上式中，Ai、Aj分别为面积要素i的面积、面积要素j的面积。βi、βj如图3中示意性
所示，分别为连接面积要素i内的微小面dAi和面积要素j内的微小面dAj的直线与微小面dAi
的法线所成的角度、以及该直线与微小面dAj的法线所成的角度。另外，r为微小面dAi和微小
面dAj之间的距离。
Tij为面积要素i与面积要素j之间的透过率。Tij使用2微小面dAi和dAj之间的光学
厚度τij并由下式计算。
[数学式2]
Tij＝exp(-τij)
另外，在面积要素i与面积要素j之间(微小面dAi与微小面dAj之间)分布有树冠的
情况下，光学厚度τij可使用树冠的分散系数(dispersioncoefficient)kext和叶面积密度a
并由下式计算。
[数学式3]
τ i j = ∫ 0 r k e x t a d r ]]>
此外，关于形态系数的具体的计算步骤将在后文中叙述，不过由上述的式(1)、式
(2)定义的形态系数满足倒易律。即，在面积Ai、面积Aj、形态系数Fij和形态系数Fji之间存在
以下的关系。
[数学式4]
AiFij＝AjFji
因此，可以根据由式(1)算出的Fij、Ai和Aj来计算Fji，并且也可以根据由式(2)算出
的Fji、Ai和Aj来计算Fij。
在参数计算处理时算出的“从面积要素i观察体积要素k的形态系数Fik”为由以下
的式(3)定义的值。
[数学式5]
F i k = 1 A i ∫ A i ∫ A k T i k cosβ i πr 2 dA i dA k ← i e f f ... ( 3 ) ]]>
该式(3)中的βi如图4中示意性所示，为连接面积要素i的微小面dAi和体积要素j内
的微小投影面dAk←i的直线与微小面dAi的法线所成的角度。另外，r为微小面dAi和微小投影
面dAk←i之间的距离。
Aeffk←i为考虑了体积要素k自身的遮挡率的、从面积要素i的方向观察的体积要素k
的有效面积。该Aeffk←i由下式计算。
[数学式6]
A k ← i e f f = ∫ A k dA k ← i e f f = ∫ A k [ 1 - exp ( - Δτ k ← i ) ] dA k ← i ... ( 4 ) ]]>
在此，Δτk←i为与微小投影面dAk←i(参照图4)垂直的方向的体积要素k的光学厚
度。在体积要素k为树木的情况下，Δτk←i使用分散系数kext、叶面积密度a、以及与dAk←i垂直
的方向的体积要素k的几何学厚度Δsk←i，并由下式计算。
[数学式7]
Δτk←i＝kextaΔsk←i
重点在于，上述的式(3)如果使用式(4)，则会成为可如下表述的式子。在参数处理
时，按照以下的式(5)，计算从面积要素i观察体积要素k的形态系数Fik。
[数学式8]
F i k = 1 A i ∫ A i ∫ A k T i k cosβ i πr 2 [ 1 - exp ( - Δτ k ← i ) ] dA i dA k ← i ... ( 5 ) ]]>
在参数计算处理时算出的“从体积要素k观察面积要素i的形态系数Fki”为由下式
定义的值。
[数学式9]
F k i = 1 A k ← i e f f ∫ A i ∫ A k T i k cosβ i πr 2 [ 1 - exp ( - Δτ k ← i ) ] dA i dA k ← i ... ( 6 ) ]]>
换言之，计算满足由下式表示的倒易律的值，作为从体积要素k观察面积要素i的
形态系数Fki。
[数学式10]
A i F i k = A k ← i e f f F k i ]]>
在参数计算处理时算出的“从体积要素k观察体积要素l的形态系数Fkl”为由下式
定义的值。
[数学式11]
F k l = 1 A k ← l e f f ∫ A l ∫ A k T k l πr 2 [ 1 - exp ( - Δτ l ← k ) ] [ 1 - exp ( - Δτ k ← l ) ] dA l ← k dA k ← l ... ( 7 ) ]]>
即，该形态系数Fkl为将体积要素k的遮挡率(“1-exp(-Δτk←l)”)、体积要素l的遮
挡率(“1-exp(-Δτl←k)”)和体积要素k、l之间的透过率Tkl纳入考虑中的形态系数。
此外，与面积要素的形态系数相同，式(7)表示的与体积要素相关的形态系数满足
下式的倒易律。
[数学式12]
A k ← l e f f F k l = A l ← k e f f F l k ]]>
因此，与其他形态系数相同，可以按照式(7)计算2体积要素之间的形态系数双方，
也可以按照式(7)计算一方的形态系数，并根据该该一方的形态系数的计算结果计算另一
方的形态系数。
在参数计算处理时，所述的各形态系数通过蒙特卡罗方法计算。
即，在计算形态系数时，根据0～1的范围的均匀随机数Rθ、并由下式计算μ、并
且基于计算结果，将生成以下所示内容的单位向量n的处理，重复与必要的形态系数的精度
相应的次数。
[数学式13]
n = ( μ c o s φ , μ s i n φ , 1 - μ 2 ) ]]>
μ2＝Rθ
φ＝2πRφ
换言之，在计算形态系数时，进行如下处理：生成按照朗伯余弦定律从面积要素
(或者体积要素)射出的大量的光子的行进方向单位向量n。
而且，在生成的各单位向量n的方向上从要素i射出具有相等的能量W0的光子时，
通过对此时射入要素j的各光子p的能量Wp进行累计，从而求出形态系数Fij。
具体地说，在从面积要素i射出的光子的个数为N、向面积要素j射入的光子的个数
为n的情况下，由下式计算从面积要素i观察面积要素j的形态系数Fij。
[数学式14]
F i j = Σ p = 1 n W p NW 0 ]]>
此外，在只有建筑物壁面等完全遮挡性的面积要素的情况下，由于向面积要素j射
入的光子的能量Wp与W0相等，因此形态系数可以由Fij＝n/N计算。在空间中分布有树木之类
的辐射透过性的要素的情况下，向面积要素j射入的光子的能量Wp在从面积要素i到达面积
要素j的期间衰减。因此，通过使光子所具有的能量衰减而考虑了能量Wp的衰减对于形态系
数的影响。即，由下式计算Wp。
[数学式15]
Wp＝Tij，pW0
在此，Tij,p为沿着光子p的路径的透过率。此外，可认为光子在透过性的体积要素
内存在被遮挡的概率，即使减少向面积要素j射入的光子的数量，也可以求出考虑了衰减影
响的形态系数。
另外，从面积要素i观察体积要素k的形态系数Fik为由式(7)定义的值。因此，从面
积要素i观察体积要素k的形态系数Fik由下式计算。
[数学式16]
F i k = Σ p = 1 n W p [ 1 - exp ( - Δτ k ← i , p ) ] NW 0 ]]>
在此，Δτk←i,p为沿着光子p的路径的、在体积要素k内部的光学厚度。N为从面积要
素i射出的光子数，n为向体积要素k射入的光子数。
在计算从体积要素k观察面积要素i的形态系数Fki时，假设来自体积要素的表面上
的各点的光子等向射出。即，从体积要素的表面上的各点按照朗伯余弦定律虚拟地射出大
量的光子。
此外，从体积要素k观察面积要素i的形态系数Fki由下式计算。
[数学式17]
A k ← i e f f F k i = Σ p = 1 m W p [ 1 - exp ( - Δτ k ← i , p ) ] ( M / S k ) W 0 ]]>
在此，通过从光子p的射出点向与光子p的行进方向相反的方向回溯来计算Δ
τk←i,p，Δτk←i,p为在体积要素k内部的光学厚度。另外，Sk为体积要素k的表面积，M为从体积
要素k射出的光子数，m为向面积要素i射入的光子数。
同样，从体积要素k观察体积要素l的形态系数Fkl由下式计算。
[数学式18]
A k ← l e f f F k l = Σ p = 1 m W p [ 1 - exp ( - Δτ k ← l , p ) ] [ 1 - exp ( - Δτ l ← k , p ) ] ( M / S k ) W 0 ]]>
此外，该式中的m为由于从体积要素k射出M个光子而向体积要素l射入的光子数。
以下，针对利用参数计算处理而算出的各体积要素k的有效表面积＜Aeff＞k、各要
素i的天空因数ωi、各面积要素i的阴影率Di、各体积要素k的有效阴影率Deffk进行说明。
在参数计算处理时算出的各体积要素k的有效表面积＜Aeff＞k为由下式定义的
值。
[数学式19]
< A e f f > k = Σ i = 1 m A k ← i e f f F k i ]]>
在该式中，m为在体积要素k的周围存在(可从体积要素k观察的)的要素(面积要素
或者体积要素)的总数，i(＝1～m)为在体积要素k的周围存在的面积要素或者体积要素的
要素编号。
在参数计算处理时，该有效表面积＜Aeff＞k由下式计算。
[数学式20]
< A e f f > k = S k M Σ p = 1 M [ 1 - exp ( - Δτ k ← p ) ] ]]>
即，有效表面积＜Aeff＞k也通过蒙特卡罗方法计算。
要素(面积要素/体积要素)i的天空因数ωi为与从要素i观察天空的形态系数相
当的值。天空因数ωi按照与从要素i观察面积要素的形态系数相同的步骤算出。
通过对从面积要素i向太阳侧射出的光子p射入其他要素时损失的能量ΔWp进行
累计而求出面积要素i的阴影率Di。更具体地说，阴影率Di使用蒙特卡罗方法并由下式计算。
[数学式21]
D i = Σ p = 1 N ΔW p NW 0 ]]>
此外，在上式中，N为从面积要素i射出的光子数。
同样，体积要素k的有效阴影率Deffk使用蒙特卡罗方法并由下式计算。
[数学式22]
D k e f f = 1 - Σ p = 1 M [ 1 - exp ( - Δτ k ← p ) ] c o s θ M / S k + Σ p = 1 M [ 1 - exp ( - Δτ k ← p ) ] cosθΔW p ( M / S k ) W 0 ]]>
所述的各面积要素i的阴影率Di、各体积要素k的有效阴影率Deffk为其值根据太阳
的位置而变化的参数。因此，当进行模拟的时间范围较长的情况下，在参数计算处理时，计
算该时间范围内的各时刻的阴影率Di和有效阴影率Deffk。
以下，对于运算部12进行的主处理(图2)的内容进行说明。
运算部12进行的主处理为如下的处理：将辐射热通量等计算处理(步骤S201)和温
度计算处理(步骤S202)重复全时间步数Nt的次数。此外，全时间步数Nt可以基于模拟对象时
刻范围和时间步长Δt决定。通过在计算条件设定文件内设定全时间步数Nt或模拟对象时
刻范围和时间步长Δt的信息，或者使用输入装置20进行输入，来进行全时间步数Nt的指
定。
在步骤S201中进行的辐射热通量等计算处理为如下的处理：使用在参数计算处理
中算出的参数(形态系数等)来计算从各要素i射出的短波辐射(可见光线)的辐射通量GS,i
[W/m2]和长波辐射(红外线)的辐射通量GL,i[W/m2]，并使用计算出的辐射通量计算由各要素
i吸收的、与短波辐射相关的净辐射热RS,i[W]和与长波辐射相关的净辐射热RL,i[W]。
具体地说，关于从各要素i射出的辐射通量GS,i、GL,i[W/m2]，下式分别成立。
[数学式23]
< A e f f > i G S , i = α S , i { A i ← S o l a r e f f S d i r e c t , i + A i ← s k y e f f S d i f f u s e , i + Σ j = 1 n A i ← j e f f F i j G S , j } ... ( 8 ) ]]>
< A e f f > i G L , i = < A e f f > i ϵ i B ( T i ) + α L , i { A i ← s k y e f f L i + Σ j = 1 n A i ← j e f f F i j G L , j } ... ( 9 ) ]]>
此外，式(8)、式(9)中的n为面积要素和体积要素的总数。<Aeff>i在要素i为体积要
素的情况下是体积要素i的有效表面积，在要素i为面积要素的情况下是面积要素i的面积。
αS,i、αL,i分别是要素i的与短波辐射、长波辐射相关的反射率，εi为要素i的射出
率。Sdirect,i为射入要素i的来自太阳的直达短波辐射通量，Sdiffuse,i为射入要素i的大气散射
短波辐射的辐射通量。Li为射入要素i的大气长波辐射的辐射通量，Aeffi←Solar、Aeffi←sky分别
为要素i的太阳方向、天空方向的有效面积。
B(Ti)为从要素i利用热辐射射出的辐射通量。在仅计算所述的GL,j作为长波辐射
的辐射通量时(在不按照波长范围计算与长波辐射相关的辐射通量时)，B(Ti)使用斯特藩-
玻尔兹曼常数σ并由下式计算。
[数学式24]
B(T)＝σT4
Sdirect,i、Sdiffuse,i使用由参数计算处理计算的天空因数ωi、阴影率Di并由下式计
算。
[数学式25]
S d i r e c t , i = c d i r e c t m a x [ 0 , S · n i ] S z ( 1 - D i ) S ↓ ]]>
Sdiffuse，i＝cdiffuseωiS↓
在此，S↓为相对于水平面向下射入的太阳辐射通量，S(＝(Sx,Sy,Sz))为太阳方位
向量。ni为面积要素i的单位法线向量，cdirect、cdiffuse为用于直散分离的系数。
此外，由上式算出的Sdirect,i为面积要素i的Sdirect,i。体积要素i的Sdirect,i使用在参
数计算处理中算出的有效阴影率Deffi并由下式计算。
[数学式26]
A i ← s o l a r e f f S d i r e c t , i = c d i r e c t ( 1 - D i e f f ) S ↓ ]]>
另外，由于射出率εi与要素i的吸收率相等，所以可使用“1-αL,i”作为射出率εi。
[数学式27]
εi＝1-αL，i
上述的式(8)和式(9)针对各个i＝1～n分别成立。即，以下的两个线性矩阵方程式
成立。
[数学式28]
D i j = A i ← j e f f F i j = A j ← i e f f F j i = D j i ]]>
在辐射通量等计算处理时，对这些线性矩阵方程式求解，由此计算从各要素i射出
的辐射通量GS,i、GL,i。
此外，由各要素i吸收的与短波辐射相关的净辐射热RS,i[W]和与长波辐射相关的
净辐射热RL,i[W]分别由以下的式(10)、式(11)计算。
[数学式29]
R S , i = A i ← S o l a r e f f S d i r e c t , i + A i ← s k y e f f S d i f f u s e , i + Σ j = 1 n A i ← j e f f F i j G S , j - < A e f f > i G S , i = < A e f f > i { G S , i / α S , i - G S , i } ... ( 10 ) ]]>
R L , i = A i ← s k y e f f L i + Σ j = 1 n A i ← j e f f F i j G L , j - < A e f f > i G L , i = < A e f f > i [ { G L , i - ϵ i B ( T i ) } / α L , i - G L , i ] ... ( 11 ) ]]>
温度计算处理(图2；步骤S202)为如下的处理：根据在辐射热通量等计算处理中计
算出的辐射热，计算模拟对象空间内的各部分的温度。除了采用将树冠视为具有透过性的
体积要素的形态系数来计算辐射热以外，该温度计算处理时的各面积要素的温度计算步骤
与一般的计算步骤相同。因此，以下仅对树冠(体积要素)的温度的计算步骤进行说明。
如图5中示意性所示，在作为树木的树冠的一部分的体积要素中，流入有短波辐射
的辐射热通量(radiationheatflux)RS和长波辐射的辐射热通量RL，从该体积要素流出显
热通量H和潜热通量LE。
因此，与作为树冠的体积要素i相关的热平衡由下式表示。
[数学式30]
Ca i V i dT l e a f , i d t = R S , i + R L , i - H i - LE i ]]>
在此，Tleaf,i为要素i内的叶的表面温度[K]，ai为要素i内的叶的面积密度[m2/m3]。
Vi为要素i的体积[m3]，C为每单位叶面积的叶的热容量[J/K/m2]。RS,i、RL,i分别为由叶吸收
的短波辐射的净辐射热(辐射热通量的强度)[W]、长波辐射的净辐射热[W]，L为蒸发潜热
[J/kg]。
Hi为从叶向大气释放的显热输送量(显热通量的强度)[W]，Ei为从叶向大气蒸腾的
水蒸气量[kg/s]。
从叶向大气释放的显热输送量Hi和从叶向大气蒸腾的水蒸气量Ei由下式计算。
[数学式31]
Hi＝aiViKh(Tleaf，i-Tair，i)
Ei＝aiViβKw(fleaf，i-fair，i)
在此，Tair,i为要素i内的大气温度[K]，fair,i为体积要素i内的大气中的水蒸气分
压[Pa]，fleaf,i为体积要素i内的叶表面的饱和水蒸气分压[Pa]，Kh为对流热传递系数[W/m2/
K]，Kw为对流水蒸气输送系数[kg/s/m2/Pa]，β为蒸发效率。
在温度计算处理时，使用时间步n的叶表面温度和热通量，来计算时间步长Δt后
的时间步n+1的叶表面温度Tleaf,i。
具体地说，从时间步n到时间步n+1之间的叶表面温度的变化量ΔTleaf,i考虑了由
于时间Δt的经过而引起的叶表面温度的变化所产生的净长波辐射、显热输送量、蒸腾量的
变化，并由下式表示。
[数学式32]
Ca i V i ΔT l e a f , i Δ t = R S , i + ( R L , i + ∂ R L , i ∂ T l e a f , i ΔT l e a f , i ) - ( H i + ∂ H i ∂ T l e a f , i ΔT l e a f , i ) - L ( E i + ∂ E i ∂ T l e a f , i ΔT l e a f , i ) ]]>
因此，可以由下式求出叶表面温度的变化量ΔTleaf,i。
[数学式33]
ΔT l e a f , i = R S , i + R L , i - H i - LE i Ca i V i Δ t - ∂ R L , i ∂ T l e a f , i + ∂ H i ∂ T l e a f , i + L ∂ E i ∂ T l e a f , i ]]>
在温度计算处理时，由该式求出叶表面温度的变化量ΔTleaf,i，然后由下式计算时
间步n+1的叶表面温度Tleaf,i(n+1)。
[数学式34]
T l e a f , i ( n + 1 ) = T l e a f , i ( n ) + ΔT l e a f , i ]]>
温度计算处理在各部分的温度的计算和计算出的各部分的温度的输出(在本实施
方式中，向存储部14存储)完成时结束。此外，在指定的次数的处理尚未完成的情况下，再次
开始辐射热通量等计算处理，在指定的次数的处理完成时，结束主处理。
如上所述，本实施方式的模拟装置10将各树冠视为一个以上的具有透过性的体积
要素，并且计算如下的形态系数作为与一个体积要素和一个面积要素相关的形态系数(参
照式(5)、式(6))，所述形态系数减少的量与透过该一个体积要素的辐射热量对应。另外，模
拟装置10计算如下的形态系数作为与两个体积要素相关的形态系数(参照式(7))，所述形
态系数减少的量与透过该两个体积要素的辐射热量对应。而且，模拟装置10使用计算出的
各形态系数计算在各要素之间交换的辐射热量。因此，根据模拟装置10，能够以无需另外计
算树干内部的状态的方式(换句话说，以低计算成本)，良好地模拟含有树冠的三维空间的
辐射热输送现象。
《变形方式》
上述实施方式的模拟装置10可进行各种变形。
例如，可将模拟装置10变形为如下装置：计算不考虑2要素之间的透过率T的形态
系数，并在计算辐射通量时等考虑2要素之间的透过率。但是，在计算形态系数时考虑2要素
之间的透过率T的方式能够得出更为准确的结果，且计算成本较低。因此，优选采用所述的
形态系数。
另外，树干的体积热容量Ca通常极小，因此可以将模拟装置10变形为以Ca＝0来计
算叶表面温度的装置。另外，还可以将模拟装置10变形为使用定义式的解析解作为一部分
或者全部的形态系数的装置。
模拟装置10为通过隐式欧拉法计算叶温度的装置，不过也可以将模拟装置10变形
为通过显式算法计算叶温度的装置或通过2次精度的Crank-Nicholson法计算叶温度的装
置。不过，隐式算法更容易得出准确的值，因此优选叶温度的计算法采用隐式算法。