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

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

Info

Publication number
WO2015072546A1
WO2015072546A1 PCT/JP2014/080198 JP2014080198W WO2015072546A1 WO 2015072546 A1 WO2015072546 A1 WO 2015072546A1 JP 2014080198 W JP2014080198 W JP 2014080198W WO 2015072546 A1 WO2015072546 A1 WO 2015072546A1
Authority
WO
WIPO (PCT)
Prior art keywords
elements
radiant heat
form factor
simulation
calculating
Prior art date
Application number
PCT/JP2014/080198
Other languages
English (en)
French (fr)
Inventor
景吾 松田
桂子 高橋
Original Assignee
独立行政法人海洋研究開発機構
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 独立行政法人海洋研究開発機構 filed Critical 独立行政法人海洋研究開発機構
Priority to US15/036,820 priority Critical patent/US10796034B2/en
Priority to DE112014005224.8T priority patent/DE112014005224T5/de
Priority to CN201480062470.9A priority patent/CN105765585B/zh
Publication of WO2015072546A1 publication Critical patent/WO2015072546A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Definitions

  • the present invention relates to a simulation apparatus, a simulation method, and a simulation program for simulating a radiant heat transport phenomenon.
  • the existing simulation technology uses the radiation related to the canopy, such as “ignoring the scattering of the radiant heat by the canopy and considering only the attenuation and absorption of the radiant heat by the canopy”.
  • the energy transfer process is simplified. For this reason, existing simulation techniques often cannot adequately simulate the amount of absorbed radiant energy near or inside the canopy.
  • An object of the present invention is to provide a technique capable of successfully simulating a radiant heat transport phenomenon in a three-dimensional space including a tree canopy at a low calculation cost.
  • a simulation apparatus for simulating a radiant heat transport phenomenon of the present invention Form factor calculation means for calculating a form factor for each two elements in a virtual three-dimensional space defined by a plurality of area elements and a plurality of volume elements, and includes two elements including one or two volume elements A form factor calculating means for calculating a form factor having a value reduced by an amount corresponding to the amount of radiant heat transmitted through the one or two volume elements, Radiant heat amount calculating means for calculating the amount of radiant heat transferred between each element using each shape factor calculated by the form factor calculating means;
  • the apparatus is configured such that the three-dimensional space is defined by the plurality of area elements and the plurality of volume elements so that crowns of a plurality of trees existing in the three-dimensional space are handled as the plurality of volume elements. Is done.
  • the simulation apparatus treats each tree crown as one or more volume elements having permeability, and the amount of radiant heat transmitted through the one volume element as a form factor related to one volume element and one area element.
  • the form factor with the value decreased by the amount corresponding to is calculated.
  • the simulation apparatus calculates a view factor whose value is decreased by an amount corresponding to the amount of radiant heat transmitted through the two volume elements as the view factor related to the two volume elements. Then, the simulation apparatus calculates the amount of radiant heat transferred between each element using each calculated form factor. Therefore, according to the simulation apparatus of the present invention, it is possible to satisfactorily simulate the radiant heat transport phenomenon in the three-dimensional space including the tree crown without requiring a separate calculation of the state inside the trunk (that is, at a low calculation cost). become.
  • the present invention can also be realized as a simulation method having the same characteristics as the above simulation apparatus, or as a simulation program for causing an information processing apparatus (computer) to function as the simulation apparatus.
  • the present invention can also be realized as a computer-readable medium in which the simulation program is recorded.
  • the radiant heat transport phenomenon in a three-dimensional space including a tree canopy can be simulated well at a low calculation cost.
  • FIG. 1 shows a configuration of a simulation system according to the embodiment.
  • the simulation system according to the present embodiment includes a simulation device 10, an input device 20, and an output device 30.
  • the simulation apparatus 10 includes a calculation unit 12, a storage unit 14, and an interface circuit (I / F) 16.
  • the interface circuit (I / F) 16 of the simulation apparatus 10 is a circuit used by the arithmetic unit 12 to communicate with other apparatuses.
  • the storage unit 14 is a nonvolatile storage device that stores a simulation program 18.
  • the storage unit 14 is also used for storing data used by the calculation unit 12 and processing results obtained by the calculation unit 12.
  • the calculation unit 12 is a unit that combines a CPU (Central Processing Unit), a RAM (Random Access Memory), and the like.
  • the arithmetic unit 12 reads out the simulation program 18 from the storage unit 14 and executes it to perform various processes (details will be described later). Further, the calculation unit 12 functions as a form factor calculation unit, a radiant heat amount calculation unit, and a temperature calculation unit according to the present invention by executing the simulation program 18.
  • the input device 20 is a device for inputting information to the simulation device 10.
  • the input device 20 includes a pointing device such as a keyboard and a mouse.
  • the output device 30 is a display such as an LCD (Liquid Crystal Display) or a CRT (Cathode-Ray Tube) for outputting information from the simulation device 10, a printer, or the like.
  • the simulation apparatus 10 is usually realized by causing a computer (such as a vector-type parallel computer) that can perform matrix operations at high speed to execute the simulation program 18. Therefore, the input device and the output device of the computer for input / output connected to the simulation device 10 (vector type parallel computer or the like) normally function as the input device 20 and the output device 30.
  • a computer such as a vector-type parallel computer
  • the input device and the output device of the computer for input / output connected to the simulation device 10 normally function as the input device 20 and the output device 30.
  • the simulation device 10 is a device for simulating a radiant heat transport phenomenon in an urban space including a plurality of trees (hereinafter referred to as a simulation target space).
  • a calculation condition setting file in which information on the altitude and direction of the sun within the simulation target time range is set is stored in the storage unit 14. Further, terrain data, building shape data, tree distribution data, ground surface temperature data, building surface temperature data, and leaf surface temperature data are input to the simulation apparatus 10.
  • Terrain data is data indicating the topography (ground shape) of the simulation target space.
  • the building shape data is data indicating the position and shape of each building existing in the simulation target space.
  • the tree distribution data is data indicating the position and shape of each tree existing in the simulation target space and the leaf area density distribution of each tree.
  • These data input to the simulation apparatus 10 may be data that can understand the structure of the simulation target space (three-dimensional city shape and three-dimensional tree distribution). Therefore, for example, two-dimensional plane data of the ground / building height (data representing the correspondence between the ground / building height and coordinates) can be used as the terrain data and the building shape data. Further, as the tree distribution data, for example, data including two-dimensional plane data of a tree index (tree identification information) and data indicating a vertical distribution of a leaf area density of a tree identified by each tree index may be used. it can.
  • the ground surface temperature data, the building surface temperature data, and the leaf surface temperature data input to the simulation apparatus 10 respectively indicate data indicating initial values of temperatures of respective parts of the ground surface and initial values of temperatures of respective parts of the respective building surfaces. It is data which shows the initial value of data and the temperature of the leaf surface of each place.
  • simulation unit 12 the contents of processing performed by the simulation apparatus 10 (calculation unit 12) will be described.
  • structure designation data terrain data, building shape data, and tree distribution data input to the simulation apparatus 10
  • initial value data initial value data.
  • the simulation apparatus 10 basically uses the various information in the calculation condition setting file and the input initial value data to calculate the temperature of each part in the simulation target space indicated by the input structure designation data. It is a device that simulates every time.
  • the processing performed by the calculation unit 12 of the simulation apparatus 10 can be classified into preprocessing and main processing.
  • the pre-processing is processing in which simulation data generation processing (step S101) and parameter calculation processing (step S102) are performed in this order.
  • the simulation data generation process performed in step S101 is based on the input structure designation data, “the walls and ground of each building in the simulation target space are handled as a plurality of area elements, and each tree in the simulation target space is processed. This is a process of generating “simulation data” for handling a tree crown as a volume element having one or more permeability.
  • the simulation data generated by the simulation data generation process only needs to know information (such as the shape of each area / volume element and the position in the simulation target space) necessary for calculating the form factor described later. . Therefore, for example, for each area / volume element, the simulation data includes “serial number, coordinate number of the corresponding calculation grid in the simulation target space, flag indicating that the area element is directed / volume element, crown Data including a flag indicating whether or not the data is included.
  • the parameter calculation process performed in step S102 is a process for calculating various parameters used in the main process based on the simulation data generated by the simulation data generation process.
  • the parameters are calculated during the parameter calculation process, the two-factor (area / volume elements) form about factor F, the effective surface area ⁇ A eff> k for each volume element k, sky view factor omega i of each element i, the area There is a shade rate D i for element i and an effective shade rate D eff k for each volume element k.
  • the “form factor F ij that looks at the area element j from the area element i” and “the form factor F ji that looks at the area element i from the area element j” calculated during the parameter calculation processing are respectively the following (1), (2) A value defined by the equation.
  • a i and A j are the area of the area element i and the area of the area element j, respectively.
  • ⁇ i and ⁇ j are respectively a straight line connecting the minute surface dA i in the area element i and the minute surface dA j in the area element j and the minute surface dA i .
  • the angle formed by the normal line is the angle formed by the normal line of the straight line and the minute surface dA j .
  • R is the distance between the minute surface dA i and the minute surface dA j .
  • T ij is the transmittance between the area element i and the area element j.
  • T ij is calculated by the following equation using the optical thickness ⁇ ij between the two minute surfaces dA i and dA j .
  • the optical thickness ⁇ ij is determined by the extinction coefficient k ext and It is calculated by the following equation using the leaf area density a.
  • the form factor defined by the above formulas (1) and (2) satisfies the reciprocity law. That is, there is the following relationship among the area A i , the area A j , the form factor F ij, and the form factor F ji .
  • F ji is calculated from F ij calculated from equation (1), A i, and A j
  • F ij is calculated from F ji calculated from equation (2), A i, and A j.
  • the “form factor F ik when the volume element k is viewed from the area element i” calculated during the parameter calculation process is a value defined by the following equation (3).
  • ⁇ i in the equation (3) is a micro surface and a straight line connecting the micro surface dA i of the area element i and the micro projection surface dA k ⁇ i in the volume element j. It is an angle formed by the normal line of dA i .
  • R is the distance between the minute surface dA i and the minute projection surface dA k ⁇ i .
  • a eff k ⁇ i is the effective area of the volume element k viewed from the direction of the area element i, taking into account the shielding factor of the volume element k itself. This A eff k ⁇ i is calculated by the following equation.
  • ⁇ k ⁇ i is the optical thickness of the volume element k in the direction perpendicular to the minute projection surface dA k ⁇ i (see FIG. 4).
  • ⁇ k ⁇ i is expressed by using the extinction coefficient k ext , the leaf area density a and the geometric thickness ⁇ s k ⁇ i of the volume element k in the direction perpendicular to dA k ⁇ i It is calculated by the following formula.
  • the above-described expression (3) is an expression that can be expressed as follows if expression (4) is used.
  • the form factor F ik in which the volume element k is viewed from the area element i is calculated according to the following equation (5).
  • the “form factor F ki obtained by viewing the area element i from the volume element k” calculated during the parameter calculation process is a value defined by the following equation.
  • a value that satisfies the reciprocity law expressed by the following equation is calculated as the form factor F ki when the area element i is viewed from the volume element k.
  • the “form factor F kl of the volume element 1 seen from the volume element k” calculated during the parameter calculation process is a value defined by the following equation.
  • the form factor F kl is determined by the shielding factor of the volume element k (“1-exp ( ⁇ k ⁇ l )”) and the shielding factor of the volume element l (“1-exp ( ⁇ l ⁇ k )”). And the shape factor taking into account the transmittance T kl between the volume elements k and l.
  • the form factor regarding the volume element represented by the equation (7) satisfies the reciprocity law of the following equation, similarly to the form factor of the area element.
  • both form factors between two volume elements can be calculated according to equation (7), or one form factor can be calculated according to equation (7) and the other form factor It can also be calculated from the calculation result of the one form factor.
  • each of the above-described form factors is calculated by the Monte Carlo method.
  • ⁇ and ⁇ are calculated from the uniform random numbers R ⁇ and R ⁇ in the range from 0 to 1 by the following formula, and the unit vector n having the following contents is generated based on the calculation result. This process is repeated a number of times according to the accuracy of the required form factor.
  • a process of generating a unit vector n of traveling directions of a large number of photons emitted from the area element (or volume element) is performed according to Lambert's cosine law.
  • the form factor F ij when the area element i is viewed from the area element i is N when the number of photons emitted from the area element i is N and the number of photons incident on the area element j is n. Is calculated by the following equation.
  • the energy W p of the photons incident on the area element j is attenuated before reaching the area element j from the area element i. Therefore, the influence of the attenuation on the form factor is taken into account by attenuating the energy of the photons. That is, W p is calculated by the following equation.
  • T ij, p is the transmittance along the path of the photon p. Even if the number of photons incident on the area element j is reduced assuming that the photons are probabilistically shielded in the transmissive volume element, the form factor can be obtained in consideration of the influence of attenuation.
  • the shape factor F ik when the volume element k is viewed from the area element i is a value defined by the equation (7). Therefore, the form factor F ik when the volume element k is viewed from the area element i is calculated by the following equation.
  • ⁇ k ⁇ i, p is the optical thickness inside the volume element k along the path of the photon p.
  • N is the number of photons emitted from the area element i, and n is the number of photons incident on the volume element k.
  • the form factor F ki obtained by viewing the area element i from the volume element k is calculated by the following equation.
  • ⁇ k ⁇ i, p is the optical thickness inside the volume element k calculated by back-tracing from the emission point of the photon p in the direction opposite to the traveling direction of the photon p.
  • S k is the surface area of the volume element k
  • M is the number of photons emitted from the volume element k
  • m is the number of photons incident on the area element i.
  • the form factor F kl when the volume element l is viewed from the volume element k is calculated by the following equation.
  • m in this equation is the number of photons incident on the volume element l as a result of M photons being emitted from the volume element k.
  • the effective shadow factor of effective surface area ⁇ A eff> k, sky factor omega i, shadow factor D i, each volume element k of each area element i of each element i of each volume element k D eff k will be described.
  • Effective surface area ⁇ A eff> k for each volume element k calculated at the parameter calculating process is a value defined by the following equation.
  • m is the total number of elements (area element or volume element) existing around (visible from the volume element k) around the volume element k
  • i is around the volume element k. This is the element number of the existing area element or volume element.
  • this effective surface area ⁇ A eff> k is calculated by the following equation.
  • the effective surface area ⁇ A eff> k is also calculated by the Monte Carlo method.
  • the sky factor ⁇ i of an element (area element / volume element) i is a value corresponding to a form factor when the element i is viewed from the sky.
  • the sky factor ⁇ i is calculated in the same procedure as the view factor when the area element is viewed from the element i.
  • the shade ratio D i of the area element i is obtained by integrating the energy ⁇ Wp that is lost when the photon p emitted from the area element i to the sun enters the other element. More specifically, the shade ratio D i is calculated by the following equation using the Monte Carlo method.
  • N is the number of photons emitted from the area element i.
  • the effective shade rate D eff k of the volume element k is calculated by the following equation using the Monte Carlo method.
  • the main processing operation unit 12 performs includes a radiant heat flux such calculation processing (step S201) and the temperature calculation process (step S202), but the number of total number of time steps N t, is a process that is repeated.
  • the total number of time steps N t may be determined based on the simulated time range and time step size Delta] t.
  • Specifying the total number of time steps N t may be the in the calculation condition setting file is set the information of all time step number N t or simulated time range and time step size ⁇ t and is entered using the input device 20 By doing that.
  • the calculation process of radiant heat flux and the like performed in step S201 uses the parameters (form factors, etc.) calculated in the parameter calculation process, and the radiant flux G S, the short wave radiation (visible light) emitted from each element i , Calculate the radiation flux G L, i [W / m 2 ] of i [W / m 2 ] and long wave radiation (infrared rays), and use the calculated radiation flux to calculate the net of shortwave radiation that is absorbed by each element i This is a process for calculating the radiant heat R S, i [W] and the net radiant heat R L, i [W] related to long wave radiation.
  • n in the equations (8) and (9) is the total number of area elements and volume elements.
  • ⁇ A eff> i is the effective surface area of the volume element i when the element i is a volume element, and is the area of the area element i when the element i is an area element.
  • ⁇ S, i and ⁇ L, i are the reflectivities of the element i regarding short-wave radiation and long-wave radiation, respectively, and ⁇ i is the emissivity of the element i.
  • S direct, i is the direct shortwave radiation flux from the sun incident on element i
  • S diffuse, i is the radiation flux of atmospheric scattered shortwave radiation incident on element i.
  • L i is the radiation flux of atmospheric long-wave radiation incident on the element i
  • a eff i ⁇ Solar and A eff i ⁇ sky are effective areas of the element i in the solar direction and in the sky direction, respectively.
  • B (T i ) is a radiation flux emitted from the element i by thermal radiation.
  • B (T i ) is expressed by the following equation using the Stefan-Boltzmann constant ⁇ : Is calculated by
  • S direct, i and S diffuse, i are calculated by the following equations using the sky rate ⁇ i and the shade rate D i calculated in the parameter calculation process.
  • S ⁇ is a solar radiation flux incident downward on the horizontal plane
  • n i is a unit normal vector of the area element i
  • c direct and c diffuse are coefficients for direct diffusion separation.
  • S direct, i calculated by the above formula is the S direct, i of surface elements i.
  • S direct, i of the volume element i is calculated by the following equation using the effective shade rate D eff i calculated by the parameter calculation process.
  • injection rate ⁇ i is equal to the absorption rate of the element i
  • “1- ⁇ L , i ” can be used as the injection rate ⁇ i .
  • the radiant flux G S, i , GL, i emitted from each element i is calculated by solving these linear matrix equations.
  • the temperature calculation process (FIG. 2; step S202) is a process for calculating the temperature of each part in the simulation target space from the radiant heat calculated by the radiant heat flux calculation process.
  • the procedure for calculating the temperature of each area element during this temperature calculation process is general except that the radiant heat is calculated using a form factor that treats the canopy as a volume element having permeability.
  • the calculation procedure is the same. Therefore, only the procedure for calculating the temperature of the crown (volume element) will be described below.
  • a radiant heat flux R S of short wave radiation and a radiant heat flux R L of long wave radiation flow into a volume element that is a part of a tree crown, and the volume element From the sensible heat flux H and the latent heat flux LE. Therefore, the heat balance regarding the volume element i which is a tree crown is expressed by the following equation.
  • T leaf, i is the leaf surface temperature [K] in the element i
  • a i is the leaf area density [m 2 / m 3 ] in the element i
  • V i is the volume [m 3 ] of the element i
  • C is the heat capacity [J / K / m 2 ] of the leaves per unit leaf area.
  • R S, i and R L, i are the net radiant heat (intensity of radiant heat flux) [W] of short wave radiation and the net radiant heat [W] of long wave radiation absorbed by the leaves, respectively. Is the latent heat of vaporization [J / kg].
  • H i is the amount of sensible heat transported from the leaves to the atmosphere (stress heat flux intensity) [W]
  • E i is the amount of water vapor evacuated from the leaves to the atmosphere [kg / s].
  • the sensible heat transport amount H i released from the leaves to the atmosphere and the water vapor amount E i transpiration from the leaves to the atmosphere are calculated by the following equations.
  • T air i is the atmospheric temperature [K] in the element i
  • f air i is the water vapor partial pressure [Pa] in the atmosphere in the volume element i
  • f leaf i is in the volume element i. saturated water vapor partial pressure of the leaf surface [Pa] of
  • K h is the convective heat transfer coefficient [W / m 2 / K]
  • K w is the convective moisture transport coefficient [kg / s / m 2 / Pa]
  • ⁇ is Evaporation efficiency.
  • the leaf surface temperature T leaf, i at time step n + 1 after time step ⁇ t is calculated using the leaf surface temperature and heat flux at time step n.
  • the amount of change ⁇ T leaf, i in the leaf surface temperature from time step n to n + 1 is the net long wave radiation, sensible heat transport amount, transpiration amount due to the change in leaf surface temperature due to the passage of time ⁇ t.
  • the leaf surface temperature variation ⁇ T leaf, i is obtained by this equation, and then the leaf surface temperature T leaf, i (n + 1) at time step n + 1 is calculated by the following equation.
  • the temperature calculation process ends when the calculation of the temperature of each part and the output of the calculated temperature of each part (storage in the storage unit 14 in this embodiment) are completed. Then, when the designated number of processes has not been completed, the calculation process of the radiant heat flux and the like is started again, and the main process ends when the designated number of processes is completed.
  • the simulation apparatus 10 treats each tree crown as one or more volume elements having permeability, and the form factor ((5)) regarding one volume element and one area element. , (6) (see formula (6)), a form factor is calculated by decreasing the value by an amount corresponding to the amount of radiant heat transmitted through the one volume element. In addition, the simulation apparatus 10 calculates a view factor whose value is decreased by an amount corresponding to the amount of radiant heat transmitted through the two volume elements as the view factor (see Expression (7)) regarding the two volume elements. And the simulation apparatus 10 calculates the amount of radiant heat transferred between each element using each calculated form factor. Therefore, according to the simulation device 10, it is possible to satisfactorily simulate the radiant heat transport phenomenon in the three-dimensional space including the tree canopy without requiring a separate calculation of the state inside the trunk (that is, at a low calculation cost). .
  • the simulation apparatus 10 can perform various modifications.
  • the simulation device 10 can be transformed into a device that calculates a form factor in which the transmittance T between the two elements is not taken into account, and that considers the transmittance between the two elements when calculating the radiation flux.
  • the transmittance T between the two elements is taken into consideration when calculating the form factor, an accurate result can be obtained and the calculation cost can be reduced. Therefore, it is preferable to adopt the above-described form factor.
  • the simulation apparatus 10 can also be transformed into an apparatus that uses an analytical solution of the definition formula as a part or all of the form factors.
  • the simulation device 10 is a device that calculates the leaf temperature by the Euler implicit method.
  • the simulation device 10 is a device that calculates the leaf temperature by the explicit method or a device that calculates the leaf temperature by the second-order accuracy Crank-Nicholson method. It can also be transformed.
  • the implicit method is more likely to obtain an accurate value, it is preferable to employ the implicit method as the leaf temperature calculation method.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム
　本発明は、放射熱輸送現象をシミュレートするためのシミュレーション装置とシミュレーション方法とシミュレーションプログラムとに関する。
　都市部の気温がその周辺部の気温より高くなる現象を、ヒートアイランド現象という。　ヒートアイランド現象を緩和するための対策として、街路樹の配置や緑地の整備が着目されている。そのため、街路樹の配置等を考慮して都市部の放射熱輸送をシミュレートする様々なシミュレーション技術（特許文献１及び非特許文献１～９参照）が開発されている。
　ただし、既存のシミュレーション技術は、計算コスト（計算量）を抑えるために、『樹冠による放射熱の散乱を無視し、樹冠による放射熱の減衰・吸収のみを考慮する』といったように、樹冠に関する放射エネルギーの伝達過程を単純化したものとなっている。そのため、既存のシミュレーション技術では、樹冠付近や樹冠内部での放射エネルギーの吸収量等を適切にシミュレートできないことが多いのが現状である。
特開２００３－０９９６９７号公報 特許５１３７０３９号公報 特開２０１２－０２１６８４号公報
吉田 伸治, 村上 周三, 持田 灯, 大岡 龍三, 富永 禎秀, 「環境緩和効果を総合的に組み込んだ新しい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、各体積要素ｋの有効日陰率Ｄeff kがある。
　まず、パラメータ算出処理時に算出される各２要素に関する形態係数Ｆについて説明する。
　パラメータ算出処理時には、２つの面積要素ｉ，ｊの組み合わせ毎に、『面積要素ｉから面積要素ｊを見た形態係数Ｆij』と『面積要素ｊから面積要素ｉを見た形態係数Ｆji』とが、算出される。また、面積要素ｉと体積要素ｋの組み合わせ毎に、『面積要素ｉから体積要素ｋを見た形態係数Ｆik』と『体積要素ｋから面積要素ｉを見た形態係数Ｆki』とが、算出される。さらに、２つの体積要素ｋ及びｌの組み合わせ毎に、『体積要素ｋから体積要素ｌを見た形態係数Ｆkl』及び『体積要素ｌから体積要素ｋを見た形態係数Ｆlk』が算出される。
　パラメータ算出処理時に算出される『面積要素ｉから面積要素ｊを見た形態係数Ｆij』、『面積要素ｊから面積要素ｉを見た形態係数Ｆji』は、それぞれ、以下の（１）、（２）式で定義される値である。
Figure JPOXMLDOC01-appb-M000001
　これらの式において、Ａi、Ａjは、それぞれ、面積要素ｉの面積、面積要素ｊの面積である。βi、βj は、図３に模式的に示してあるように、それぞれ、面積要素ｉ内の微小面ｄＡiと面積要素ｊ内の微小面ｄＡjとを結ぶ直線と微小面ｄＡiの法線のなす角度、当該直線と微小面ｄＡjの法線のなす角度である。また、ｒは、微小面ｄＡi・微小面ｄＡj間の距離である。
　Ｔijは、面積要素ｉと面積要素ｊの間の透過率である。Ｔijは、２微小面ｄＡi及びｄＡj 間の光学的厚さτijを用いて次式により算出される。
Figure JPOXMLDOC01-appb-M000002
　また、面積要素ｉと面積要素ｊとの間（微小面ｄＡiと微小面ｄＡjとの間）に樹冠が分布している場合、光学的厚さτijは、樹冠の消散係数ｋext及び葉面積密度ａを用いて次式により算出される。
Figure JPOXMLDOC01-appb-M000003
　尚、形態係数の具体的な算出手順については後述するが、上記した（１）、（２）式で定義される形態係数は、相反則を満たしている。すなわち、面積Ａiと面積Ａjと形態係数Ｆijと形態係数Ｆjiとの間には、以下の関係がある。
Figure JPOXMLDOC01-appb-M000004
　従って、Ｆjiを、（１）式から算出したＦijとＡiとＡjとから算出することも、Ｆijを、（２）式から算出したＦjiとＡiとＡjとから算出することもできる。
　パラメータ算出処理時に算出される『面積要素ｉから体積要素ｋを見た形態係数Ｆik』は、以下の（３）式により定義される値である。
Figure JPOXMLDOC01-appb-M000005
　この（３）式におけるβiは、図４に模式的に示してあるように、面積要素ｉの微小面ｄＡiと体積要素ｊ内の微小投影面ｄＡk←iとを結ぶ直線と微小面ｄＡiの法線のなす角度である。また、ｒは、微小面ｄＡi・微小投影面ｄＡk←i間の距離である。
　Ａeff k←iは、体積要素ｋ自身の遮蔽率を考慮した，面積要素ｉの方向から見た体積要素ｋの有効面積である。このＡeff k←iは、次式により算出される。
Figure JPOXMLDOC01-appb-M000006
　ここで、Δτk←iは、微小投影面ｄＡk←i(図４参照)に垂直な方向の体積要素ｋの光学的厚さである。体積要素ｋが樹木の場合、Δτk←iは、消散係数ｋext、葉面積密度ａおよびｄＡk←iに垂直な方向の体積要素ｋの幾何学的厚さΔｓk←iを用いて、次式により算出される。
Figure JPOXMLDOC01-appb-M000007
　要するに、上記した（３）式は、（４）式を用いれば、以下のように表記できる式となっている。パラメータ処理時には、以下の（５）式に従って、面積要素ｉから体積要素ｋを見た形態係数Ｆikが算出される。
Figure JPOXMLDOC01-appb-M000008
　パラメータ算出処理時に算出される『体積要素ｋから面積要素ｉを見た形態係数Ｆki』は、次式により定義される値である。
Figure JPOXMLDOC01-appb-M000009
　換言すれば、体積要素ｋから面積要素ｉを見た形態係数Ｆkiとしては、以下の式で表される相反則を満たす値が算出される。
Figure JPOXMLDOC01-appb-M000010
　パラメータ算出処理時に算出される『体積要素ｋから体積要素ｌを見た形態係数Ｆkl』は、次式により定義される値である。
Figure JPOXMLDOC01-appb-M000011
　すなわち、この形態係数Ｆklは、体積要素ｋの遮蔽率（“1-exp(-Δτk←l)”）と体積要素ｌの遮蔽率（“1-exp(-Δτl←k)”）と体積要素ｋ、ｌ間の透過率Ｔklとを考慮に入れた形態係数となっている。
　尚、（７）式が表している体積要素に関する形態係数は、面積要素の形態係数と同様に次式の相反則を満たす。
Figure JPOXMLDOC01-appb-M000012
　従って、他の形態係数と同様に、２体積要素間の形態係数の双方を、（７）式に従って算出することも、一方の形態係数を（７）式に従って算出し、他方の形態係数を、当該一方の形態係数の算出結果から算出することも出来る。
　パラメータ算出処理時には、上記した各形態係数が、モンテカルロ法により算出される。
　すなわち、形態係数の算出時には、０から１までの範囲の一様乱数Ｒθ、Ｒφから、以下の式によりμ、φを算出し、算出結果に基づき、以下に示した内容の単位ベクトルｎを生成する処理が、必要とされる形態係数の精度に応じた回数、繰り返される。
Figure JPOXMLDOC01-appb-M000013
　換言すれば、形態係数の算出時には、Lambertの余弦則にしたがって面積要素（又は体積要素）から射出される多数の光子の進行方向単位ベクトルｎを生成する処理が行われる。
　そして、生成した各単位ベクトルｎの方向に要素ｉから等しいエネルギーＷ0を持った光子が射出された場合に要素ｊに入射されることになる各光子ｐのエネルギーＷpを積算することにより、形態係数Ｆijが求められる。
　具体的には、面積要素ｉから面積要素ｊをみた形態係数Ｆijは、面積要素ｉから射出させた光子の個数がＮ、面積要素ｊに入射された光子の個数がｎであった場合には、以下の式により算出される。
Figure JPOXMLDOC01-appb-M000014
　なお、建物壁面などの完全遮蔽性の面積要素のみの場合には、面積要素ｊに入射する光子のエネルギーＷpはＷ0と等しいため、形態係数は、Ｆij＝ｎ／Ｎにより算出できる。樹木のような放射透過性の要素が空間中に分布する場合には、面積要素ｊに入射する光子のエネルギーＷpは、面積要素ｉから面積要素ｊに到達するまでの間に減衰する。そのため、その減衰の形態係数に対する影響を光子の持つエネルギーを減衰させることで考慮する。すなわち、次式により、Ｗpが計算される。
Figure JPOXMLDOC01-appb-M000015
　ここで、Ｔij,pは、光子ｐの経路に沿った透過率である。なお、光子が透過性の体積要素内で確率的に遮蔽されると考えて、面積要素ｊに入射する光子の数を減少させても、減衰の影響を考慮した形態係数を求めることが出来る。
　また、面積要素ｉから体積要素ｋをみた形態係数Ｆikは、（７）式により定義される値である。従って、面積要素ｉから体積要素ｋをみた形態係数Ｆikは、次式により算出される。
Figure JPOXMLDOC01-appb-M000016
　ここで、Δτk←i,pは、光子ｐの経路に沿った、体積要素ｋの内部での光学的厚さである。Ｎは、面積要素ｉから射出された光子数であり、ｎは、体積要素ｋに入射された光子数である。
　体積要素ｋから面積要素ｉをみた形態係数Ｆkiの算出時には、体積要素の表面上の各点からの光子の等方的な射出が仮定される。すなわち、体積要素の表面上の各点からLambertの余弦則に従って多数の光子が仮想的に射出される。
　そして、体積要素ｋから面積要素ｉをみた形態係数Ｆkiが、次式により算出される。
Figure JPOXMLDOC01-appb-M000017
　ここで、Δτk←i,pは、光子ｐの射出点から光子ｐの進行方向とは逆の方向へバックトレースすることによって計算される，体積要素ｋの内部での光学的厚さである。また、Ｓkは、体積要素ｋの表面積であり、Ｍは、体積要素ｋから射出された光子数であり、ｍは、面積要素ｉに入射された光子数である。
　同様に、体積要素ｋから体積要素ｌをみた形態係数Ｆklは、次式により計算される。
Figure JPOXMLDOC01-appb-M000018
　尚、この式におけるｍは、Ｍ個の光子を体積要素ｋから射出した結果として体積要素ｌに入射された光子数である。
　以下、パラメータ算出処理で算出される、各体積要素ｋの有効表面積＜Aeff＞k、各要素ｉの天空率ωi、各面積要素ｉの日陰率Ｄi、各体積要素ｋの有効日陰率Ｄeff kについて説明する。
　パラメータ算出処理時に算出される各体積要素ｋの有効表面積＜Aeff＞kは、以下の式により定義される値である。
Figure JPOXMLDOC01-appb-M000019
 　この式において、ｍは、体積要素ｋの周囲に存在する（体積要素ｋから見える）要素（面積要素又は体積要素）の総数であり、ｉ（＝１～ｍ）は、体積要素ｋの周囲に存在する面積要素又は体積要素の要素番号である。
　パラメータ算出処理時には、この有効表面積＜Ａeff＞kが、次式により算出される。
Figure JPOXMLDOC01-appb-M000020
　すなわち、有効表面積＜Ａeff＞kもモンテカルロ法により算出される。
　要素（面積要素／体積要素）ｉの天空率ωiは、要素ｉから空をみた形態係数に相当する値である。天空率ωiは、要素ｉから面積要素をみた形態係数と同様の手順で算出される。
　面積要素ｉの日陰率Ｄiは、面積要素iから太陽側へ射出された光子ｐが、他の要素に入射した際に失うエネルギーΔＷｐを積算することにより求められる。より具体的には、日陰率Ｄiは、モンテカルロ法を用いて以下の式により算出される。
Figure JPOXMLDOC01-appb-M000021
　尚、上記式において、Ｎは、面積要素iから射出された光子数である。
　同様に、体積要素ｋの有効日陰率Ｄeff kは、モンテカルロ法を用いて以下の式により算出される。
Figure JPOXMLDOC01-appb-M000022
　上記した各面積要素ｉの日陰率Ｄi、各体積要素ｋの有効日陰率Ｄeff kは、太陽の位置により値が変わるパラメータである。従って、シミュレーションを行う時間範囲が比較的に長い場合には、パラメータ算出処理時に、当該時間範囲内の各時刻における日陰率Ｄi及び有効日陰率Ｄeff kが算出される。
　以下、演算部１２が行う主処理（図２）の内容を説明する。
　演算部１２が行う主処理は、放射熱フラックス等算出処理（ステップＳ２０１）と温度算出処理（ステップＳ２０２）とが、全時間ステップ数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] については、それぞれ、以下の式が成立する。
Figure JPOXMLDOC01-appb-M000023
　尚、（８）、（９）式におけるｎは、面積要素及び体積要素の総数である。<Ａeff>ｉは、要素ｉが体積要素である場合には、体積要素ｉの有効表面積であり、要素ｉが面積要素である場合には、面積要素ｉの面積である。
　αS,i、αL,iは、それぞれ、要素iの、短波放射、長波放射に関する反射率であり、εiは、要素ｉの射出率である。Ｓdirect,iは、要素ｉに入射する太陽からの直達短波放射フラックスであり、Ｓdiffuse,iは、要素iに入射する大気散乱短波放射の放射フラックスである。Ｌiは、要素ｉに入射する大気長波放射の放射フラックスであり、Ａeff i←Solar、Ａeff i←skyは、それぞれ、要素iの、太陽の方向、空の方向の有効面積である。
　Ｂ(Ｔi)は、要素ｉから熱放射により射出される放射フラックスである。長波放射の放射フラックスとして、上記したＧL,jのみを算出する場合（長波放射に関する放射フラックスを波長範囲別に算出しない場合）、Ｂ(Ｔi)は、ステファン・ボルツマン定数σを用いて次式により算出される。
Figure JPOXMLDOC01-appb-M000024
　Ｓdirect,i、Ｓdiffuse,iは、パラメータ算出処理で算出された天空率ωi、日陰率Ｄiを用いて以下の式により算出される。
Figure JPOXMLDOC01-appb-M000025
　ここで，Ｓ↓は水平面に下向きに入射する太陽放射フラックスであり、Ｓ（＝(Sx,Sy,Sz)）は太陽方位ベクトルである。ｎiは、面積要素ｉの単位法線ベクトルであり、cdirect、cdiffuseは、直散分離のための係数である。
　尚、上記式で算出されるＳdirect,iは、面積要素ｉのＳdirect,iである。体積要素ｉのＳdirect,iは，パラメータ算出処理で算出された有効日陰率Ｄeff iを用いて次式により算出される。
Figure JPOXMLDOC01-appb-M000026
　また、射出率εiは、要素ｉの吸収率と等しいため、射出率εiとして、“１－αL,ｉ”を使用することが出来る。
Figure JPOXMLDOC01-appb-M000027
　上記した（８）式及び（９）式は、ｉ＝１～ｎのそれぞれについて成立する。すなわち、以下の２つの線形行列方程式が成立する。
Figure JPOXMLDOC01-appb-M000028
　放射フラックス等算出処理時には、これらの線形行列方程式を解くことで、各要素iから射出される放射フラックスＧS,i，ＧL,iが、算出される。
　そして、各要素ｉに吸収される、短波放射に関する正味の放射熱ＲS,i[W]、及び長波放射に関する正味の放射熱ＲL,i[W]が、それぞれ、以下の（１０）、（１１）式により算出される。
Figure JPOXMLDOC01-appb-M000029
　温度算出処理（図２；ステップＳ２０２）は、放射熱フラックス等算出処理にて算出された放射熱からシミュレーション対象空間内の各部の温度を算出する処理である。この温度算出処理時における各面積要素の温度の算出手順は、放射熱が、樹冠を透過性を有する体積要素として取り扱った形態係数を用いて算出されたものであることを除けば、一般的な算出手順と同様のものである。そのため、以下では、樹冠（体積要素）の温度の算出手順のみを説明することにする。
　図５に模式的に示してあるように、樹木の樹冠の一部である体積要素には、短波放射の放射熱フラックスＲSと長波放射の放射熱フラックスＲLとが流入し、当該体積要素からは、顕熱フラックスＨと潜熱フラックスＬＥとが流出する。
　従って、樹冠である体積要素ｉに関する熱収支は、次式により表されることになる。
Figure JPOXMLDOC01-appb-M000030
　ここで，Ｔ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は、以下の式により算出される。
Figure JPOXMLDOC01-appb-M000031
　ここで、Ｔair,iは、要素ｉ内の大気温度[K]，fair,iは、体積要素ｉ内の大気中の水蒸気分圧[Pa]、ｆleaf,iは、体積要素ｉ内の葉表面の飽和水蒸気分圧[Pa]，Ｋhは、対流熱伝達係数[W/m2/K]，Ｋwは、対流水蒸気輸送係数[kg/s/m2/Pa]、βは蒸発効率である。
　温度算出処理時には、時間ステップｎにおける葉表面温度及び熱フラックスを用いて、時間刻みΔt後の時間ステップｎ＋１における葉表面温度Ｔleaf,iが、計算される。
　具体的には、時間ステップｎからｎ＋１までの間の葉表面温度の変化量ΔＴleaf,iは，時間Δtの経過に起因する葉表面温度の変化による正味長波放射，顕熱輸送量，蒸散量の変化を考慮すると、次式により表される。
Figure JPOXMLDOC01-appb-M000032
　従って，葉表面温度の変化量ΔＴleaf,iを、次式により求めることができる。
Figure JPOXMLDOC01-appb-M000033
　温度算出処理時には、この式により、葉表面温度の変化量ΔＴleaf,iが求められてから、時間ステップｎ+１における葉表面温度Ｔleaf,i (n+1)が次式により、算出される。
Figure JPOXMLDOC01-appb-M000034
　温度算出処理は、各部の温度の算出と、算出した各部の温度の出力（本実施形態では、記憶部１４への記憶）とが完了したときに終了する。そして、指定された回数の処理が完了していなかった場合には、再び、放射熱フラックス等算出処理が開始され、指定された回数の処理が完了したときに、主処理が終了する。
　以上、説明したように、本実施形態に係るシミュレーション装置１０は、各樹冠を１つ以上の透過性を有する体積要素として取り扱って、１つの体積要素と１つの面積要素に関する形態係数（（５）、（６）式参照）として、当該１つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。また、シミュレーション装置１０は、２つの体積要素に関する形態係数（（７）式参照）として、当該２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する。そして、シミュレーション装置１０は、算出した各形態係数を用いて各要素間で授受される放射熱量を算出する。従って、シミュレーション装置１０によれば、樹幹内部の状態を別途計算する必要がない形で（つまり、低計算コストで）、樹冠を含む３次元空間における放射熱輸送現象を良好にシミュレートできることになる。
　《変形形態》
　上記した実施形態に係るシミュレーション装置１０は、各種の変形を行えるものである。
　例えば、シミュレーション装置１０を、２要素間の透過率Ｔが考慮されていない形態係数を算出し、放射フラックスの算出時等に２要素間の透過率を考慮する装置に変形することが出来る。ただし、形態係数の算出時に２要素間の透過率Ｔを考慮しておいた方が正確な結果が得られるし、計算コストも低くなる。従って、上記した形態係数を採用しておくことが好ましい。
　また、樹幹の体積熱容量Caは、通常、かなり小さいので、シミュレーション装置１０を、葉表面温度を、Ca=0として算出する装置に変形することも出来る。また、シミュレーション装置１０を、一部又は全ての形態係数として、その定義式の解析解を用いる装置に変形することも出来る。
　シミュレーション装置１０は、オイラー陰解法により葉温度を算出する装置であるが、シミュレーション装置１０を、陽解法により葉温度を算出する装置や、２次精度のクランク・ニコルソン法により葉温度を算出する装置に変形することも出来る。ただし、陰解法の方が正確な値を得やすいので、葉温度の算出法には、陰解法を採用しておくことが好ましい。
１０ シミュレーション装置
１２ 演算部
１４ 記憶部
１６ インタフェ－ス（Ｉ／Ｆ）回路

Claims (6)

  1. 　放射熱輸送現象をシミュレートするシミュレーション装置において、
    　複数の面積要素と複数の体積要素とで定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出手段であって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出手段と、
    　前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出手段と、
    　を備え、
    　前記３次元空間内に存在する複数の樹木の樹冠が前記複数の体積要素として取り扱われるように、前記３次元空間が前記複数の面積要素と前記複数の体積要素とにより定義されている
    　ことを特徴とするシミュレーション装置。
  2. 　前記形態係数算出手段は、各２要素に関する形態係数として、各２要素間におけう放射熱量の透過率を考慮に入れた形態係数を算出する
    　ことを特徴とする請求項１に記載のシミュレーション装置。
  3. 　前記放射熱量算出手段により算出された放射熱量に基づき、各要素の温度を算出する温度算出手段を、さらに備える
    　ことを特徴とする請求項１又は２に記載のシミュレーション装置。
  4. 　前記温度算出手段は、各要素の温度を陰解法により算出する
    　ことを特徴とする請求項３に記載のシミュレーション装置。
  5. 　放射熱輸送現象をシミュレートするシミュレーション方法において、
    　コンピュータが、
    　複数の樹木の樹冠が複数の体積要素として取り扱われるように、複数の面積要素と複数の体積要素とが定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出ステップと、
    　前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
    　を行うことを特徴とするシミュレーション方法。
  6. 　放射熱輸送現象をシミュレートするためのシミュレーションプログラムにおいて、
    　コンピュータに、
    　複数の樹木の樹冠が複数の体積要素として取り扱われるように、複数の面積要素と複数の体積要素とが定義された仮想的な３次元空間内の各２要素に関する形態係数を算出する形態係数算出ステップであって、１つ又は２つの体積要素を含む２要素に関する形態係数として、当該１つ又は２つの体積要素を透過する放射熱量に応じた分だけ値を減少させた形態係数を算出する形態係数算出ステップと、
    　前記形態係数算出手段により算出された各形態係数を用いて各要素間で授受される放射熱量を算出する放射熱量算出ステップと、
    　を行わせることを特徴とするシミュレーションプログラム。
PCT/JP2014/080198 2013-11-14 2014-11-14 放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム WO2015072546A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
US15/036,820 US10796034B2 (en) 2013-11-14 2014-11-14 Simulation apparatus, simulation method, and simulation program relating to radiation heat transport phenomenon
DE112014005224.8T DE112014005224T5 (de) 2013-11-14 2014-11-14 Simulationsvorrichtung, Simulationsverfahren und Simulationsprogramm für das Strahlungswärme-Transportphänomen
CN201480062470.9A CN105765585B (zh) 2013-11-14 2014-11-14 关于辐射热输送现象的模拟装置和模拟方法

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013-236190 2013-11-14
JP2013236190A JP6238439B2 (ja) 2013-11-14 2013-11-14 放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム

Publications (1)

Publication Number Publication Date
WO2015072546A1 true WO2015072546A1 (ja) 2015-05-21

Family

ID=53057480

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/080198 WO2015072546A1 (ja) 2013-11-14 2014-11-14 放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム

Country Status (5)

Country Link
US (1) US10796034B2 (ja)
JP (1) JP6238439B2 (ja)
CN (1) CN105765585B (ja)
DE (1) DE112014005224T5 (ja)
WO (1) WO2015072546A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114896854A (zh) * 2022-05-12 2022-08-12 大连理工大学 一种用于托卡马克中等离子体辐射演化的模拟方法
CN117455977A (zh) * 2023-09-27 2024-01-26 杭州市交通工程集团有限公司 一种基于三维激光扫描的堆料体积计算方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009057189A1 (ja) * 2007-10-29 2009-05-07 Japan Agency For Marine-Earth Science And Technology 熱放射エネルギーのシミュレーション装置、及び、方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3686931B2 (ja) 2001-09-21 2005-08-24 財団法人理工学振興会 熱環境の予測方法、およびプログラム
JP3964799B2 (ja) * 2003-02-07 2007-08-22 東京瓦斯株式会社 放射伝熱解析方法及び装置
US8781801B2 (en) * 2007-10-29 2014-07-15 Japan Agency For Marine-Earth Science And Technology Meteorological phenomena simulation device and method
CN101320098B (zh) * 2008-07-11 2012-04-18 重庆大学 基于数字图像分析的城镇热岛特性预测方法及系统
JP5648784B2 (ja) 2010-07-13 2015-01-07 国立大学法人横浜国立大学 加熱装置

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2009057189A1 (ja) * 2007-10-29 2009-05-07 Japan Agency For Marine-Earth Science And Technology 熱放射エネルギーのシミュレーション装置、及び、方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
AYA HAGISHIMA: "Numerical Analysis on Cooling Effect of a Row of Trees in an Urban Canyon", JOURNAL OF ARCHITECTURE, PLANNING AND ENVIRONMENTAL ENGINEERING, 30 November 1999 (1999-11-30), pages 83 - 90 *
SHINJI YOSHIDA: "Development of New Plant Canopy Model for Coupled Simulation of Convection, Radiation and Moisture Transport", MONTHLY JOURNAL OF INSTITUTE OF INDUSTRIAL SCIENCE, vol. 51, no. 1, January 1999 (1999-01-01), pages 11 - 16 *

Also Published As

Publication number Publication date
DE112014005224T5 (de) 2016-08-18
US20160267206A1 (en) 2016-09-15
CN105765585A (zh) 2016-07-13
JP6238439B2 (ja) 2017-11-29
JP2015095239A (ja) 2015-05-18
US10796034B2 (en) 2020-10-06
CN105765585B (zh) 2019-10-22

Similar Documents

Publication Publication Date Title
Mirhosseini et al. View factor calculation using the Monte Carlo method for a 3D strip element to circular cylinder
Peterson et al. Simulation of astronomical images from optical survey telescopes using a comprehensive photon Monte Carlo approach
WO2015072545A1 (ja) シミュレーション装置、シミュレーション方法及びシミュレーションプログラム
Jensen et al. On various modeling approaches to radiative heat transfer in pool fires
Sun et al. An implicit unified gas kinetic scheme for radiative transfer with equilibrium and non-equilibrium diffusive limits
Zhang et al. The influence of different black carbon and sulfate mixing methods on their optical and radiative properties
WO2015072546A1 (ja) 放射熱輸送現象に関するシミュレーション装置、シミュレーション方法及びシミュレーションプログラム
Battaglia et al. Forward Monte Carlo computations of fully polarized microwave radiation in non-isotropic media
Wang et al. Discontinuous finite element method for vector radiative transfer
Trovalet et al. Modified finite-volume method based on a cell vertex scheme for the solution of radiative transfer problems in complex 3D geometries
US8306792B2 (en) Simulator and simulating method of heat radiation energy
JP2016149976A (ja) 熱環境シミュレーション方法及び装置並びにプログラム
Li et al. The point spread function reconstruction by using Moffatlets—I
CN112347708A (zh) 一种基于物理过程的实时光学烟幕的渲染方法
Viskanta Computation of radiative transfer in combustion systems
Sánchez et al. A three-dimensional atmospheric radiative transfer model based on the discrete-ordinates method
Fu et al. Sky polarization pattern under multi-layer environment of atmosphere and sea fog
Shang et al. A computational approach for hypersonic nonequilibrium radiation utilizing space partition algorithm and Gauss quadrature
Ignat’Ev et al. A new version of the discrete ordinate method for the calculation of the intrinsic radiation in horizontally homogeneous atmospheres
Hilton et al. Dynamic modelling of radiant heat from wildfires
CN104819774A (zh) 一种基于微透镜阵列的火焰光场探测泛尺度分析方法
Zhang et al. Monte Carlo method of radiative transfer applied to a turbulent flame modeling with LES
Sun et al. Radiance and polarization in the diffusion region with an arbitrary scattering phase matrix
Wu et al. Three-dimensional mesoscopic simulation of radiative transfer in graded index media
Nair A high-order multiscale global atmospheric model

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 14862406

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 15036820

Country of ref document: US

WWE Wipo information: entry into national phase

Ref document number: 1120140052248

Country of ref document: DE

Ref document number: 112014005224

Country of ref document: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14862406

Country of ref document: EP

Kind code of ref document: A1