図面を参照して、本発明の実施形態につい� ��説明する。実施形態の説明は、例示であり� ��本発明は、実施形態の構成に限定されない� ��
〔実施形態〕
 〈システム構成〉
 図1は、本実施形態のシステム構成を示す図 である。本実施形態のシステムは、シミュレ ーション装置100、入力装置200、出力装置300、 記憶装置400を含む。

 シミュレーション装置100は、主に計算(演 算)を行う演算部102、初期データやプログラ� �等を記憶する記憶部104、外部装置との入出� �を行うインタフェース(I/F)106、演算部102で実 行されるプログラム及びデータが展開される メモリ108、外部との通信機能を有する通信部 110を含む。

 演算部102は、地形データ取得手段、気象� ��ータ取得手段、3次元データ取得手段、気象 データ算出手段、グリッド属性データ取得手 段、熱放射エネルギー算出手段、境界算出手 段、境界気象データ算出手段、出力手段とし て機能し得る。また、これらの手段のうちの 1つ、または、複数の手段が、演算部102とは� �の独立の演算部で機能してもよい。これら� �各構成要素は、シミュレーション装置100の� �算部102上で実行されるプログラムとして構� �することもできる。また、各構成要素の一� �、または、すべてをハードウェア回路で構� �してもよい。

 入力装置200は、キーボード、マウス等の� ��インティングデバイス等を含む。出力装置3 00は、LCDやCRT等のディスプレイ、プリンタ等� ��含むコンピュータである。記憶装置400は、� ��ードディスクドライブ(HDD)、フレキシブル� �ィスクドライブ、CDドライブ、DVDドライブ等 を含む。記憶装置400は、記憶部104に格納され るデータ等を、記憶装置400に装着する記録媒 体に記録することもできる。記憶装置400は、 記憶部104の代替の手段として、使用すること ができる。

 〈処理フロー〉
 図2は、本実施形態の気象シミュレーション の処理フローを示す図である。

 最初に、シミュレーションの対象となる� ��域の地形データ、気象データ、都市空間の3 次元形状のデータを入力する(S1002)。

 地形データには、シミュレーションの対� ��となる領域内の、緯度、経度、標高、地質� ��が含まれ得る。ここでいう地質は、地表面� ��属性(裸地、芝生、舗装面等の地面の特性)� �表す。地質は、地表面の放射率に影響を与� �る。放射率は、地表面の属性の1つである。� ��質の相違によって、地表面の放射率が相違� ��る。気象データには、シミュレーションの� ��象となる領域内の全部または一部の地点の� ��風(風向、風速等)、気温、気圧、湿度、水� �気量等が含まれ得る。都市空間の3次元形状� ��、人工の構造物(ビルの壁面、屋上面、道路 の路面など)、自然の構造物(木、崖、地表面� ��ど。雲などを含むこともできる)などを含む 立体物の現実に即した3次元座標を指定する� �とによって、定義される。3次元座標は、直� ��座標系、極座標系等、どのような座標系に� ��って与えられてもよい。入力される都市空� ��の3次元形状は、シミュレーションの対象と なる領域の一部でもよい。

 地形データ、気象データ、3次元形状デー タは、ユーザが入力装置200を使用してデータ を入力すること等により、定義することが可 能である。入力されたデータは、記憶部104に 格納される。

 地形データ、気象データ、3次元形状デー タは、入力装置から入力されたものが記憶部 104やメモリ108に格納されるようにしてもよく 、記憶装置400にあらかじめ格納されたものが 記憶部104、メモリ108に格納されるようにして もよく、通信部110で受信したものが記憶部104 や、メモリ108に格納されるようにしてもよい 。

 演算部102は、シミュレーションに必要な� ��地形データ、気象データ、立体物の3次元形 状のデータ、定数、変数等を、メモリに取得 し、シミュレーションに使用する。取得とは 、入力装置200から入力されるものをメモリに 格納すること、記憶部104、記憶装置400にあら かじめ格納されているものをメモリに読み込 むこと、通信装置で受信されるものをメモリ に読み込むこと等を含み得る。

 次に、形態係数、天空率、日陰率を算出� ��る(S1004)。

 シミュレーション装置100の演算部102は、S 1002で定義された3次元形状を、所定の数のグ� ��ッド(面要素)に分割する。また、所定の面� �を有するグリッドに分割してもよい。3次元� ��状のグリッドのデータは、記憶部104に格納� ��れる。分割数、グリッドの面積等を、ユー� ��が指定するようにすることもできる。分割� ��、グリッドの面積等を制御することにより� ��計算結果の精度を制御することができる。� ��た、各グリッドの面積が3次元形状から求め られ、記憶部104に格納される。

 また、演算部102は、あらかじめ分割され� ��状態の立体物の3次元形状のデータを、取得 してもよい。取得とは、入力装置200から入力 されるものをメモリに格納すること、記憶部 104、記憶装置400にあらかじめ格納されている ものをメモリに読み込むこと、通信装置で受 信されるものをメモリに読み込むこと等を含 み得る。

 演算部102は、複数のグリッドに分割され� ��3次元形状のグリッド毎に、すべてのグリッ ドとの間の形態係数を計算する。また、演算 部102は、複数のグリッドに分割された3次元� �状のグリッド毎に、天空率、日陰率を計算� �る。形態係数及び形態係数の計算、天空率� �日陰率の計算については、後に詳述する。� �リッド毎に算出された形態係数は、記憶部10 4に格納される。

 シミュレーション装置100は、他の装置等� ��よって算出された形態係数、天空率、日陰� ��を、取得してもよい。

 次に、大気中の気象データを算出する(S10 06)。演算部102は、地形データ、気象データを 使用して、気象モデルに基づいて、シミュレ ーションの対象となる領域の大気中の風、気 温、湿度、気圧等の気象データを算出する。 気象モデルでは、シミュレーションの対象と なる領域を、任意の大きさのグリッドに分割 することができる。気象モデルで使用するグ リッドは、都市空間の3次元形状が与えられ� �いる領域では、都市空間の3次元形状のグリ� ��ドと同じスケールとすることができる。ま� ��、それ以外の領域では、都市空間の3次元形 状のグリッドよりも大きなグリッドを使用す ることができる。気象モデルについては、後 に詳述する。

 次に、熱放射エネルギーの計算をする(S10 08)。演算部102は、都市空間の熱放射エネルギ ーを計算するため、行列要素のデータ作成を する。各グリッドの温度は、S1006で算出され� ��、温度が境界条件として使用される。各グ� ��ッドの放射率は、ユーザが入力装置200を使� ��して入力する等により、与えられる。各グ� ��ッドの面積、形態係数は、S1004で求められ� �いる。各グリッドの放射率、温度、面積、� �態係数等の各データは、記憶部104に格納さ� �る。演算部102は、これら記憶部104に格納さ� �た各データから、熱放射エネルギー計算の� �めの行列要素データを作成する。

 演算部102は、作成した行列要素データを� ��に、各グリッドの熱放射エネルギーの計算� ��する。熱放射エネルギーの計算については� ��後に詳述する。演算部102は、グリッド毎の� ��放射エネルギーの計算結果を、所定の出力� ��置300、記憶装置400、記憶部104などに出力す� ��。

 次に、都市空間のグリッドのフラックスの� ��算をする(S1010)。演算部102は、S1008で算出し� ��、熱放射エネルギーG _i を使用して、壁面表面へのフラックス、室内 壁面へのフラックス、地表面下向きフラック ス、地表面上向きフラックス、壁面対流熱伝 達フラックスを計算する。また、演算部102は 、都市空間の圧力変動を計算する(S1012)。圧� �変動は、気温の関数である。よって、圧力� �動の代わりに、温度変化として表すことも� �きる。各フラックス、圧力変動については� �後に詳述する。

 次に、演算部102は、都市空間の気象デー� ��の更新を行う(S1014)。演算部102は、算出した 各フラックス、圧力変動を使用して、所定の 微小時間ステップ後の都市空間の建物及び地 表面等と大気との境界の気温を含む気象デー タの算出を行う。これにより、所定の微小時 間ステップ後の都市空間の3次元形状等を反� �した気温を含む気象データを得ることがで� �る。演算部102は、気温とは異なる気温と等� �の気象データを算出してもよい。

 その後、演算部102は、計算を終了するか� ��かの判断(S1016)をする。演算部102は、計算を 終了しない場合(S1016;NO)、時刻の更新を行い� �所定の微小時間ステップ後の都市空間の影� �を反映した気温を含む気象データを使用し� �、当該所定の微小時間ステップ後の大気中� �気象データの計算を行う(S1006)。

 演算部102により算出された気象データ等� ��、適宜、出力装置300から出力することがで� ��る。また、演算部102により算出された気象� ��ータ等は、適宜、記憶部104、記憶装置400に� ��納することができ、通信部110を使用して、� ��部の装置に送信することもできる。

 データの作成、計算は、演算部102が記憶� ��104に記憶されたプログラムをメモリ108に展� ��して実行することにより行われる。このと� ��、入力装置200から入力されたデータ、記憶� ��104、記憶装置400にあらかじめ格納されたデ� ��タ、あるいは、通信部110で受信されるデー� ��が適宜使用される。

 (形態係数)
 図3は、形態係数を説明する図である。

 図3のように、面積がA _i 、A _j の面i、jがあるとする。面iから全方位に放射 される熱放射エネルギーは、E _i A _i である。ここで、面iの放射照度をEiとする。 面iから面jに直接入射するエネルギーをQ _ij とする。面iから全方位に放射されるエネル� �ーのうち、直接面jに入射するエネルギーの� ��合をF _ij とすると、

となる。このF _ij を形態係数という。

 図3のように、面iと面jとの距離はrであると する。ここで、面i、面j上に、微小な有限面� ��素dA _i 、dA _j をとる。面iと面jとを結ぶ直線と、面i、面j� �法線とのなす角を、それぞれβ _i 、β _j とする。有限面要素dA _i からdA _j をのぞむ立体角dω _j 内に単位時間当たり射出される放射エネルギ ーdQ _ij は、面iの放射輝度I _bi を用いて、

で求められる。ここで立体角の定義より、

であるから、放射エネルギーdQ _ij は、

となる。また、逆にdAjからdAiに対して単位時 間当たり射出される放射エネルギーdQ _ji は、

となる。dA _i とdA _j との間の正味の放射伝熱量dQ _net,ij は、

を用いて、面i、面jの温度が、それぞれT _i 、T _j の場合、

となる。面i、面jとの間の正味の放射伝熱� ��は、上式を積分することにより得られる。

ここで、F _ij 、F _ji は、面iと面jとの間の形態係数であり、

で与えられる。面iと面jとの間の形態係数� ��、面iと面jとの間の幾何学的な相対位置関� �で決まる。

 (形態係数の計算)
 上記の形態係数を求める際、モンテカルロ� ��を採用する。すなわち、壁面等の各グリッ� ��から、光子をランダム方向へ飛ばして光子� ��挙動を模擬する。光子を射出したグリッド� ��、個々の光子を追跡する過程で遮蔽される� ��リッドとの関係を検出する。

 建物や道路などの壁面上のあるグリッド� ��着目する。以下、任意のグリッドi(i=1~n)が� �表する壁面のことを、グリッド面i(i=1~n)とす る。

 モンテカルロ法では計算精度を向上させる� ��は膨大な粒子(光子)数が必要である。例え� �、各グリッドから10 ⁶ 光子を飛ばした場合、解析解との誤差は10 ^-4 オーダーである。形態係数の必要な精度に応 じて、飛ばす光子の数を増減することができ る。モンテカルロ法を用いることにより、建 物等の複雑な3次元形状であっても、容易に� �精度の形態係数を算出することができる。

 形態係数は、3次元形状に基づいて解析解 によって求められてもよい。

 図4は、放射方向の単位ベクトルを示す図 である。方位角θ、天頂角η方向にビームを� �射する場合、放射方向の単位ベクトルは次� �ように表される。

 このとき放出点を中心にした単位球面上� ��面積要素は、

と表される。

 Lambertの余弦則が成り立つ場合には、θに� ��いては一様であるが、ηについてはその方� �への放射確率がcosηに比例する。

 したがって、積分要素は、

となる。θ、ηに対しては、一様乱数R _θ、 R _η を用いてそれぞれ、

とすればよい。cは定数であり、両辺を積� �することにより、

 が導かれ、

から、

が導かれる。

 以上より、Lambertの余弦則が成り立つ場合 に、グリッドから光子を放射する頻度は、以 下のように書ける。

 グリッド毎に、この光子の放射頻度に従� ��て、光子を所定数ずつ放射する。光子の1つ 1つが、所定のエネルギーに相当する。よっ� �、各グリッドに直接入射した光子数を計数� �ることにより、形態係数を求めることがで� �る。放射する光子数を増やすことで、精度� �高い形態係数を求めることができる。

 図5は、グリッド面Aから放射した光子の様� �を示す図である。グリッド面Aから上記の放� ��頻度に従って放射した光子の一部は、グリ� ��ド面Bに直接到達している。よって、グリッ ド面Aから飛ばした全光子の数をNとし、グリ� ��ド面Aからグリッド面Bに直接到達した光子� �数をnとすると、形態係数F _AB は、

となる。

 (天空率の計算)
 天空率は、グリッド面から空が見える割合� ��いう。つまり、天空率は、グリッド面から� ��射される光子が他の壁面や地面に遮蔽され� ��に大気中に抜ける光子の割合である。

 従って、各グリッド面の天空率f _sky は、次のように求められる。

 (日陰率の計算)
 日射率は、計算対象のグリッドが、計算対� ��の時から所定時間前までの時間に対する当� ��所定時間の間に日陰であった時間の割合で� ��る。日陰率の計算は、形態係数の計算と同� ��に、各グリッド面から光子を所定時間の間� ��太陽の位置が想定される方向に放射し、建� ��等に遮蔽されるか否かを調べることにより� ��計算できる。形態係数の計算とは異なり、� ��陰率の計算の際、各グリッド面から放射す� ��光子は、所定時間の間に想定される太陽の� ��置の方向のみである。

 日陰率は、次の式で表される。

 (気象モデル)
 ここでは、気象モデルの例として、非静力� ��・大気大循環モデル(MSSG-A、Atmospheric componen t of the Multi-Scale Simulator for the Geoenvironment )について説明する。気象モデルとしては、� �れに限らず、どのような気象モデルを使用� �てもよい。

 非静力学・大気大循環モデルは、メソス� ��ールの気象現象を対象とする全球および領� ��に適用可能な気象モデルである。

 MSSG-Aは、主に次の特徴を有する。
(1)3次元完全圧縮性のナヴィエ・ストークス� �程式を使用する。
(2)全球はYin-Yang格子系を採用し、領域との1way -結合あるいは2way-結合が可能である。
(3)離散化手法は、5次風上差分法、3次のルン� ��・クッタ法、HE-VI(Horizontally Explicit - Vertica lly Implicit)型スプリッティング法を採用する� ��
(4)乱流モデルとして、LES(Large Eddy Simulation)� �デルを採用し、乱流スケールの詳細流れが� �現可能である。
(5)複数の雲微物理モデルや放射モデルが使用 可能であり、陸面モデル、海氷モデルも導入 され得る。
(6)大気流れの基本方程式として、完全圧縮の 非静力学方程式、連続の式(質量保存式)、状� ��方程式を使用する。

 非静力学・大気大循環モデルでは、初期� ��件として地形データを与え、パラメータと� ��て、ある時刻の、風、気温、湿度、圧力等� ��気象データを入力すると、地形データに基� ��いた、当該時刻の計算領域全体の風、気温� ��湿度、圧力がグリッド毎に算出することが� ��きる。非静力学・大気大循環モデルでは、� ��れらの気象データに基づいて、さらに、大� ��中の雲、降雨、降雪、雹などの予測を算出� ��ることができる。気象データ算出の際、す� ��てのグリットの気象データを、入力データ� ��して与える必要はない。非静力学・大気大� ��環モデルでは、データのないグリッドにつ� ��ては、適切にデータが補間されるからであ� ��。

 また、非静力学・大気大循環モデルでは� ��グリッド毎に出力された気象データに基づ� ��て、任意の時間経過後の計算領域全体のグ� ��ッド毎の気象データの算出を行うことがで� ��る。

 気象データを入出力する単位であるグリ� ��ドの大きさは、領域により制御することが� ��きる。よって、都市空間ではグリッドの密� ��を高くしてより詳細な気象データを与え、� ��市空間以外ではグリッドの密度を小さくし� ��気象データを与えることができる。

 MSSG-Aの基礎方程式系には、非静力学の条� ��のもとで、完全圧縮方程式系を採用してい� ��。この方程式系を次に示す。

 風のλ成分の方程式(予測式)は、次式で表 される。

 風のφ成分の方程式(予測式)は、次式で表 される。

 風のz成分の方程式(予測式)は、次式で表� ��れる。

 圧力変化の方程式(予測式)は、次式で表� �れる。

 また、連続の式は、次式で表される。

 気体の状態方程式は、次式で表される。

 ただし、ある変数xに対してA(x)は次式で� �義される。

である。

 ここで、

であり、

である。

 基本場とは、水平一様でかつ静力学平衡� ��態にある大気場である。密度の基本場とは� ��水平一様でかつ静力学平衡状態にある大気� ��における密度の鉛直分布である。また、気� ��の基本場とは水平一様でかつ静力学平衡状� ��にある大気場における気圧の鉛直分布であ� ��。外力項とは、計算対象とするシステムの� ��の影響によって計算対象とするシステム内� ��与えられる力を現す項である。加熱項とは� ��計算対象とするシステムの外の影響によっ� ��計算対象とするシステム内に与えられる熱� ��を表す項である。

 計算される大気中の気象データはρu、ρv� ��ρw、ρ’、p’であり、これらの変数の時間� ��化を計算すれば大気の流れを予測すること� ��できる。

 計算初期条件の気象データ、及び、都市� ��間における数値計算で求めた都市空間の気� ��データ等を、上記の式に代入することで、� ��期時間から任意時間経過後の都市空間を含� ��計算領域全体の都市空間の影響を含む気象� ��ータを得ることができる。

 (熱放射エネルギー)
 任意のグリッド面である非黒体面における� ��射伝熱過程は、他の任意のグリッド面から� ��熱放射線の反射を考慮する必要がある。単� ��面積から単位時間内に、空間に放出される� ��放射エネルギーG _i は、次の式で書ける。

 右辺第1項は、自己より放射される熱放射 エネルギーである。右辺第2項は、任意のグ� �ッド面からの熱放射線の反射による熱放射� �ネルギーである。

 ここで、任意のグリッド面i(i=1~n)につい� �、

かつ、

が成立する。放射率(ε _i )は、地表面の地質、壁面等の材質などによ� �て決まる。放射率は、黒体であれば1である� ��

 グリッド面iの温度(T _i )、放射率(ε _i )は、あらかじめ与えられる。グリッド面iの� ��積(A _i )は、建物等の3次元形状から求めることがで� ��る。また、形態係数(F _ij )は、上述したように求めることができる。

 よって、n個の熱放射エネルギーG _i の条件式についての次の連立方程式を解くと 、各グリッド面i(i=1~n)の熱放射エネルギーG _i が求まる。この式の行列要素データは、上記 のパラメータ(放射率、温度等)から作成する� ��とができる。

 また、面iから放出される正味の放射熱フラ ックスq _Ri は面iから放出される正味の放射伝熱量Q _Ri を用いて次の式で求めることができる。

 (壁面表面へのフラックス)
 壁面のグリッド面の番号をi _s 、x方向の格子番号をi _x 、y方向の格子番号をj _y とすると、壁面のグリッド面の表面へのフラ ックスは、次の式(1)で得られる。

 ここで、

である。

 式(1)の右辺の第1項は、太陽からの短波長 の直達放射フラックスである。第2項は、太� �の短波長の全天からの散乱による放射フラ� �クスである。第3項は、壁面からの長波長放� ��フラックス(熱放射エネルギー:射度)である� ��第4項は、周辺の壁面からこの壁面に入射す る長波長フラックスである。第5項は、大気� �ら壁面に入射する長波入射フラックスであ� �。第6項は、壁面表面から大気への対流熱伝� ��によるフラックスである。

 太陽から地表面へ入射する短波入射フラ� ��クス、および、大気から壁面に入射する長� ��入射フラックスは、他の手段に求められた� ��が、あらかじめ与えられる。

 (室内壁面のフラックス)
 室内壁面から室内への熱フラックスは、次� ��式(2)で与えられる。

 ここで、

である。

 (地表面下向きフラックス)
 地表面下向きフラックスQ _down は、大気から地表面への下向き長波フラック スq _long と、周辺壁面から地表面に入射する長波フラ ックスδq _long の寄与を加算する必要がある。周辺壁面から の寄与は、次の計算で得られる。

 ここで、式(4)の第2の等号は、次の射度の 式から導かれる。

となる。

 (地表面上向きフラックス)
 地表面上向きフラックスQ _up は、地表面から大気への上向き長波フラック スに加えて、周辺壁面から地表面に入射する 長波フラックスが地表面で反射される寄与を 加算して、熱放射エネルギーを計算する必要 がある。地表面の射度は、両者の寄与を加え たものであるので、熱放射エネルギーQ _up は、次のようになる。

 (壁面対流熱伝達フラックス)
 壁面対流熱伝達フラックスは、次の式で与� ��られる。

 ここで、

である。

 (圧力変動)
 圧力変動は、大気のエネルギー変化を圧力� ��化に換算した量である。壁面の温度変化に� ��って変化する単位時間当たりの圧力変動は� ��次の通りである。

 ここで、

である。

 (外壁と内壁との間の熱伝導)
 図6は、外壁と内壁との間の熱伝導を説明す る図である。左側が外壁側、右側が内壁側を 表す。各壁面の任意のグリッドにおける温度 、および、壁面から放出される圧力変動を求 めるために、壁面を多層のスラブでモデル化 した。

 ここで、次の変数を使用する。

 次の式が、熱伝導モデルの基本方程式で� ��る。

 建物は、屋根、壁から構成される中空な� ��体とし、建物の屋根、方位ごとの壁につい� ��熱容量を考慮した1次元熱伝導の多層モデル を解く。室内温度は、境界条件として与え、 壁面において均一かつ一定とする。建物、地 盤の熱伝導は、この式を解くことにより求め ることができる。

 各層の熱フラックスの入出量の差が層の� ��度上昇に寄与することから次の式で解ける� ��

ここで、

である。

 層間の熱フラックスJ _i は、次のように計算する。

これから、次のようになる。

従って、

となる。

 よって、層間の温度T _i ^* は、

となる。これより、

であるから、層iから層i+1への熱フラック� �は、次のようになる。

 これらから、次の3重対角行列で示される 連立方程式を得る。

ここで、T ^new は現在の時刻からdtだけ進んだ時刻における� ��度を表す。また、

である。

 〈本実施形態の作用効果〉
 本実施形態によると、建物等や地表面の反� ��等を反映した、都市空間の建物等の表面(壁 面と大気との境界)の温度(圧力変動)を算出す ることができる。また、算出された都市空間 の現実に即した温度等の気象データを気象モ デルに入力することにより、都市空間の影響 を反映した、都市空間を含む計算領域全体の 気象データを得ることができる。