CN105765585A - 关于辐射热输送现象的模拟装置、模拟方法和模拟程序 - Google Patents

Abstract

本发明提供关于辐射热输送现象的模拟装置、模拟方法和模拟程序。模拟装置能以低计算成本良好地模拟含有树冠的三维空间的辐射热输送现象。模拟装置包括：形态系数计算单元，计算与由多个面积要素和多个体积要素定义的虚拟三维空间内的各2要素相关的形态系数，并且计算如下的形态系数作为与包括一个或者两个体积要素的2要素相关的形态系数，该形态系数减少的值与透过该一个或者两个体积要素的辐射热量对应；以及辐射热量计算单元，使用由所述形态系数计算单元计算出的各形态系数，计算在各要素之间交换的辐射热量，以将存在于所述三维空间内的多颗树木的树冠视为所述多个体积要素的方式，由所述多个面积要素和所述多个体积要素定义所述三维空间。

Description

关于辐射热输送现象的模拟装置、模拟方法和模拟程序

技术领域

本发明涉及用于模拟辐射热输送现象的模拟装置、模拟方法和模拟程序。

背景技术

将城市部的气温高于周边部的气温的现象称为热岛现象。作为用于缓解热岛现象的对策，着眼于路旁树的配置、绿地的配备。因此，研发出考虑路旁树的配置等而模拟城市部的辐射热输送的各种模拟技术(参照专利文献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(LiquidCrystalDisplay)、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的有效表面积＜A^eff＞_k、各要素i的天空因数(skyfactor)ω_i、各面积要素i的阴影率D_i、各体积要素k的有效阴影率D^eff _k。

首先，针对与参数计算处理时算出的各2要素相关的形态系数F进行说明。

在参数计算处理时，对于两个面积要素i、j的每个组合，计算“从面积要素i观察面积要素j的形态系数F_ij”和“从面积要素j观察面积要素i的形态系数F_ji”。另外，对于面积要素i和体积要素k的每个组合，计算“从面积要素i观察体积要素k的形态系数F_ik”和“从体积要素k观察面积要素i的形态系数F_ki”。进而，对于两个体积要素k和l的每个组合，计算“从体积要素k观察体积要素l的形态系数F_kl”和“从体积要素l观察体积要素k的形态系数F_lk”。

在参数计算处理时算出的“从面积要素i观察面积要素j的形态系数F_ij”、“从面积要素j观察面积要素i的形态系数F_ji”分别是由以下的式(1)、式(2)定义的值。

[数学式1]

$F_{ij}=\frac{1}{A_i}{\∫}_{Ai}{\∫}_{Aj}\frac{T_{ij}{\mathrm{cos\β}}_i{\mathrm{cos\β}}_j}{{\mathrm{\πr}}^2}{\mathrm{dA}}_i{\mathrm{dA}}_j\mathrm{...}\left(1\right)$

$F_{ji}=\frac{1}{A_j}{\∫}_{Ai}{\∫}_{Aj}\frac{T_{ij}{\mathrm{cos\β}}_i{\mathrm{cos\β}}_j}{{\mathrm{\πr}}^2}{\mathrm{dA}}_i{\mathrm{dA}}_j\mathrm{...}\left(2\right)$

在上式中，A_i、A_j分别为面积要素i的面积、面积要素j的面积。β_i、β_j如图3中示意性所示，分别为连接面积要素i内的微小面dA_i和面积要素j内的微小面dA_j的直线与微小面dA_i的法线所成的角度、以及该直线与微小面dA_j的法线所成的角度。另外，r为微小面dA_i和微小面dA_j之间的距离。

T_ij为面积要素i与面积要素j之间的透过率。T_ij使用2微小面dA_i和dA_j之间的光学厚度τ_ij并由下式计算。

[数学式2]

T_ij＝exp(-τ_ij)

另外，在面积要素i与面积要素j之间(微小面dA_i与微小面dA_j之间)分布有树冠的情况下，光学厚度τ_ij可使用树冠的分散系数(dispersioncoefficient)k_ext和叶面积密度a并由下式计算。

[数学式3]

${\mathrm{\τ}}_{ij}={\∫}_0^r{k}_{ext}adr$

此外，关于形态系数的具体的计算步骤将在后文中叙述，不过由上述的式(1)、式(2)定义的形态系数满足倒易律。即，在面积A_i、面积A_j、形态系数F_ij和形态系数F_ji之间存在以下的关系。

[数学式4]

A_iF_ij＝A_jF_ji

因此，可以根据由式(1)算出的F_ij、A_i和A_j来计算F_ji，并且也可以根据由式(2)算出的F_ji、A_i和A_j来计算F_ij。

在参数计算处理时算出的“从面积要素i观察体积要素k的形态系数F_ik”为由以下的式(3)定义的值。

[数学式5]

$F_{ik}=\frac{1}{A_i}{\∫}_{Ai}{\∫}_{Ak}\frac{T_{ik}{\mathrm{cos\β}}_i}{{\mathrm{\πr}}^2}{\mathrm{dA}}_i{\mathrm{dA}}_{k\←i}^{eff}\mathrm{...}\left(3\right)$

该式(3)中的β_i如图4中示意性所示，为连接面积要素i的微小面dA_i和体积要素j内的微小投影面dA_k←i的直线与微小面dA_i的法线所成的角度。另外，r为微小面dA_i和微小投影面dA_k←i之间的距离。

A^eff _k←i为考虑了体积要素k自身的遮挡率的、从面积要素i的方向观察的体积要素k的有效面积。该A^eff _k←i由下式计算。

[数学式6]

$A_{k\←i}^{eff}={\∫}_{Ak}{\mathrm{dA}}_{k\←i}^{eff}={\∫}_{Ak}\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←i})\]{\mathrm{dA}}_{k\←i}\mathrm{...}\left(4\right)$

在此，Δτ_k←i为与微小投影面dA_k←i(参照图4)垂直的方向的体积要素k的光学厚度。在体积要素k为树木的情况下，Δτ_k←i使用分散系数k_ext、叶面积密度a、以及与dA_k←i垂直的方向的体积要素k的几何学厚度Δs_k←i，并由下式计算。

[数学式7]

Δτ_k←i＝k_extaΔs_k←i

重点在于，上述的式(3)如果使用式(4)，则会成为可如下表述的式子。在参数处理时，按照以下的式(5)，计算从面积要素i观察体积要素k的形态系数F_ik。

[数学式8]

$F_{ik}=\frac{1}{A_i}{\∫}_{Ai}{\∫}_{Ak}\frac{T_{ik}{\mathrm{cos\β}}_i}{{\mathrm{\πr}}^2}\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←i})\]{\mathrm{dA}}_i{\mathrm{dA}}_{k\←i}\mathrm{...}\left(5\right)$

在参数计算处理时算出的“从体积要素k观察面积要素i的形态系数F_ki”为由下式定义的值。

[数学式9]

$F_{ki}=\frac{1}{A_{k\←i}^{eff}}{\∫}_{Ai}{\∫}_{Ak}\frac{T_{ik}{\mathrm{cos\β}}_i}{{\mathrm{\πr}}^2}\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←i})\]{\mathrm{dA}}_i{\mathrm{dA}}_{k\←i}\mathrm{...}\left(6\right)$

换言之，计算满足由下式表示的倒易律的值，作为从体积要素k观察面积要素i的形态系数F_ki。

[数学式10]

$A_i{F}_{ik}=A_{k\←i}^{eff}{F}_{ki}$

在参数计算处理时算出的“从体积要素k观察体积要素l的形态系数F_kl”为由下式定义的值。

[数学式11]

$F_{kl}=\frac{1}{A_{k\←l}^{eff}}{\∫}_{Al}{\∫}_{Ak}\frac{T_{kl}}{{\mathrm{\πr}}^2}\[1-\exp (-{\mathrm{\Δ\τ}}_{l\←k})\]\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←l})\]{\mathrm{dA}}_{l\←k}{\mathrm{dA}}_{k\←l}\mathrm{...}\left(7\right)$

即，该形态系数F_kl为将体积要素k的遮挡率(“1-exp(-Δτ_k←l)”)、体积要素l的遮挡率(“1-exp(-Δτ_l←k)”)和体积要素k、l之间的透过率T_kl纳入考虑中的形态系数。

此外，与面积要素的形态系数相同，式(7)表示的与体积要素相关的形态系数满足下式的倒易律。

[数学式12]

$A_{k\←l}^{eff}{F}_{kl}=A_{l\←k}^{eff}{F}_{lk}$

因此，与其他形态系数相同，可以按照式(7)计算2体积要素之间的形态系数双方，也可以按照式(7)计算一方的形态系数，并根据该该一方的形态系数的计算结果计算另一方的形态系数。

在参数计算处理时，所述的各形态系数通过蒙特卡罗方法计算。

即，在计算形态系数时，根据0～1的范围的均匀随机数R_θ、并由下式计算μ、并且基于计算结果，将生成以下所示内容的单位向量n的处理，重复与必要的形态系数的精度相应的次数。

[数学式13]

$n=(\mathrm{\μ}cos\mathrm{\φ},\mathrm{\μ}sin\mathrm{\φ},\sqrt{1-{\mathrm{\μ}}^2})$

μ²＝R_θ

φ＝2πR_φ

换言之，在计算形态系数时，进行如下处理：生成按照朗伯余弦定律从面积要素(或者体积要素)射出的大量的光子的行进方向单位向量n。

而且，在生成的各单位向量n的方向上从要素i射出具有相等的能量W₀的光子时，通过对此时射入要素j的各光子p的能量W_p进行累计，从而求出形态系数F_ij。

具体地说，在从面积要素i射出的光子的个数为N、向面积要素j射入的光子的个数为n的情况下，由下式计算从面积要素i观察面积要素j的形态系数F_ij。

[数学式14]

$F_{ij}=\frac{\underset{p=1}{\overset{n}{\Σ}}{W}_p}{{\mathrm{NW}}_0}$

此外，在只有建筑物壁面等完全遮挡性的面积要素的情况下，由于向面积要素j射入的光子的能量W_p与W₀相等，因此形态系数可以由F_ij＝n/N计算。在空间中分布有树木之类的辐射透过性的要素的情况下，向面积要素j射入的光子的能量W_p在从面积要素i到达面积要素j的期间衰减。因此，通过使光子所具有的能量衰减而考虑了能量W_p的衰减对于形态系数的影响。即，由下式计算W_p。

[数学式15]

W_p＝T_ij，pW₀

在此，T_ij,p为沿着光子p的路径的透过率。此外，可认为光子在透过性的体积要素内存在被遮挡的概率，即使减少向面积要素j射入的光子的数量，也可以求出考虑了衰减影响的形态系数。

另外，从面积要素i观察体积要素k的形态系数F_ik为由式(7)定义的值。因此，从面积要素i观察体积要素k的形态系数F_ik由下式计算。

[数学式16]

$F_{ik}=\frac{\underset{p=1}{\overset{n}{\Σ}}{W}_p\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←i,p})\]}{{\mathrm{NW}}_0}$

在此，Δτ_k←i,p为沿着光子p的路径的、在体积要素k内部的光学厚度。N为从面积要素i射出的光子数，n为向体积要素k射入的光子数。

在计算从体积要素k观察面积要素i的形态系数F_ki时，假设来自体积要素的表面上的各点的光子等向射出。即，从体积要素的表面上的各点按照朗伯余弦定律虚拟地射出大量的光子。

此外，从体积要素k观察面积要素i的形态系数F_ki由下式计算。

[数学式17]

$A_{k\←i}^{eff}{F}_{ki}=\frac{\underset{p=1}{\overset{m}{\Σ}}{W}_p\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←i,p})\]}{(M/S_k)W_0}$

在此，通过从光子p的射出点向与光子p的行进方向相反的方向回溯来计算Δτ_k←i,p，Δτ_k←i,p为在体积要素k内部的光学厚度。另外，S_k为体积要素k的表面积，M为从体积要素k射出的光子数，m为向面积要素i射入的光子数。

同样，从体积要素k观察体积要素l的形态系数F_kl由下式计算。

[数学式18]

$A_{k\←l}^{eff}{F}_{kl}=\frac{\underset{p=1}{\overset{m}{\Σ}}{W}_p\[1-\exp (-{\mathrm{\Δ\τ}}_{k\←l,p})\]\[1-\exp (-{\mathrm{\Δ\τ}}_{l\←k,p})\]}{(M/S_k)W_0}$

此外，该式中的m为由于从体积要素k射出M个光子而向体积要素l射入的光子数。

以下，针对利用参数计算处理而算出的各体积要素k的有效表面积＜A^eff＞_k、各要素i的天空因数ω_i、各面积要素i的阴影率D_i、各体积要素k的有效阴影率D^eff _k进行说明。

在参数计算处理时算出的各体积要素k的有效表面积＜A^eff＞_k为由下式定义的值。

[数学式19]

$<A^{eff}{>}_k=\underset{i=1}{\overset{m}{\&Sigma;}}{A}_{k\&LeftArrow;i}^{eff}{F}_{ki}$

在该式中，m为在体积要素k的周围存在(可从体积要素k观察的)的要素(面积要素或者体积要素)的总数，i(＝1～m)为在体积要素k的周围存在的面积要素或者体积要素的要素编号。

在参数计算处理时，该有效表面积＜A^eff＞_k由下式计算。

[数学式20]

$<A^{eff}{>}_k=\frac{S_k}{M}\underset{p=1}{\overset{M}{\&Sigma;}}\&lsqb;1-\exp (-{\mathrm{\&Delta;\&tau;}}_{k\&LeftArrow;p})\&rsqb;$

即，有效表面积＜A^eff＞_k也通过蒙特卡罗方法计算。

要素(面积要素/体积要素)i的天空因数ω_i为与从要素i观察天空的形态系数相当的值。天空因数ω_i按照与从要素i观察面积要素的形态系数相同的步骤算出。

通过对从面积要素i向太阳侧射出的光子p射入其他要素时损失的能量ΔW_p进行累计而求出面积要素i的阴影率D_i。更具体地说，阴影率D_i使用蒙特卡罗方法并由下式计算。

[数学式21]

$D_i=\frac{\underset{p=1}{\overset{N}{\&Sigma;}}{\mathrm{\&Delta;W}}_p}{{\mathrm{NW}}_0}$

此外，在上式中，N为从面积要素i射出的光子数。

同样，体积要素k的有效阴影率D^eff _k使用蒙特卡罗方法并由下式计算。

[数学式22]

$D_k^{eff}=1-\frac{\underset{p=1}{\overset{M}{\&Sigma;}}\&lsqb;1-\exp (-{\mathrm{\&Delta;\&tau;}}_{k\&LeftArrow;p})\&rsqb;cos\mathrm{\&theta;}}{M/S_k}+\frac{\underset{p=1}{\overset{M}{\&Sigma;}}\&lsqb;1-\exp (-{\mathrm{\&Delta;\&tau;}}_{k\&LeftArrow;p})\&rsqb;{\mathrm{cos\&theta;\&Delta;W}}_p}{(M/S_k)W_0}$

所述的各面积要素i的阴影率D_i、各体积要素k的有效阴影率D^eff _k为其值根据太阳的位置而变化的参数。因此，当进行模拟的时间范围较长的情况下，在参数计算处理时，计算该时间范围内的各时刻的阴影率D_i和有效阴影率D^eff _k。

以下，对于运算部12进行的主处理(图2)的内容进行说明。

运算部12进行的主处理为如下的处理：将辐射热通量等计算处理(步骤S201)和温度计算处理(步骤S202)重复全时间步数N_t的次数。此外，全时间步数N_t可以基于模拟对象时刻范围和时间步长Δt决定。通过在计算条件设定文件内设定全时间步数N_t或模拟对象时刻范围和时间步长Δt的信息，或者使用输入装置20进行输入，来进行全时间步数N_t的指定。

在步骤S201中进行的辐射热通量等计算处理为如下的处理：使用在参数计算处理中算出的参数(形态系数等)来计算从各要素i射出的短波辐射(可见光线)的辐射通量G_S,i[W/m²]和长波辐射(红外线)的辐射通量G_L,i[W/m²]，并使用计算出的辐射通量计算由各要素i吸收的、与短波辐射相关的净辐射热R_S,i[W]和与长波辐射相关的净辐射热R_L,i[W]。

具体地说，关于从各要素i射出的辐射通量G_S,i、G_L,i[W/m²]，下式分别成立。

[数学式23]

$<A^{eff}{>}_i{G}_{S,i}={\mathrm{\&alpha;}}_{S,i}\left\{A_{i\&LeftArrow;Solar}^{eff}{S}_{direct,i}+A_{i\&LeftArrow;sky}^{eff}{S}_{diffuse,i}+\underset{j=1}{\overset{n}{\&Sigma;}}{A}_{i\&LeftArrow;j}^{eff}{F}_{ij}{G}_{S,j}\right\}\mathrm{...}\left(8\right)$

$<A^{eff}{>}_i{G}_{L,i}=<A^{eff}{>}_i{\mathrm{\&epsiv;}}_{i}B\left(T_i\right)+{\mathrm{\&alpha;}}_{L,i}\left\{A_{i\&LeftArrow;sky}^{eff}{L}_i+\underset{j=1}{\overset{n}{\&Sigma;}}{A}_{i\&LeftArrow;j}^{eff}{F}_{ij}{G}_{L,j}\right\}\mathrm{...}\left(9\right)$

此外，式(8)、式(9)中的n为面积要素和体积要素的总数。<A^eff>i在要素i为体积要素的情况下是体积要素i的有效表面积，在要素i为面积要素的情况下是面积要素i的面积。

α_S,i、α_L,i分别是要素i的与短波辐射、长波辐射相关的反射率，ε_i为要素i的射出率。S_direct,i为射入要素i的来自太阳的直达短波辐射通量，S_diffuse,i为射入要素i的大气散射短波辐射的辐射通量。L_i为射入要素i的大气长波辐射的辐射通量，A^eff _i←Solar、A^eff _i←sky分别为要素i的太阳方向、天空方向的有效面积。

B(T_i)为从要素i利用热辐射射出的辐射通量。在仅计算所述的G_L,j作为长波辐射的辐射通量时(在不按照波长范围计算与长波辐射相关的辐射通量时)，B(T_i)使用斯特藩-玻尔兹曼常数σ并由下式计算。

[数学式24]

B(T)＝σT⁴

S_direct,i、S_diffuse,i使用由参数计算处理计算的天空因数ω_i、阴影率D_i并由下式计算。

[数学式25]

$S_{direct,i}=c_{direct}\frac{max\&lsqb;0,S\&CenterDot;n_i\&rsqb;}{S_z}(1-D_i)S^{\&DownArrow;}$

Sd_iffuse，i＝c_diffuseω_iS^↓

在此，S^↓为相对于水平面向下射入的太阳辐射通量，S(＝(Sx,Sy,Sz))为太阳方位向量。n_i为面积要素i的单位法线向量，c_direct、c_diffuse为用于直散分离的系数。

此外，由上式算出的S_direct,i为面积要素i的S_direct,i。体积要素i的S_direct,i使用在参数计算处理中算出的有效阴影率D^eff _i并由下式计算。

[数学式26]

$A_{i\&LeftArrow;solar}^{eff}{S}_{direct,i}=c_{direct}(1-D_i^{eff})S^{\&DownArrow;}$

另外，由于射出率ε_i与要素i的吸收率相等，所以可使用“1-α_L,i”作为射出率ε_i。

[数学式27]

ε_i＝1-α_L，i

上述的式(8)和式(9)针对各个i＝1～n分别成立。即，以下的两个线性矩阵方程式成立。

[数学式28]

$D_{ij}=A_{i\&LeftArrow;j}^{eff}{F}_{ij}=A_{j\&LeftArrow;i}^{eff}{F}_{ji}=D_{ji}$

在辐射通量等计算处理时，对这些线性矩阵方程式求解，由此计算从各要素i射出的辐射通量G_S,i、G_L,i。

此外，由各要素i吸收的与短波辐射相关的净辐射热R_S,i[W]和与长波辐射相关的净辐射热R_L,i[W]分别由以下的式(10)、式(11)计算。

[数学式29]

$\begin{array}{c}{R}_{S,i}=A_{i\&LeftArrow;Solar}^{eff}{S}_{direct,i}+A_{i\&LeftArrow;sky}^{eff}{S}_{diffuse,i}+\underset{j=1}{\overset{n}{\&Sigma;}}{A}_{i\&LeftArrow;j}^{eff}{F}_{ij}{G}_{S,j}-<A^{eff}{>}_i{G}_{S,i}\\ =<A^{eff}{>}_i\left\{G_{S,i}/{\mathrm{\&alpha;}}_{S,i}-G_{S,i}\right\}\end{array}\mathrm{...}\left(10\right)$

$\begin{array}{c}{R}_{L,i}=A_{i\&LeftArrow;sky}^{eff}{L}_i+\underset{j=1}{\overset{n}{\&Sigma;}}{A}_{i\&LeftArrow;j}^{eff}{F}_{ij}{G}_{L,j}-<A^{eff}{>}_i{G}_{L,i}\\ =<A^{eff}{>}_i\&lsqb;\left\{G_{L,i}-{\mathrm{\&epsiv;}}_{i}B\left(T_i\right)\right\}/{\mathrm{\&alpha;}}_{L,i}-G_{L,i}\&rsqb;\end{array}\mathrm{...}\left(11\right)$

温度计算处理(图2；步骤S202)为如下的处理：根据在辐射热通量等计算处理中计算出的辐射热，计算模拟对象空间内的各部分的温度。除了采用将树冠视为具有透过性的体积要素的形态系数来计算辐射热以外，该温度计算处理时的各面积要素的温度计算步骤与一般的计算步骤相同。因此，以下仅对树冠(体积要素)的温度的计算步骤进行说明。

如图5中示意性所示，在作为树木的树冠的一部分的体积要素中，流入有短波辐射的辐射热通量(radiationheatflux)R_S和长波辐射的辐射热通量R_L，从该体积要素流出显热通量H和潜热通量LE。

因此，与作为树冠的体积要素i相关的热平衡由下式表示。

[数学式30]

${\mathrm{Ca}}_i{V}_i\frac{{\mathrm{dT}}_{leaf,i}}{dt}=R_{S,i}+R_{L,i}-H_i-{\mathrm{LE}}_i$

在此，T_leaf,i为要素i内的叶的表面温度[K]，a_i为要素i内的叶的面积密度[m²/m³]。V_i为要素i的体积[m³]，C为每单位叶面积的叶的热容量[J/K/m²]。R_S,i、R_L,i分别为由叶吸收的短波辐射的净辐射热(辐射热通量的强度)[W]、长波辐射的净辐射热[W]，L为蒸发潜热[J/kg]。

H_i为从叶向大气释放的显热输送量(显热通量的强度)[W]，E_i为从叶向大气蒸腾的水蒸气量[kg/s]。

从叶向大气释放的显热输送量H_i和从叶向大气蒸腾的水蒸气量E_i由下式计算。

[数学式31]

H_i＝a_iV_iK_h(T_leaf，i-T_air，i)

E_i＝a_iV_iβK_w(f_leaf，i-f_air，i)

在此，T_air,i为要素i内的大气温度[K]，f_air,i为体积要素i内的大气中的水蒸气分压[Pa]，f_leaf,i为体积要素i内的叶表面的饱和水蒸气分压[Pa]，K_h为对流热传递系数[W/m²/K]，K_w为对流水蒸气输送系数[kg/s/m²/Pa]，β为蒸发效率。

在温度计算处理时，使用时间步n的叶表面温度和热通量，来计算时间步长Δt后的时间步n+1的叶表面温度T_leaf,i。

具体地说，从时间步n到时间步n+1之间的叶表面温度的变化量ΔT_leaf,i考虑了由于时间Δt的经过而引起的叶表面温度的变化所产生的净长波辐射、显热输送量、蒸腾量的变化，并由下式表示。

[数学式32]

${\mathrm{Ca}}_i{V}_i\frac{{\mathrm{\&Delta;T}}_{leaf,i}}{\mathrm{\&Delta;}t}=R_{S,i}+\left(R_{L,i}+\frac{\&part;R_{L,i}}{\&part;T_{leaf,i}}{\mathrm{\&Delta;T}}_{leaf,i}\right)-\left(H_i+\frac{\&part;H_i}{\&part;T_{leaf,i}}{\mathrm{\&Delta;T}}_{leaf,i}\right)-L\left(E_i+\frac{\&part;E_i}{\&part;T_{leaf,i}}{\mathrm{\&Delta;T}}_{leaf,i}\right)$

因此，可以由下式求出叶表面温度的变化量ΔT_leaf,i。

[数学式33]

${\mathrm{\&Delta;T}}_{leaf,i}=\frac{R_{S,i}+R_{L,i}-H_i-{\mathrm{LE}}_i}{\frac{{\mathrm{Ca}}_i{V}_i}{\mathrm{\&Delta;}t}-\frac{\&part;R_{L,i}}{\&part;T_{leaf,i}}+\frac{\&part;H_i}{\&part;T_{leaf,i}}+L\frac{\&part;E_i}{\&part;T_{leaf,i}}}$

在温度计算处理时，由该式求出叶表面温度的变化量ΔT_leaf,i，然后由下式计算时间步n+1的叶表面温度T_leaf,i ⁽ⁿ⁺¹⁾。

[数学式34]

$T_{leaf,i}^{\left(n+1\right)}=T_{leaf,i}^{\left(n\right)}+{\mathrm{\&Delta;T}}_{leaf,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法计算叶温度的装置。不过，隐式算法更容易得出准确的值，因此优选叶温度的计算法采用隐式算法。

Claims (6)

1.一种模拟装置，用于模拟辐射热输送现象，其特征在于，包括：

形态系数计算单元，计算与由多个面积要素和多个体积要素定义的虚拟三维空间内的各2要素相关的形态系数，并且计算如下的形态系数作为与含有一个或者两个体积要素的2要素相关的形态系数，该形态系数减少的值与透过该一个或者两个体积要素的辐射热量对应；以及

辐射热量计算单元，使用由所述形态系数计算单元计算出的各形态系数，计算在各要素之间交换的辐射热量，

以将存在于所述三维空间内的多颗树木的树冠视为所述多个体积要素的方式，由所述多个面积要素和所述多个体积要素定义所述三维空间。

2.根据权利要求1所述的模拟装置，其特征在于，所述形态系数计算单元计算将各2要素之间的辐射热量的透过率纳入考虑中的形态系数，作为与各2要素相关的形态系数。

3.根据权利要求1或2所述的模拟装置，其特征在于，还包括温度计算单元，所述温度计算单元基于由所述辐射热量计算单元计算出的辐射热量计算各要素的温度。

4.根据权利要求3所述的模拟装置，其特征在于，所述温度计算单元利用隐式算法计算各要素的温度。

5.一种模拟方法，用于模拟辐射热输送现象，其特征在于，

由计算机进行如下的步骤：

形态系数计算步骤，以将多颗树木的树冠视为多个体积要素的方式，来计算与由多个面积要素和多个体积要素定义的虚拟三维空间内的各2要素相关的形态系数，并且计算如下的形态系数作为与含有一个或者两个体积要素的2要素相关的形态系数，该形态系数减少的值与透过该一个或者两个体积要素的辐射热量对应；以及

辐射热量计算步骤，使用由所述形态系数计算单元计算出的各形态系数，来计算在各要素之间交换的辐射热量。

6.一种模拟程序，用于模拟辐射热输送现象，其特征在于，

使计算机进行如下的步骤：

形态系数计算步骤，以将多颗树木的树冠视为多个体积要素的方式，来计算与由多个面积要素和多个体积要素定义的虚拟三维空间内的各2要素相关的形态系数，并且计算如下的形态系数作为与含有一个或者两个体积要素的2要素相关的形态系数，该形态系数减少的值与透过该一个或者两个体积要素的辐射热量对应；以及

辐射热量计算步骤，使用由所述形态系数计算单元计算出的各形态系数，来计算在各要素之间交换的辐射热量。