JP2007292465A - シミュレーション方法、シミュレーションプログラムおよびシミュレーション装置 - Google Patents

Abstract

【課題】計算時間を短縮でき、様々な自然現象を高い精度で予測できるシミュレーション方法、シミュレーションプログラムおよびシミュレーション装置を提供する。
【解決手段】観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、実水滴が、任意個数の属性と初期時刻における全体空間を分割した分割空間内の位置座標とで表され、予め設定した所定の同属性を有する実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が所定確率を基準に多重度に応じた確率で衝突し、衝突して併合した際に多重度が変化する場合に、超水滴に関するデータを演算することで、任意時間後の実水滴に関するデータを出力するシミュレーション装置１であって、入力手段３と、データ記憶手段５と、微物理モデル演算手段７と、流体力学モデル演算手段９と、出力手段１１とを備えた。
【選択図】図１

Description

本発明は、粒子、液滴または水滴の時刻経過に伴う変化をシミュレーションするシミュレーション方法、シミュレーションプログラムおよびシミュレーション装置に関する。

従来、雲の形成、降雨、降雪、落雷等の自然現象を支配する極めて複雑な物理法則を解析するために、コンピュータを用いた数値シミュレーションが行われている。この数値シミュレーションによる解析の精度を向上させることで、実際に発生した自然現象をコンピュータ上で正確に再現することや、これから発生する自然現象を正確に予測することができる。

一般に、このような数値シミュレーションでは、自然現象を支配する極めて複雑な物理法則の解析を行う際に、当該自然現象を２つの過程に分けて、コンピュータによって演算を行っている。２つの過程の一方の過程は、大気の流れを扱う雲の力学過程であり、他方の過程は、雲や雨の構成要素である水滴の移動や状態の変化を扱う雲の微物理過程である。なお、これらの過程は相互に影響しあっている。

雲の力学過程のシミュレーションに関しては、従来手法である流体力学モデルを使って行われており、コンピュータの劇的な進歩により、計算の精度が急速に向上している。
雲の微物理過程のシミュレーションに関しては、雲が１立方メートル当たり、約１０^９個の膨大な数の水滴から構成されていることから、当該雲のすべての微物理過程を厳密に計算することは、現存するいかなるコンピュータを用いても実現不可能であり、今後実現される見通しもない。

それゆえ、現状では、雲の微物理過程のシミュレーションに関しては、大胆な近似モデルを用いて、数値シミュレーションを行っている。ここで、雲の微物理過程のシミュレーションに関して、具体的な従来手法（厳密なモンテカルロ法、改善したモンテカルロ法、ビン法、バルクパラメタリゼーション法）について説明する。

厳密なモンテカルロ法（以下、厳密モンテカルロ法という、非特許文献１参照）は、雲に含まれている水滴同士の衝突確率に、ランダムに発生させた数値を用いており、原理的には正確に雲の微物理過程をシミュレーションすることができるが、膨大なデータ記憶領域と計算量が必要となる。この厳密モンテカルロ法を改善したモンテカルロ法（以下、改善モンテカルロ法という、非特許文献２参照）は、膨大なデータ記憶領域を必要とせず、大幅な改善を成したものであるが、依然として膨大な計算量を必要とするものである。

ここで、これら厳密モンテカルロ法および改善モンテカルロ法による膨大な計算量について、概略の計算時間について説明する。非特許文献２によれば、５０［ｍ^３］の空間における２０分間の現象をシミュレーションするのに、文献開示当時のコンピュータで５．５時間かかっている。

そうすると、雲の形成、降水現象を計算するのに少なくとも１０^３［ｋｍ^３］＝１０^１２［ｍ^３］程度の空間における２時間程度の現象をシミュレーションする必要があるとすれば、文献開示当時のコンピュータで６．６×１０^１１時間＝７．５×１０^７年かかることになる。仮に、コンピュータの性能向上率が現状のまま１０年で１００倍速くなるとすると、後５０年経過しないと、実用的な計算コストでシミュレーションできないことになる。

ビン法（非特許文献３および非特許文献４参照）は、雲が形成される空間に存在する水滴を個別に取り扱うのではなく、分布関数として取り扱うこととし、この分布関数の変化を、水滴の属性（性質）ごとに考慮したビンモデルというモデルにモデル化して計算を行うものである。そして、このビン法は、現存するコンピュータによって、雲の形成について、十分な規模で微物理過程の数値シミュレーションを行うことが可能である。このビン法では、水滴を個別に取り扱っていないので、当該水滴の粒子性に起因する現象の正確な記述が保証されていないこととなる。

また、ビン法では、水滴を分布関数として取り扱う結果、既存のビンモデルを精密化し、考慮する水滴の属性の種類を増加させると、分布関数の次元が増加し、計算コストとデータ記憶領域とが膨大になることが予測される。仮に、水滴の属性として、水滴の半径Ｒ［ｍ］のみを対象とする場合、１次元の分布関数を取り扱えば、シミュレーションすることができる。ちなみに、分布関数とは、例えば、数密度分布関数ｆ（Ｒ）であり、このｆ（Ｒ）において、ｆ（Ｒ）ｄＲは半径ＲとＲ＋ｄＲとの間の大きさを持つ水滴の数と定義されるものである。

ここで、ビン法によって、水滴の複数の属性を取り扱って、より精密にシミュレーションする場合について説明する。例えば、水滴の半径Ｒの他に、水滴の速度（ｘ，ｙ，ｚ方向の３成分）、水滴に溶けているＮａＣｌ等の凝結核の質量、水滴の温度、水滴に帯電している帯電荷の７つの属性について取り扱う状況を想定する。ビン法では、属性数を７に増加させると、原理的に７次元の分布関数を取り扱う必要が生じる。この７次元の分布関数を取り扱うことは、１次元の分布関数に比べ、データ記憶領域が６乗倍、計算時間が１２乗倍になることを意味している。

一般的にビン法では、ｄ次元の分布関数を取り扱う場合を想定する。そして、分布関数の１次元当たりのビンの幅に比例する微小量パラメータをεとし、この微小量パラメータεは、シミュレーションの精度を表すものであり、この微小量パラメータεの値が小さい方がシミュレーションの精度が高くなる。そして、ビン法では、データ記憶領域は（１／ε）^ｄに比例し、計算時間は（１／ε）^２ｄに比例することになる。これにより、分布関数の次元ｄの増加と共に急速に計算時間が増加し、シミュレーションが困難になることが予測される。さらに、ビン法では、水滴の状態のみならず、水滴（水）の固体相である雪、あられ、雹（ひょう）等の状態も取り扱うとすれば、属性の数はさらに増加し、困難さが一層増してしまうことになる。

バルクパラメタリゼーション法（非特許文献５および非特許文献６参照）は、雲の力学過程と雲の微物理過程とを連結させて、雲の形成、降雨等の自然現象をシミュレーションする現在主流の方法である。このバルクパラメタリゼーション法の特徴は、雲の微物理過程を大幅に簡略化したパラメータとして表現し、現象を近似的に再現するように当該パラメータを調整して、雲の力学過程に雲の微物理過程を取り込むことである。それゆえ、このバルクパラメタリゼーション法では、雲の状態変化を直接計算することができず、多様に変化する自然現象（想定外の気象状況）について、高い精度で予測することは困難である。

このように、雲の微物理過程の扱いについては、気象学、気候学の研究において、重要な課題の一つとなっている。さらに、ビン法による雲の微物理過程と、雲の力学過程とを連結する方法（非特許文献７参照）も模索されている。
D.T.Gillespie,"An Exact Method for Numerically Simulating the Stochastic C-oalescence Process in a Cloud,"J.Atoms.Sci.,32,1977(1975) M.SeeB「ドイツ語のエスツェット」elberg,T.Trautmann,and M.Thorn,"Stochastic simulations as a benchmark for mathematical methods solving the coalescence equation,"Atmos.Res.,40,33(1996) A.Bott,"A Flux Method for the Numerical Solution of the Stochastic Collection Equation,"J.Atoms.Sci.,55,2284(1998) A.Bott,"A Flux Method for the Numerical Solution of the Stochastic Collection Equation:Extension to Two-Dimensional Particle Distributions,"J.Atoms.Sci.,57,284(2000) E.Kessler,"On the Distribution and Continuity of Water Substance in Atomspheric Circulations,"Met.Monograph,Vol.10,No.32,American Meteorological Society,Boston,84 pp M.Murakami,"Numerical Modeling of Dynamical and Microphysical Evolution of an Isolated Convective Cloud,"J.Meteor.Soc.Japan,68,107(1990) B.H.Lynn,et al.,"Spectral(Bin) Microphysics Coupled with a Mesoscale Model(MM5).Part I:Model Description and First Results,"Mon.Wea.Rev.,133,44(2005)

しかしながら、前記した従来技術にはそれぞれ問題がある。すなわち、厳密モンテカルロ法および改善モンテカルロ法では、計算時間がかかりすぎるという問題があり、ビン法では、水滴を分布関数として取り扱うことによる不確かさや、属性数の増加に伴う計算時間の増加による拡張性の低さ等、多くの問題がある。また、バルクパラメタリゼーション法には、雲の微物理過程を大幅に簡略化してしまうため、計算時間を短縮できるものの、多様に変化する自然現象（想定外の気象状況）について、高い精度で予測することができないという問題がある。

そこで、本発明では、前記した問題を解決し、計算時間を短縮でき、水滴等の対象を分布関数として取り扱うことなく、当該対象の属性数が増加しても計算時間の増加を抑制することができ、様々な自然現象を高い精度で予測することができるシミュレーション方法、シミュレーションプログラムおよびシミュレーション装置を提供することを目的とする。

前記課題を解決するため、請求項１に記載のシミュレーション方法は、観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するシミュレーション方法であって、入力ステップと、演算ステップと、出力ステップと、を含む手順とした。

かかる手順によれば、シミュレーション方法は、入力ステップにおいて、初期時刻と超粒子の属性と超粒子の総個数と体積と超粒子の速度と超粒子の位置座標と流体場変数とを初期変数として入力する。続いて、シミュレーション方法は、演算ステップにおいて、入力ステップにて入力された初期変数に基づいて、経過時間に伴った実粒子の運動が体積と速度と位置座標と流体場変数とに従って属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する。そして、シミュレーション方法は、出力ステップにおいて、演算ステップにより超粒子の属性、速度、位置座標、個数および多重度の演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する。

請求項２に記載のシミュレーションプログラムは、観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するために、コンピュータを、入力手段、演算手段、出力手段、として機能させる構成とした。

かかる構成によれば、シミュレーションプログラムは、入力手段によって、初期時刻と超粒子の属性と超粒子の総個数と体積と超粒子の速度と超粒子の位置座標と流体場変数とを初期変数として入力する。続いて、シミュレーションプログラムは、演算手段によって、入力手段で入力された初期変数に基づいて、経過時間に伴った実粒子の運動が体積と速度と位置座標と流体場変数とに従って属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する。そして、シミュレーションプログラムは、出力手段によって、演算手段により超粒子の属性、速度、位置座標、個数および多重度の演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する。

請求項３に記載のシミュレーション装置は、観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するシミュレーション装置であって、入力手段と、演算手段と、出力手段と、を備える構成とした。

かかる構成によれば、シミュレーション装置は、入力手段によって、初期時刻と超粒子の属性と超粒子の総個数と体積と超粒子の速度と超粒子の位置座標と流体場変数とを初期変数として入力する。続いて、シミュレーション装置は、演算手段によって、入力手段で入力された初期変数に基づいて、経過時間に伴った実粒子の運動が体積と速度と位置座標と流体場変数とに従って属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する。そして、シミュレーション装置は、出力手段によって、演算手段により超粒子の属性、速度、位置座標、個数および多重度の演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する。

請求項４に記載のシミュレーション方法は、観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、前記超水滴に関するデータを演算することで、任意時間後の前記実水滴に関するデータを出力するシミュレーション方法であって、変数入力ステップと、微物理モデル演算ステップと、流体力学モデル演算ステップと、出力ステップと、を含む手順とした。

かかる手順によれば、シミュレーション方法は、変数入力ステップにおいて、初期時刻と超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積と超水滴の位置座標と当該分割空間ごとの前記実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する。続いて、シミュレーション方法は、微物理モデル演算ステップにおいて、超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積とに基づいて、全体空間内における超水滴の運動による位置変化、超水滴の凝結成長による水量変化および超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る。なお、相互作用量は、空間単位体積当たり液体の水の質量と、単位時間当たりに蒸発した空間単位体積当たりの水の質量とを指している。また、シミュレーション方法は、流体力学モデル演算ステップにおいて、微物理モデル演算ステップにて得られた相互作用量および周辺環境データに基づいて、実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を微物理モデル演算ステップにフィードバックする。そして、シミュレーション方法は、出力ステップにおいて、流体力学モデル演算ステップと微物理モデル演算ステップとによる演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を任意時間後の実水滴に関するデータとして出力すると共に、任意時間後の周辺環境データを出力する。

請求項５に記載のシミュレーションプログラムは、観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、前記超水滴に関するデータを演算することで、任意時間後の前記実水滴に関するデータを出力するために、コンピュータを、変数入力手段、微物理モデル演算手段、流体力学モデル演算手段、出力手段、として機能させる構成とした。

かかる構成によれば、シミュレーションプログラムは、変数入力手段によって、初期時刻と超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積と超水滴の位置座標と当該分割空間ごとの実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する。続いて、シミュレーションプログラムは、微物理モデル演算手段によって、超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積とに基づいて、全体空間内における超水滴の運動による位置変化、超水滴の凝結成長による水量変化および超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る。また、シミュレーションプログラムは、流体力学モデル演算手段によって、微物理モデル演算手段で得られた相互作用量および周辺環境データに基づいて、実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を微物理モデル演算手段にフィードバックする。シミュレーションプログラムは、出力手段によって、流体力学モデル演算手段と微物理モデル演算手段とによる演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を任意時間後の実水滴に関するデータとして出力すると共に、任意時間後の周辺環境データを出力する。

請求項６に記載のシミュレーションプログラムは、請求項５に記載のシミュレーションプログラムにおいて、前記微物理モデル演算手段が、超水滴運動演算手段と、超水滴凝結成長演算手段と、超水滴衝突併合演算手段と、を有することを特徴とする。

かかる構成によれば、シミュレーションプログラムは、微物理モデル演算手段の超水滴運動演算手段によって、超水滴にかかる重力および空気抵抗が釣り合った状態で、超水滴の運動が風速によって変化し、当該風速に対して一定の相対速度である終端速度で運動するものとして、当該終端速度を演算する。また、シミュレーションプログラムは、微物理モデル演算手段の超水滴凝結成長演算手段によって、超水滴に含まれる水量が周辺環境データに含まれる湿度によって変化するものとして、当該水量を演算する。そして、シミュレーションプログラムは、微物理モデル演算手段の超水滴衝突併合演算手段によって、超水滴同士が衝突する任意のペアを設定し、当該ペアが衝突併合する確率を衝突併合確率とし、ペアの数を所定数に削減すると共に衝突併合確率を所定幅上昇させて、超水滴同士の衝突併合による超水滴の属性、多重度および個数を演算する。

請求項７に記載のシミュレーションプログラムは、請求項６に記載のシミュレーションプログラムにおいて、超水滴衝突併合演算手段が、前記衝突併合の過程をモンテカルロ法による数値シミュレーションに従って演算することを特徴とする。

かかる構成によれば、シミュレーションプログラムは、超水滴衝突併合演算手段によって演算する衝突併合の過程にモンテカルロ法を用いることで、分割空間内の超水滴について、実際に乱数を発生させて衝突併合させる。

請求項８に記載のシミュレーション装置は、観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、前記超水滴に関するデータを演算することで、前記任意時間後の実水滴に関するデータを出力するシミュレーション装置であって、変数入力手段と、微物理モデル演算手段と、流体力学モデル演算手段と、出力手段と、を備える構成とした。

かかる構成によれば、シミュレーション装置は、変数入力手段によって、初期時刻と超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積と超水滴の位置座標と当該分割空間ごとの実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する。続いて、シミュレーション装置は、微物理モデル演算手段によって、超水滴の属性と超水滴の総個数と全体空間の体積と分割空間の体積とに基づいて、全体空間内における超水滴の運動による位置変化、超水滴の凝結成長による水量変化および超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る。また、シミュレーション装置は、流体力学モデル演算手段によって、微物理モデル演算手段で得られた相互作用量および周辺環境データに基づいて、実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を微物理モデル演算手段にフィードバックする。シミュレーション装置は、出力手段によって、流体力学モデル演算手段と微物理モデル演算手段とによる演算を、任意時間に達するまで繰り返した後に、繰り返した後の結果を任意時間後の実水滴に関するデータとして出力すると共に、任意時間後の周辺環境データを出力する。

請求項９に記載のシミュレーション装置は、請求項８に記載のシミュレーション装置において、前記微物理モデル演算手段が、超水滴運動演算手段と、超水滴凝結成長演算手段と、超水滴衝突併合演算手段と、を有することを特徴とする。

かかる構成によれば、シミュレーション装置は、微物理モデル演算手段の超水滴運動演算手段によって、超水滴にかかる重力および空気抵抗が釣り合った状態で、超水滴の運動が風速によって変化し、当該風速に対して一定の相対速度である終端速度で運動するものとして、当該終端速度を演算する。また、シミュレーション装置は、微物理モデル演算手段の超水滴凝結成長演算手段によって、超水滴に含まれる水量が周辺環境データに含まれる湿度によって変化するものとして、当該水量を演算する。そして、シミュレーション装置は、微物理モデル演算手段の超水滴衝突併合演算手段によって、超水滴同士が衝突する任意のペアを設定し、当該ペアが衝突併合する確率を衝突併合確率とし、ペアの数を所定数に削減すると共に衝突併合確率を所定幅上昇させて、超水滴同士の衝突併合による超水滴の属性、多重度および個数を演算する。

請求項１０に記載のシミュレーション装置は、請求項９に記載のシミュレーション装置において、超水滴衝突併合演算手段が、前記衝突併合の過程をモンテカルロ法による数値シミュレーションに従って演算することを特徴とする。

かかる構成によれば、シミュレーション装置は、超水滴衝突併合演算手段によって演算する衝突併合の過程にモンテカルロ法を用いることで、分割空間内の超水滴について、実際に乱数を発生させて衝突併合させる。

請求項１、２、３に記載の発明によれば、空間内のある体積内で同じ属性を有している複数の実粒子を超粒子として取り扱って演算することで、実粒子を分布関数として取り扱うことなく、当該実粒子の属性数が増加しても計算時間の増加を抑制することができ、様々な自然現象を高い精度で予測することができる。

請求項４、５、８に記載の発明によれば、分割空間内で同じ属性を有している複数の実水滴を超水滴として取り扱って、超水滴モデルを雲の微物理モデルとして用いることで、実水滴を分布関数として取り扱うことなく、当該実水滴の属性数が増加しても、雲の形成、降雨等の計算時間の増加を抑制することができ、高い精度で予測することができる。

請求項６、９に記載の発明によれば、超水滴について、超水滴の運動、超水滴の凝結成長および超水滴の衝突併合の演算を行っているので、より高い精度で実水滴の経過時間に伴う変化を予測することができる。

請求項７、１０に記載の発明によれば、超水滴の衝突併合の演算について、モンテカルロ法を用いることで、計算時間を大幅に短縮することができる。

次に、本発明の実施形態について、適宜、図面を参照しながら詳細に説明する。
（シミュレーション装置の構成）
図１はシミュレーション装置のブロック図である。この図１に示すように、シミュレーション装置１は、実粒子の一形態である水滴に着目し、同属性を有する複数の水滴を超水滴として取り扱う超水滴モデルを用いて、観測される全体空間における経過時間に伴う当該水滴の変化をシミュレーションすることで、雲の形成、降雨等の自然現象を予測するもので、入力手段（変数入力手段）３と、データ記憶手段５と、微物理モデル演算手段７と、流体力学モデル演算手段９と、出力手段１１とを備えている。

まず、このシミュレーション装置１で取り扱う超水滴および超水滴モデルについて説明する。超水滴は、１個１個の実水滴が持つ物理量、例えば、実水滴の大きさ、凝結核の量および種類、電荷、速度、温度等を「実水滴の属性」とし、この場合に、同じ属性を持つ任意数の実水滴の集合としている。つまり、１個の超水滴には、任意数の実水滴が含まれており、この任意数を多重度ｎとして表現している。なお、このｎは、整数０，１，２，３，・・・である。例えば、この多重度ｎが５、１００、４００，０００となれば、これらの多重度ｎの超水滴を計算するということは、実水滴について個々に計算するよりも、単純に言えば、計算時間が１／５、１／１００、１／４００，０００となることを意味している。つまり、多重度ｎを大きくとればとるほど計算時間は短くすることができる。

そして、各超水滴は、実水滴と同様に、風や重力に従って自由に全体空間内を移動し続ける。ある超水滴が時刻ｔに全体空間上のある点「ｒ」（ｔ）［ｍ］に存在するとし、全体空間が適当な大きさの格子（分割空間）に分割されているとする。なお、この「ｒ」（ｔ）［ｍ］は、ｘｙｚ座標の３成分を指定することで決定され、ベクトルｒ＝（ｘ，ｙ，ｚ）となるが、ここでは、ベクトルｒを、カギ括弧を用いて「ｒ」で表すこととする。超水滴で表現されるｎ個の実水滴は、「ｒ」（ｔ）［ｍ］を含む分割空間内に一様ランダムに分布していると解釈する。これらの対応関係を模式的に図２に示す。

図２（ａ）に示すように、多重度３の超水滴が分割空間ａから分割空間ｂに移動した場合、図２（ｂ）に示すように、３個の実水滴が分割空間ａから分割空間ｂに移動したことと同じことになる。ただし、図２（ａ）と図２（ｂ）とを比較してもわかるように、実水滴の場合は、３個の実水滴が移動前には分割空間ａに、移動後には分割空間ｂに収まっていればよく、分割空間ａ，ｂ内であれば、各実水滴の実際の座標位置は問題がない。

超水滴モデルは、このような超水滴が全体空間内を自由に移動（運動）し、ある確率（詳細は後記する）で衝突して結合（衝突併合）することで、新たな超水滴となることを少なくとも想定したものである。この実施の形態では、この超水滴モデルに、超水滴同士が衝突併合する過程のみならず、超水滴が凝結成長する過程を加味しており、さらに超水滴が運動する過程について計算時間を短縮する工夫が施されている。ここで、これら超水滴モデルにおける超水滴の運動、超水滴の凝結成長、超水滴同士の衝突併合について概略を説明する（詳細な説明は微物理モデル演算手段７において行う）。

超水滴の運動は、全体空間内に存在する実水滴が重力および大気から空気抵抗を受けており、通常の運動方程式に従うのと同様に、当該運動方程式に従って時間発展する。
超水滴の凝結成長は、全体空間内に存在する実水滴が凝結核の量と周辺の湿度に応じて大気中の水蒸気を吸収または放出する法則に従うのと同様に、この法則に従って直接計算することができる。

超水滴同士の衝突併合は、多重度ｎ_ｊを持つ超水滴と多重度ｎ_ｋ（＞ｎ_ｊ）を持つ超水滴とが衝突併合する場合を想定する。この場合、多重度ｎ_ｋ−ｎ_ｊを持つ超水滴と多重度ｎ_ｊを持つ超水滴とが生成されると解釈する。

さらに、多重度ｎ_ｊを持つ超水滴と多重度ｎ_ｋ（＝ｎ_ｊ）を持つ超水滴とが衝突併合する場合を想定する。この場合、多重度ｎ_ｋ−ｎ_ｊ＝０を持つ超水滴と多重度ｎ_ｊを持つ超水滴とが生成されると解釈するのではなく、多重度［ｎ_ｊ／２］を持つ超水滴と多重度ｎ_ｊ−［ｎ_ｊ／２］を持つ超水滴とが生成されると解釈する。ここで、［ｎ_ｊ／２］はｎ_ｊ／２個を越えない最大整数である。以下、この明細書中“［］”と記載した場合、この［］はガウス記号を表すこととする（数式中は“”を付加せず［］のまま使用する）。

例を挙げると、多重度６を持つ超水滴同士が衝突併合すると、６個の実水滴ができるので、これを半分にして、多重度３を持つ超水滴を２個生成すると解釈する。また、多重度５を持つ超水滴同士が衝突併合すると、多重度３を持つ超水滴と多重度２を持つ超水滴とが生成されると解釈する。

この超水滴同士が衝突併合する様子を図３に示す。この図３は、（ａ）の領域に衝突前の実水滴（２個の実水滴と３個の実水滴とが衝突）、（ｂ）の領域に衝突前の超水滴（多重度２の超水滴と多重度３の超水滴とが衝突）、（ｃ）の領域に衝突後の実水滴、（ｄ）の領域に衝突後の超水滴を模式的に示している。

超水滴同士の衝突併合は、同一の分割空間内における実水滴同士の衝突する確率に基づき、乱数を発生させて、当該衝突併合を起こす過程を計算する。

なお、この超水滴は、任意属性に関し実水滴の数の分布を示す分布関数と当該超水滴の数の分布を示す分布関数とが同じ分布となるように、入力された超水滴の総個数から計算して、それぞれの超水滴の多重度が設定されたものである。

ここで、超粒子およびこの超粒子を用いた超粒子法についても説明しておく。超粒子も超水滴と同様に、任意属性に関し実粒子の数の分布を示す分布関数と当該超粒子の数の分布を示す分布関数とが同じ分布となるように、入力された超粒子の総個数から計算して、それぞれの超粒子の多重度が設定されたものである。この多重度を持つ超粒子は、実際に存在している実粒子の分布関数を、近似して求めたものであるので、当該超粒子を用いて、当該実粒子の分布関数（元の分布関数という）を完全に再現することができない。仮に、元の分布関数を再現しようとすると、元の分布関数は任意属性の属性値に関するε程度の幅を持って再現されることになる。このεは、ビン法におけるビンの幅に対応した量である。このεを小さくするために、各超粒子の多重度を減少させ、超粒子の個数を増加させる必要がある。

なお、すべての超粒子の多重度が１になった時点で、実際に存在している実粒子と同じ状態となり、超粒子法は、真に厳密な自然現象のシミュレーションを行うことが可能となる。この場合、εは有限の値ではあるが、０に限りなく近づくことになる。このεが有限の値に留まることは、実粒子の個数が非常に多くなったとしても当該実粒子の個数が無限には多くないことを意味している。

次に、このシミュレーション装置１でシミュレーションする空間である全体空間について説明する。
全体空間は、自然現象が観測される空間、自然現象の動向を予測したい空間の約数キロ範囲を設定することができる。この全体空間は、当該全体空間の体積を特定可能なｘ軸、ｙ軸およびｚ軸方向の長さを備えており、任意の大きさを有することができる。なお、この任意の大きさとは、全体空間が数キロ範囲に限定されるものではなく、例えば、数ミリ四方の空間であるとか、数千キロの広範囲に亘る空間であってもよい。また、全体空間の形状は、直方体に限定されず、任意の形状であってもよい。

さらに、この全体空間は、便宜上、水平方向の空間の端が連続してつながっているとみなす仮想的な空間を想定している。つまり、全体空間の右側の境界から出て行った実粒子は、左側の境界から入ってくる空間である。しばしば、このような仮想的な空間の規定は、周期境界条件と呼ばれる。

そして、この全体空間を任意の大きさに分割した空間を分割空間（空間格子ともいう）としている。この分割空間も、全体空間と同様に当該分割空間の体積を特定可能なｘ軸、ｙ軸およびｚ軸方向の長さを備えている。なお、この実施の形態では、分割空間を直方体としており、当該直方体の１つの頂点を格子点としている。

なお、超粒子および超粒子モデルは、超水滴および超水滴モデルの適用範囲を実水滴に限定せずに拡張したもので、超粒子同士が衝突した際に結合するものであれば、任意組成の実液滴や実粒子であってもよい。そして、シミュレーション装置１は、任意組成の実液滴や実粒子の経過時間に伴う変化をシミュレーションすることができる。

入力手段３は、超水滴に関するデータと、全体空間に関するデータと、計算時間に関するデータと、周辺環境に関するデータ（周辺環境データ）と、記憶域に関するデータとを、初期変数として入力するものである。

超水滴に関するデータは、超水滴の属性（最も簡略化した場合、超水滴の半径Ｒと凝結核の質量のみ）と、超水滴の総個数Ｎとである。この実施の形態では、超水滴の属性は、超水滴の半径と超水滴の凝結核の質量とに限定している。なお、ここでいう超水滴の属性には、超水滴の速度ｖと、超水滴の位置ｒと、超水滴の多重度ｎとは含まれていないこととしている。

具体的には、超水滴に関するデータは、各超水滴に含まれる水滴の半径Ｒと、各超水滴に含まれる水滴の速度「ｖ」＝（ｖｘ，ｖｙ，ｖｚ）と、各超水滴中の水滴に含まれる凝結核の質量Ｍと、各超水滴の位置「ｒ」＝（ｘ，ｙ，ｚ）と、各超水滴の多重度ｎである。ただし、この実施の形態では、これらの超水滴に関するデータを全て入力する必要はなく、超水滴の総個数のみを設定すれば、これらのデータが微物理モデル演算手段７において、設定されるように構成されている。

全体空間に関するデータは、水平方向の領域の長さ（全体空間の幅、全体空間の奥行き）［ｋｍ］と、垂直方向の領域の長さ（全体空間の高さ）［ｋｍ］と、水平方向の空間格子間隔（分割空間の幅、分割空間の奥行き）［ｍ］と、垂直方向の空間格子間隔（分割空間の高さ）［ｍ］とである。なお、これらの全体空間に関するデータを用いて計算することができるが、これらの他に、水平方向の格子の数（水平方向の分割空間の数）と、垂直方向の格子の数（垂直方向の分割空間の数）とを入力してもよい。

計算時間に関するデータは、シミュレーションを開始する時刻である初期時刻［ｓｅｃ］と、経過時間に伴う変化を出力可能な時間ステップ（刻み時刻）［ｓｅｃ］と、シミュレーションを実行する時間である総計算時間［ｓｅｃ］とである。

周辺環境に関するデータは、水滴が存在している全体空間の状態、すなわち、大気の状態を表す変数である。具体的には、風速「Ｕ」＝（Ｕ，Ｖ，Ｗ）と、相対湿度Ｓ（Ｓ＝１で湿度１００％を表す）と、湿潤大気密度ρ＝ρ_ｄ＋ρ_ｖ（ρ_ｄ：乾燥大気の密度、ρ_ｖ：水蒸気の密度）と、水蒸気の混合比ｑ_ｖ＝ρ_ｖ／ρと、気温Ｔと、温位θと、Π＝（Ｐ／Ｐ_０）^{（Ｒｄ／ｃｐ）}（Π：エクスナー関数、Ｐ_０：基準気圧１０００［ｈＰａ］、Ｒ_ｄ：乾燥空気の気体定数、ｃ_ｐ：定圧比熱）と、空間単位体積当たり液体の水の質量ρ_ωと、単位時間当たりに蒸発した空間単位体積当たりの水の質量ρＳ_ｖとである。

記憶域に関するデータは、データ記憶手段５において、入力手段３で入力したデータや計算途中であるデータを記憶する記憶域を指定するものである。記憶域に関するデータで指定されるものには、例えば、Ｎ個の個数を持つ配列として、超水滴の終端速度を保存する記憶域、超水滴のｘ座標を保存する記憶域、超水滴のｙ座標を保存する記憶域、超水滴のｚ座標を保存する記憶域、超水滴の半径を保存する記憶域、超水滴の凝結核の質量を保存する記憶域、超水滴の多重度を保存する記憶域がある。

また、記憶域に関するデータで指定されるものには、流体力学モデル（大気流体場を一般的な流体力学モデルとしてもの）における風速ｘ成分（ｘ軸方向の成分）（Ｕを指す）を保存する記憶域、流体力学モデルにおける風速ｙ成分（ｙ軸方向の成分）（Ｖを指す）を保存する記憶域、流体力学モデルにおける風速ｚ成分（ｚ軸方向の成分）（Ｗを指す）を保存する記憶域、流体力学モデルにおける温位θを保存する記憶域、流体力学モデルにおける大気密度ρを保存する記憶域、流体力学モデルにおける水蒸気混合比ｑ_ｖを保存する記憶域、流体力学モデルにおける圧力ｐを保存する記憶域、流体力学モデルにおける温度Ｔを保存する記憶域、流体力学モデルにおける相対湿度Ｓを保存する記憶域がある。

データ記憶手段５は、入力手段３で入力された初期変数（各種のデータ）と、微物理モデル演算手段７および流体力学モデル演算手段９で計算途中であるデータとを記憶し、必要に応じて適宜出力するものである。このデータ記憶手段５は、一般的なメモリやハードディスク等によって構成されている。

微物理モデル演算手段７は、入力手段３で入力された初期変数の超水滴に関するデータに基づいて、予め設定されている超水滴モデルを雲の微物理モデルとして用いて、時間ステップ（刻み時刻）ごとに、雲の微物理過程をシミュレーションするもので、超水滴運動演算手段７ａと、超水滴凝結成長演算手段７ｂと、超水滴衝突併合演算手段７ｃとを備えている。なお、この微物理モデル演算手段７に、超水滴に関するデータが入力されると、超水滴モデルに基づいて、各超水滴の多重度が決定され、各超水滴を識別する番号が振られる。

超水滴運動演算手段７ａは、水滴が重力と大気からの空気抵抗とを受けて一般的な運動方程式に従って運動するのと同様に、超水滴も当該運動方程式に従って運動することとして当該超水滴の運動を演算するものである。この運動方程式は、次に示す数式（１）のように与えられる。

この数式（１）の前段において、「ｇ」（ここでは、ｇ＝（０，０，ｇ））は重力定数、「Ｆ_Ｄ」は空気抵抗、「Ｕ」は風速を示している。また、数式（１）の後段において、「ｒ」は水滴の位置座標を、「ｖ」は水滴の速度を示している。まず、この数式（１）を、水滴に適用した場合について説明する。実際には、水滴は、当該水滴の半径が非常に小さいため、当該水滴が大気から受ける空気抵抗と重力との大きさが直ちに釣り合って打ち消し合い、風速に対して一定の相対速度で運動することになる。この相対速度を終端速度とし、水滴は常にこの終端速度で運動するものとし、さらに、この終端速度を参考文献１で開示されていることに従って、水滴の半径の関数として与えることとする。参考文献１は、K.V.Beard,”Terminal Velocity and Shape of Cloud and Precipitation Drops Aloft,”J.Atoms.Sci.,33,851(1976).である。

そうすると、水滴の速度ｖ［ｍ／ｓ］は、当該水滴の半径Ｒ［ｍ］と風速「Ｕ」とを与えることで決定することになる。図４に水滴の半径と終端速度（正確には、終端速度のｚ成分、なお、終端速度のｘ，ｙ成分は０［ｍ／ｓ］）との関係を示した。この図４に示したように、１０［μｍ］程度の半径を持つ水滴の終端速度はほとんど０［ｍ／ｓ］に等しいために、結果として落ちてこない（雨にならない）こと、すなわち、雲となっていることがわかる。また、１［ｍｍ］程度の半径を持つ水滴の終端速度は１０［ｍ／ｓ］弱になっているために、雨となって落ちてくることがわかる。正確には、上昇気流、すなわち、風速「Ｕ」のｚ成分Ｗが水滴の終端速度よりも小さいならば、当該水滴は雨となって落ちてくることになる。

超水滴の速度も、この水滴の場合と同様になり、半径Ｒ［ｍ］と風速「Ｕ」とで与えられる。さらに、超水滴の位置座標も、この水滴の場合と同様になり、超水滴の速度「ｖ」で与えられる。また、数式（１）において、ｍは水滴に含まれる水の質量であり、この水の質量は、（４π／３）Ｒ^３ρ_ｌｉｑ（ρ_ｌｉｑ＝１［ｇ／ｃｃ］は水の密度を表す）で与えられる。また、水滴に含まれる水の質量ｍは、水滴の凝結核の質量Ｍ（Ｍ＝１．０×１０^−１６［ｇ］）に比べて十分に大きいので、この凝結核の質量Ｍが水滴の運動に及ぼす影響を無視しても問題がない。なお、終端速度とは、重力ｍｇと空気抵抗Ｆ_ｄとが釣り合った状態の速度である。つまり、−ｍｇ＝Ｆ_ｄが成立している場合の速度である。

そして、超水滴の終端速度をｖｔ、風速「Ｕ」＝（Ｕ，Ｖ，Ｗ）とすると、超水滴の速度「ｖ」＝「Ｕ」−（０，０，ｖｔ）となる。つまり、ｖｘ＝Ｕ、ｖｙ＝Ｖ、ｖｚ＝Ｗ−ｖｔとなる。なお、この超水滴の終端速度ｖｔをＲの関数としてグラフ化したものが図４である。

超水滴凝結成長演算手段７ｂは、水滴に含まれる水量（ここでは、当該水量を求めるのに水滴の半径に着目している）が周辺の湿度に応じて大気中の水蒸気を吸収または放出して変化するのと同様に、超水滴に含まれる水量（超水滴の半径）の変化、すなわち、凝結成長の演算を行うものである。一般に、経過時間に伴う水滴の凝結成長は、次に示す数式（２）のように与えられる。

この数式（２）において、Ｓは相対湿度であり、Ｆ_ｋは水滴の熱拡散に関する項であり、このＦ_Ｋは気温に依存し、Ｆ_ｄは水滴の水蒸気拡散に関する項であり、このＦ_ｄも気温に依存する。また、ａの項は水滴の表面張力によって凝結成長が抑制される効果を表しており、このａの項も気温に依存し、ｂの項は凝結核（ＣＣＮ）の溶解効果によって凝結成長が促進される効果を表しており、このｂの項は凝結核の質量Ｍに依存する。なお、この数式（２）は参考文献２によって得られる。この参考文献２は、A Short Course in Cloud Physics,Third Edition,R.R.Rogers & M.K.Yan(Butter worth Heinemann)である。

そして、この超水滴凝結成長演算手段７ｂは、この数式（２）に従って、超水滴の凝結成長を演算する。

超水滴衝突併合演算手段７ｃは、超水滴同士がする衝突して結合する（衝突併合）による超水滴の属性、多重度および個数を演算するものである。この超水滴衝突併合演算手段７ｃでは、後記するモンテカルロ法による数値シミュレーションによって、超水滴の衝突併合の演算を実行している。

まず、ある空間にある水滴同士（粒子同士）が確率的に衝突することについて説明する。体積ΔＶの空間中を２つの水滴が飛んでいることを想定し、水滴ｉ＝１，２とし、速度「ｖ_ｉ」と半径Ｒ_ｉとを有しているとする。この場合時間Δｔの間に２つの水滴が実効的にスウィープする体積（接触する可能性がある範囲、近くを通過する可能性がある範囲）は、次に示す数式（３）のように与えられる。

この数式（３）で表される状態を、図５に模式的に示す。この図５に示すように、ΔＶにおける２つの水滴の位置が一様にランダム、すなわち、ΔＶのどこに存在する確率も等しいとするならば、２つの水滴がΔｔの間に衝突する確率は、実効的にスウィープする体積中に他方の水滴が存在する確率に等しくなり、次に示す数式（４）のように与えられる。

そして、十分に小さい体積中に、複数の水滴が存在している場合には、当該水滴同士が同様の確率で衝突すると想定される。しかし、この確率は、水滴の体積（または半径）が小さくなると、単純に求めることができなくなる場合が生じる。水滴は元々微小であり、凝結成長か衝突併合しなければ、微小のままであり、慣性が小さいため、衝突併合の対象となる他の水滴を回り込むようにして動いたり、衝突したとしても併合せずに跳ね返したりする場合があるからである。

このような状態を考慮して、水滴ｊと水滴ｋとの衝突併合確率Ｐ_ｊｋは次に示す数式（５）のように与えられる。

この数式（５）において、Ｅ_ｊｋは衝突効率と呼ばれており、［０，１］の間の実数値をとるＲ_ｊとＲ_ｋの関数である。なお、この衝突効率Ｅ_ｊｋについては、理論的・実験的に詳細に求められており、参考文献群３に従って与えることとしている。参考文献群３は、M.H.Davis,J.Atoms.Sci.,29,911(1972). W.D.HALL,”A DETAILED MICROPHYSICAL MODEL WITHIN A TWO-DIMENSIONAL DYNAMIC FRAMEWORK-MODEL DESCRIPTION AND PRELIMINARY-RESULTS,”37,2486(1980) P.R.Jonas,”COLLISION EFFICIENCY OF SMALL DROPS,”QUART.J.
ROY.METEOR.SOC.98,681(1972)

このように水滴同士の衝突併合確率が明確化したので、超水滴同士の衝突併合確率への適用について説明する。まず、超水滴同士の衝突併合については、前記したように多重度に応じて、新たな超水滴が生成されることになる。この場合の超水滴同士の衝突併合確率を求めることとする。

まず、同一の属性を有しているｎ_ｊ個の実水滴と、別の同一の属性を有しているｎ_ｋ個の実水滴とが衝突すると想定する。そうすると、全部でｎ_ｊｎ_ｋ個のペア（組み合わせ）を作ることができ、それぞれのペアは、数式（５）に示した衝突併合確率Ｐ_ｊｋで衝突併合する。この場合の実水滴の総衝突併合数の期待値は、ｎ_ｊｎ_ｋＰ_ｊｋである。

このことを前提に、ｎ_ｊ個の多重度を持つ超水滴と、ｎ_ｋ個の多重度を持つ超水滴とが衝突併合する状況を想定する。そうすると、衝突併合後、多重度｜ｎ_ｋ−ｎ_ｊ｜の超水滴と、多重度ｍｉｎ（ｎ_ｋ，ｎ_ｊ）の超水滴とが生成されることとなり、これは、ｍｉｎ（ｎ_ｋ，ｎ_ｊ）回の実水滴の衝突併合が起こったことを意味する。そして、超水滴の衝突併合確率をＰ_ｊｋ ^（ｓ）で表すと、次に示す数式（６）のように与えられる。なお、ｎ_ｋ＝ｎ_ｊの場合［ｎ_ｊ／２］とｎ_ｊ−［ｎ_ｊ／２］の超水滴が生成されるが、実水滴の衝突併合回数は、やはりｍｉｎ（ｎ_ｋ，ｎ_ｊ）である。それゆえ、ｎ_ｋ＝ｎ_ｊの場合にも数式（６）の形は変わらないこととなる。

そして、この超水滴の衝突併合確率Ｐ_ｊｋ ^（ｓ）によって、実水滴の総衝突併合数の期待値が再現されることが、次に示す数式（７）により確認される。

超水滴衝突併合演算手段７ｃでは、超水滴の衝突併合確率をＰ_ｊｋ ^（ｓ）が求まったので、これに基づいて、超水滴の衝突併合過程をモンテカルロ法による数値シミュレーションを実行する。すなわち、乱数を発生させて、この乱数に従って同一分割空間内の超水滴同士を衝突併合させる。

超水滴衝突併合演算手段７ｃにおいて、時刻ｔにおける超水滴の衝突併合の演算手法を以下に説明する。
超水滴衝突併合演算手段７ｃは、まず、時刻ｔにおける各超水滴の位置を調べ、各超水滴を分割空間ごとに予め分類しておく。ある分割空間に超水滴がＮ個存在しているとし、各超水滴ｌ＝１，２，・・・，Ｎはそれぞれｎ_ｌ個（ｎは任意の整数）の多重度を持つとする。

超水滴衝突併合演算手段７ｃは、これら全ての超水滴について、２つの対（ｊとｋ）を、乱数を用いて選択する。この対を衝突併合候補対（ペア）とする。分割空間内に奇数の超水滴が存在しているとすると、衝突併合候補対に含まれない超水滴が１個残ることになるが、この場合Ｎ／２個を越えない最大整数 “［Ｎ／２］”個のペアができる。

この衝突併合対の作り方についてより詳細に説明する。まず、Ｓ＝（１，２，・・・，Ｎ）という順列をランダムに並べ替え、並べ替えた順列を、Ｓ´＝（ａ１，ａ２，・・・，ａＮ）、ａ_ｋ＝１，２，・・・，Ｎとする。この並べ替えた順列Ｓ´中の数字は重複しない。つまり、ｋ≠ｋ´ならば、ａ_ｋ≠ａ_ｋ´である。そうすると、Ｓ´は全部でＮ！通りの種類が考えられる。そして、Ｎ個の要素からなる順列Ｓをランダムに並べ替えることを、Ｎ！個の順列が等しい確率１／Ｎ！で出現する操作であると定義する。

なお、この超水滴衝突併合演算手段７ｃでは、数値計算上、超水滴の個数Ｎが増加しても、当該個数Ｎに比例した計算コストでこの操作を実現することができる。そして、超水滴衝突併合演算手段７ｃは、並べ替えた順列Ｓ´を使って、並べ替えた順に衝突併合候補対を作っていく。すなわち、衝突併合候補対は、（ａ１，ａ２）、（ａ３，ａ４）、（ａ５，ａ６）・・・となる。

また、分割空間内の超水滴同士は、本来、当該超水滴同士の全ての組み合わせについて衝突する可能性があるため、Ｎ（Ｎ−１）／２個の全てのペアの組み合わせに対して、衝突併合する可能性を判定する必要がある。しかし、これではペアの組み合わせ数が多くなりすぎて、計算時間がかかりすぎるので、効率が悪い。そこで、超水滴衝突併合演算手段７ｃは、ペアの組み合わせ、つまり、一方の超水滴が他方の超水滴にそれぞれ衝突することが一通りに定まったことを、“［Ｎ／２］”個のペアで代表させることとしている。

そうすることで、超水滴衝突併合演算手段７ｃは、超水滴同士の衝突併合にかかる計算時間の短縮を図り、単位時間当たりの計算コストがＮ^２に比例するのではなく、Ｎに比例するようにしている。

超水滴衝突併合演算手段７ｃでは、ペアの組み合わせ数を減らした分、衝突併合確率Ｐ_ｊｋ ^（ｓ）を上昇させることで、分割空間内の総衝突数の期待値が再現できるように調整している。

つまり、超水滴衝突併合演算手段７ｃでは、超水滴のｉ番目のペア（ｊ，ｋ）に対し、衝突併合確率はｐｉ：＝（Ｎ−１）Ｐ_ｊｋ ^（ｓ）＝（Ｎ−１）ｍａｘ（ｎ_ｊ，ｎ_ｋ）Ｐ_ｊｋとしている。そうすることで、実際に分割空間内の総衝突併合数の期待値が次に示す数式（８）のように再現されることとなる。

そして、超水滴衝突併合演算手段７ｃは、全ての衝突併合候補対ｉ＝１，・・・，“［Ｎ／２］”について以下の演算を行う。乱数をＲａｎ（０，１）とし、超水滴のペア（ｊ，ｋ）が時刻ｔからｔ＋Δｔまでに実際に衝突する回数をｑとすると、次に示す数式（９）のようになる。

なお、この数式（９）において、ｑは０（衝突併合しない）または１（衝突併合する）のどちらかになるべきであるが、数値計算の効率も考えて複数回衝突する場合も考慮している。このｑ＝０の場合、ｉ番目の超水滴のペアには、何の演算も実行せず、ｑ≠０の場合、次に示す数式（１０）のような演算を実行する。なお、この数式（１０）が超水滴におけるモンテカルロ法を示したものである。

なお、この数式（１０）において、多重度ｎ_ｊ≦多重度ｎ_ｋの場合、ｉの操作において、ｊとｋと逆転させる。
また、数式（１０）において“´”の付いた量は、衝突併合後に更新された値である。そして、超水滴衝突併合演算手段７ｃでは、数式（８）および数式（９）に示した演算を全ての分割空間について行って、多重度が０になったら超水滴は除外していく。これによって、時刻ｔにおける衝突併合の演算は終了し、時刻ｔ＋Δｔの状態が得られることとなる。
以上説明したように、超水滴衝突併合演算手段７ｃは、超水滴同士の衝突併合について演算を実行している。

流体力学モデル演算手段９は、入力手段３で入力された周辺環境データおよび微物理モデル演算手段７で演算された結果に基づき、雲の流体力学過程をモデル化した流体力学モデル（非静力学モデル）を用いた演算を行うものである。そして、この流体力学モデル演算手段９では、演算結果を微物理モデル演算手段７にフィードバックしている。
この流体力学モデル演算手段９では、次に示す数式（１１）から数式（１５）までを用いて、雲の流体力学過程の演算を行っている。

これらの数式（１１）から数式（１５）において、Ｄ／Ｄｔ＝∂／∂ｔ＋「ｖ」・∇はラグランジュ微分を、ρ＝ρ_ｄ＋ρ_ｖは湿潤大気密度を、ｑ_ｖ＝ρ_ｖ／ρは水蒸気の混合比を、「ｖ」は風速を、Ｔは気温を、θは温位を、Π＝（Ｐ／Ｐ_０）^{（Ｒｄ／ｃｐ）}はエクスナー関数を、Ｐ_０は基準気圧１０００［Ｐａ］を、ρ_ωは空間単位体積当たり液体の水の質量を、ρＳ_ｖは単位時間当たりに水になった水蒸気の質量を、「ｇ」は重力定数を、λとκは乱流による輸送係数を、Ｒ_ｄは乾燥空気の気体定数を、ｃ_ｐは定圧比熱を、Ｌは水蒸気の潜熱を、それぞれ示している。

これらの中で、ρ_ωとＳ_ｖとは、雲の微物理過程、すなわち、微物理モデル演算手段７による演算によって決定される量である。数式（１１）のρ_ω「ｇ」の項は、水滴により空気が引きずられる効果を表している。

また、数式（１３）のＳ_ｖＬ／ｃ_ｐΠの項は、水蒸気が液体の水になる時に潜熱を放出して大気が温められる効果を表している。さらに、数式（１５）のＳ_ｖの項は、水蒸気が水に凝結したり、水滴が蒸発したりすることにより、大気中の水蒸気の量が増減する元となることを表している。

そして、流体力学モデル演算手段９は、演算結果を微物理モデル演算手段７にフィードバックさせる。ここで、流体力学モデル演算手段９から微物理モデル演算手段７に、また、微物理モデル演算手段７から流体力学モデル演算手段９に送られる情報について明確にしておく。

流体力学モデル演算手段９から微物理モデル演算手段７には、風速、気圧、気温、湿度の情報が送られる。微物理モデル演算手段７において、これらの情報を通して、超水滴モデルを用いた演算が実行される。すなわち、超水滴の運動や凝結成長が、当該超水滴の位置における風速、気圧、気温、湿度に依存している（影響を受ける）からである。

また、微物理モデル演算手段７から流体力学モデル演算手段９には、ρ_ωとＳ_ｖとが送られる。すなわち、流体力学モデルは、これらρ_ωとＳ_ｖという項を通して、超水滴モデルの影響を受ける。具体的にこれらの量は、分割空間ごとに超水滴を分別し、各分割空間内において和をとって計算されている。つまり、ある分割空間におけるρ_ωとＳ_ｖとは、次に示す数式（１６）および数式（１７）に示すように求められる。

なお、これら数式（１６）および数式（１７）において、ΔＶは分割空間の体積であり、和は分割空間内の全ての超水滴に対してとることとしている。また、ｍ_ｉ：＝ｎ_ｉ（４π／３）Ｒ_ｉ ^３ρ_ｌｉｑは、各超水滴が表現している水の質量である。ちなみに、ρ_ｌｉｑ＝１［ｇ／ｃｃ］は水の密度を表している。

出力手段１１は、入力手段３に入力された計算時間に関するデータ、時間ステップ（刻み時刻）と、総計算時間とに従って、微物理モデル演算手段７および流体力学モデル演算手段９で演算された結果を出力するものである。

このシミュレーション装置１によれば、微物理モデル演算手段７において、分割空間内で同じ属性を有している複数の実水滴の超水滴として取り扱って、超水滴モデルを雲の微物理モデルとして用いることで、実水滴を分布関数として取り扱うことなく、当該実水滴の属性数が増加しても、雲の形成、降雨等の計算時間の増加を抑制することができ、高い精度で予測することができる。

また、このシミュレーション装置１によれば、微物理モデル演算手段７の超水滴運動演算手段７ａ、超水滴凝結成長演算手段７ｂおよび超水滴衝突併合演算手段７ｃによって、超水滴について、超水滴の運動、超水滴の凝結成長および超水滴の衝突併合の演算を行っているので、より高い精度で実水滴の経過時間に伴う変化を予測することができる。

さらに、シミュレーション装置１によれば、超水滴衝突併合演算手段７ｃによる超水滴の衝突併合の演算について、モンテカルロ法を用いることで、計算時間を大幅に短縮することができる。

（シミュレーション装置の動作）
次に、図６に示すフローチャートを参照して、シミュレーション装置１の動作について説明する（適宜、図１参照）。
まず、シミュレーション装置１は、入力手段３に、初期変数が入力され、これにより、当該装置１で演算に必要とされる種々の変数が設定される（ステップＳ１）。続いて、シミュレーション装置１は、入力手段３によって入力された初期変数を、データ記憶手段５に記憶する（ステップＳ２）。

そうすると、シミュレーション装置１は、微物理モデル演算手段７の超水滴衝突併合演算手段７ｃによって、超水滴の衝突併合の演算を実行する（ステップＳ３）。また、シミュレーション装置１は、微物理モデル演算手段７の超水滴運動演算手段７ａによって、超水滴の速度の演算を実行して更新する（ステップＳ４）。さらに、シミュレーション装置１は、微物理モデル演算手段７の超水滴凝結成長演算手段７ｂによって、超水滴の凝結成長の演算を実行する（ステップＳ５）。

そして、シミュレーション装置１は、流体力学モデル演算手段９によって、流体力学モデルを用いた演算を実行し、水蒸気の加減を行う（ステップＳ６）。そして、シミュレーション装置１は、この水蒸気を加減した結果を、微物理モデル演算手段７にフィードバックする。

そうすると、シミュレーション装置１は、微物理モデル演算手段７の超水滴運動演算手段７ａによって、超水滴の移動距離の演算を実行する（ステップＳ７）。そして、シミュレーション装置１は、微物理モデル演算手段７の超水滴凝結成長演算手段７ｂによって、分割空間単位体積当たり液体の水の質量ρ_ωを計算する（ステップＳ８）。そして、シミュレーション装置１は、この分割空間単位体積当たり液体の水の質量ρ_ωを計算した結果を、流体力学モデル演算手段９にフィードバックする。

そして、シミュレーション装置１は、フィードバックされた結果に従って、流体力学モデル演算手段９によって新たな流体力学モデルの計算を実行する（ステップＳ９）。そして、シミュレーション装置１は、現在の時刻ｔ（右辺のｔｉｍｅ）に時間ステップΔｔを加算して経過時間（左辺のｔｉｍｅ）を求め（ステップＳ１０）、この経過時間が総計算時間（ｔｉｍｅ＿ｍａｘ）に達したか否かを判定する（ステップＳ１１）。

この後、シミュレーション装置１は、経過時間が総計算時間（ｔｉｍｅ＿ｍａｘ）に達したと判定しない場合（ステップＳ１１、Ｎｏ）には、ステップＳ２に戻り、動作を継続し、達したと判定した場合（ステップＳ１１、Ｙｅｓ）には、動作を終了する。

（シミュレーション装置による具体的な計算例、シミュレーション結果）
次に、シミュレーション装置１による具体的な計算例と、シミュレーション結果について説明する。
まず、シミュレーション装置１に入力される初期変数を、以下のように設定する。初期時刻ｔｉｍｅ＝０［ｓｅｃ］、時間ステップΔｔ＝０．２［ｓｅｃ］、総計算時間ｍａｘ＿ｔｉｍｅ＝５４００．０［ｓｅｃ］、分割空間の高さΔｚ＝４０［ｍ］、分割空間の幅Δｘ＝５０［ｍ］、全体空間の高さｍａｘ＿ｚ＝１０［ｋｍ］、全体空間の幅ｍａｘ＿ｘ＝１２［ｋｍ］、垂直方向の分割空間の個数ｉｚ＿ｍａｘ＝ｍａｘ＿ｚ／Δｚ＝２５０［個］、水平方向の分割空間の個数ｉｘ＿ｍａｘ＝ｍａｘ＿ｘ／Δｘ＝２４０［個］、超水滴の総個数Ｎ＝２．４×１０^７［個］を初期変数として設定する。

また、全ての分割空間に、風速、気圧、温度、湿度等（温位θ、大気密度ρ、水蒸気混合比ｑ_ｖ、圧力Ｐ、温度Ｔ、相対湿度Ｓ）（以下、総括する場合には、流体変数という）を初期変数として設定している。なお、これらの初期変数は、実際に観測した観測値に基づいて設定している。

なお、ここでは、全体空間および分割空間の奥行き方向の構造は一様としている。つまり、全体空間および分割空間が奥行き方向に同じような状態で連なっているとして、奥行きを無視した空間（平面な領域）に対してシミュレーションを行って、当該シミュレーションを簡略化し、計算時間の短縮を図っている。

そして、各分割空間での風速、気圧、温度、湿度を設定し、シミュレーション装置１に対して、全体空間の中のある分割空間（例えば、上空５ｋｍ以下）に不安定で湿った空気があるという初期条件を与える。さらに、この初期条件を与える分割空間の一部分（例えば、数百メートル四方）の温度を周囲の分割空間よりも１度上げることで、不安定化を促進する。

そして、シミュレーション装置１は、全ての超水滴の変数を初期化し、超水滴のｘ座標を保存する記憶域ＰＴＬ＿ＸＸに記憶されている座標（以下、超水滴のｘ座標ＰＴＬ＿ＸＸという）と超水滴のｚ座標を保存する記憶域ＰＴＬ＿ＺＺに記憶されている座標（以下、超水滴のｚ座標ＰＴＬ＿ＺＺという）とに基づいて、各超水滴を全体空間内にランダムに配置する。

それから、実際の大気中（全体空間）には、１．０×１０^８［個／ｍ^３］個の凝結核が飛んでいるとする（正確には、大気中に多数の凝結核が飛んでいて、これらそれぞれに水蒸気が凝結することで水滴ができる。よって、ここでいう凝結核は、半径がほとんど０の水滴のことを指す）と、シミュレーション装置１は、超水滴の多重度の初期値を全て同じ値として、データ記憶手段５における超水滴の多重度を保存する記憶域ＰＴＬ＿ＮＮに１．０×１０^８［個／ｍ^３］×１０［ｋｍ］×１２［ｋｍ］／Ｎ＝１．２×１０^１６／２．４×１０^７＝５．０×１０^８［個］を記憶する。なお、ここで、任意の多重度を持つ超水滴が、任意属性の属性値に関するε程度の幅を持って、実水滴の分布関数により再現されるものであるということを言及するのであれば、シミュレーション結果によって、超水滴の総個数および多重度が適当であるか否かを確認する必要がある。しかし、ここでは、超水滴の総個数および多重度が適当である（精度が十分である）という前提で以下の説明をする。

また、シミュレーション装置１は、各超水滴の持つ凝結核の質量を０〜１．０×１０^−１６［ｇ］の一様ランダムな値を持つものとして取り扱っているので、データ記憶手段５における超水滴の凝結核の質量を保存する記憶域ＰＴＬ＿ＣＣＮ（凝結核：Cloud Condensation Nucleiの略）（ｉＮ）（超水滴の凝結核の質量ＰＴＬ＿ＣＣＮ（ｉＮ）という、以下同様）に１．０×１０^−１６×ＲＡＮＤＯＭ［ｇ］を記憶する。なお、ＲＡＮＤＯＭは［０，１］の一様な乱数である。

このように、各超水滴の持つ凝結核の質量と、周辺の湿度が決定すると、凝結核が吸収する水蒸気の量が定まり、シミュレーション装置１では、超水滴の半径が決定されることとなる。つまり、超水滴の凝結核の質量ＰＴＬ＿ＣＣＮ（ｉＮ）と、超水滴のｘ座標ＰＴＬ＿ＸＸおよび超水滴のｚ座標ＰＴＬ＿ＺＺにおける湿度ＦＬＤ＿Ｓ（ｉｘ，ｉｚ）によって、超水滴の半径ＰＴＬ＿ＲＲ（ｉＮ）が決定される。

この初期状態の段階では、どの分割空間においても湿度は飽和しておらず、凝結核に凝結している水蒸気の量は極微小である。それゆえ、水滴の半径は約１０^−８［ｍ］程度の極めて小さなものであり、視認できないこと（つまり、大気は透明に見えること）を示している。また、各超水滴の半径が定まるので、終端速度も決定されることになり、超水滴の終端速度ＰＴＬ＿ＶＺ（ｉＮ）とする。この超水滴の終端速度ＰＴＬ＿ＶＺ（ｉＮ）は、ほとんど０［ｍ／ｓ］と等しくなり、超水滴は風に流されて動くだけで落ちてこないことを意味している。

そして、シミュレーション装置１では、現在時刻ｔにおける全ての超水滴について演算した結果および現在時刻ｔにおける流体力学モデルについて演算した結果を、データ記憶手段５に記憶する。

その後、シミュレーション装置１は、各分割空間に、どんな超水滴があるのかを分別する。分割空間が全部で（ｉｘ＿ｍａｘ）×（ｉｚ＿ｍａｘ）個ある。例えば、ｉｘ＝２、ｉｚ＝３の分割空間には、ｉｎ＝３，８，１００，５１１，１２３４の超水滴のあることが複数の超水滴を分別した結果わかったとする。

そして、シミュレーション装置１は、すべての分割空間（ｉｘ，ｉｚ）、ｉｘ＝１，・・・，ｉｘ＿ｍａｘ、ｉｚ＝１，・・・，ｉｚ＿ｍａｘに対して、以下の処理を実行する。

まず、シミュレーション装置１は、計算対象となる分割空間にＮｓ個の超水滴があるとし、ランダムに［Ｎｓ／２］個のペアを生成する。例えば、（ｉｘ，ｉｚ）＝（２，３）の分割空間では、Ｎｓ＝５なので、［Ｎｓ／２］＝２となり、これに従って、２つペアをランダムに生成したものを（３，１００）、（８，１２３４）とする。正確には、ランダムに（３，８，１００，５１１，１２３４）を並べ代えて（３，１００，８，５１１，１２３４）とし、この順で２つずつペアにしていく。

そして、各ペアの衝突確率が超水滴の半径で決定することになる。ここで、ｉ番目のペアの衝突確率をｐ_ｉとする。この時に、シミュレーション装置１は、微物理モデル演算手段７の超水滴衝突併合演算手段７ｃによって、［０，１］の乱数ＲＡＮＤＯＭを発生させて、ＲＡＮＤＯＭ＜ｐ_ｉの場合、ｉ番目のペアは衝突させて、ＲＡＮＤＯＭ＞ｐ_ｉの場合、ｉ番目のペアは衝突させないこととしている。

なお、時間ステップΔｔが十分に小さく無いと、衝突確率ｐ_ｉが１よりも大きくなってしまう場合があり、この衝突確率ｐ_ｉという確率が１よりも大きくなってしまうのは、通常、理論の適用範囲外となる。しかし、衝突確率ｐ_ｉが１よりも大きくなってしまってもシミュレーション結果には大きな影響を与えないことが予測される（時間ステップΔｔがある程度小さければ、衝突確率ｐ_ｉが１より大きくなることは滅多に生じないので）。そこで、シミュレーション装置１では、ｐ_ｉ＞１の場合も想定して、［０，１］の乱数ＲＡＮＤＯＭを衝突併合の判定をする度に発生させて、ＲＡＮＤＯＭ＜ｐ_ｉ−［ｐ_ｉ］の場合、ｉ番目のペアは［ｐ_ｉ］＋１回衝突させて、ＲＡＮＤＯＭ＞ｐ_ｉ−［ｐ_ｉ］の場合、ｉ番目のペアは［ｐ_ｉ］回衝突させる。つまり、ｐ_ｉ＜１の場合、すなわち、［ｐ_ｉ］＝０の場合、このペアｉの衝突併合がｐ_ｉという確率で実現されることがわかる。

例えば、（ｉｘ，ｉｚ）＝（２，３）の分割空間の１番目のペア（３，１００）の衝突確率が０．４（４０％の確率で衝突する）だとすると、ＲＡＮＤＯＭが０．４より小さければ、超水滴［３］と超水滴［１００］とを衝突させ、衝突後の属性を更新する。逆にＲＡＮＤＯＭが０．４より大きければ、超水滴［３］と超水滴［１００］とを衝突させない。例えば、衝突確率が１．３の場合（１３０％の確率で衝突する）、１．３−［１．３］＝１．３−１＝０．３である。この場合に、ＲＡＮＤＯＭが０．３より小さければ、１＋１＝２回衝突させ、ＲＡＮＤＯＭが０．３より大きければ、１回だけ衝突させる。

そして、超水滴の終端速度は超水滴の半径の大きさにより決定されるので、シミュレーション装置１は、微物理モデル演算手段７の超水滴運動演算手段７ａによって、時刻ｔにおける超水滴の半径ＰＴＬ＿ＲＲ（ｉＮ）により、超水滴の終端速度ＰＴＬ＿ＶＺ（ｉＮ）を計算し、これをｔｉｍｅ＝ｔｉｍｅ＋Δｔの終端速度とする。

また、超水滴は、当該超水滴の周辺の大気の湿度や温度に応じて、水蒸気を吸収したり、水蒸気が蒸発したりすることで、当該超水滴の半径の大きさが変わるため、シミュレーション装置１は、微物理モデル演算手段７の超水滴凝結成長演算手段７ｂによって、超水滴の凝結成長を計算する。ここで、超水滴の半径ＰＴＬ＿ＲＲ（ｉＮ）を時刻ｔｉｍｅの状態から時刻ｔｉｍｅ＋Δｔの状態に更新する。

そして、シミュレーション装置１は、微物理モデル演算手段７の超水滴凝結成長演算手段７ｂにより、超水滴がどれくらい水蒸気を吸収或いは放出したかがわかるので、全体空間（ｉｘ，ｉｚ）における水蒸気量を加減する。つまり、シミュレーション装置１は、流体力学モデル演算手段９における流体力学モデルの水蒸気混合比ＦＬＤ＿ＱＶ（ｉｘ，ｉｚ）の値を加減する。

そして、超水滴は風速の影響を受けて移動し、超水滴の速度は、周辺の風速と当該超水滴の終端速度との差によって決定するので、シミュレーション装置１は、超水滴のｘ座標ＰＴＬ＿ＸＸ（ｉＮ）と超水滴のｚ座標ＰＴＬ＿ＺＺ（ｉＮ）との時刻ｔｉｍｅ＋Δｔにおける値を、超水滴の終端速度ＰＴＬ＿ＶＺ（ｉＮ）と、流体力学モデルにおける風速のｘ成分ＦＬＤ＿ＵＵと、流体力学モデルにおける風速のｚ成分ＦＬＤ＿ＷＷとに基づいて更新する。

そして、シミュレーション装置１は、空間単位体積当たり液体の水の質量ρ_ωを求める。この空間単位体積当たり液体の水の質量ρ_ωは、各分割空間にある超水滴の質量を集計すれば求められる。例えば、分割空間にｉＮ＝１，３，１０，１１の超水滴が存在しているとすると、ρ_ω＝Σａ×ＰＴＬ＿ＮＮ（ｉＮ）×４π／３｛ＰＴＬ＿ＲＲ（ｉＮ）｝^３（ｉＮ＝１，３，１０，１１、ａは水の比重）によって求めることができる。

その後、シミュレーション装置１は、流体力学モデル演算手段９において、流体変数を更新する。この流体変数は、一般的な非静力学方程式に従って経過時間に伴って変化する。全ての分割空間における流体力学モデルにおける変数ＦＬＤ＿任意（ｉｘ、ｉｚ）を、非静力学方程式に従って、時刻ｔｉｍｅ＋Δｔの状態を更新する。そして、シミュレーション装置１は、総計算時間に達するまで、これらの処理を繰り返し、実行する。

（シミュレーション装置による演算結果）
次に、シミュレーション装置１による演算結果（シミュレーション結果）について、図７から図１５を参照して説明する（適宜、図１参照）。
図７から図１５までは、横軸に全体空間の幅［ｍ］をとり、縦軸に全体空間の高さ［ｍ］をとり、各時刻における雲の生成状態を示したものである。

図７は初期時刻０［ｓｅｃ］から６００［ｓｅｃ］後、図８は初期時刻０［ｓｅｃ］から１２００［ｓｅｃ］後、図９は初期時刻０［ｓｅｃ］から１８００［ｓｅｃ］後、図１０は初期時刻０［ｓｅｃ］から２４００［ｓｅｃ］後、図１１は初期時刻０［ｓｅｃ］から３０００［ｓｅｃ］後、図１２は初期時刻０［ｓｅｃ］から３６００［ｓｅｃ］後、図１３は初期時刻０［ｓｅｃ］から４２００［ｓｅｃ］後、図１４は初期時刻０［ｓｅｃ］から４８００［ｓｅｃ］後、図１５は初期時刻０［ｓｅｃ］から５４００［ｓｅｃ］後を示している。

これら図７から図１５までに示したように、何も存在していない空間から、目に見えない水滴が衝突併合或いは凝結成長して、雲を形成し、ある時刻（約１８００ｓｅｃ）から雨が降り始め、ある時刻（約５４００ｓｅｃ）まで雨が降り続くことがわかる。

（超水滴モデルから超粒子モデルに拡張する場合について）
これまでのシミュレーション装置１の説明は、実水滴の経過時間における変化を、超水滴モデルを用いてシミュレーションすることについて述べてきた。これより、シミュレーションの対象を実水滴から実粒子に拡張する場合、すなわち、超水滴モデルを用いた方法から超粒子モデルを用いた方法（超粒子法）に拡張する場合について説明する。

この超粒子法は、流体中で多数の実粒子が衝突併合を繰り返す場合に、時間経過による変化を数値シミュレーションする場合に適用できる。なお、一般に数値シミュレーションの対象となる「系」が以下に述べる要件を満たせば、超粒子法を適用することができる。

この「系」には、時間ｔと空間「ｒ」とによる“時空間”の概念がある。この空間「ｒ」の次元をｄｒとし、この次元ｄ_ｒを１以上の任意の整数で表すこととする。

まず、この空間「ｒ」にＮ個の実粒子が存在するとし、Ｎは時間ｔの関数であり、Ｎ＝Ｎ（ｔ）とする。そして、個々の実粒子はｍ_ｐ個の属性を有し、これを「Ａ」＝（Ａ^（１），Ａ^（２），・・・，Ａ^（ｍｐ））と表すこととする。そして、実粒子は空間内を運動するので、少なくとも属性には、実粒子の速度「ｖ」が含まれていて、１番目からｄ_ｒ番目の属性は速度を表すこととし、「ｖ」＝（Ａ^（１），Ａ^（２），・・・，Ａ^（ｄｒ））となる。ただし、超水滴の場合と同様に、「ｖ」が他の属性の関数として一意に決まる場合は顕わに属性として扱わない場合もある。

また、個々の実粒子は空間上の位置座標「ｑ」を持っており、この位置座標「ｑ」＝（ｑ^（１），ｑ^（２），・・・，ｑ^（ｄｒ））と表すこととする。そうすると、全ての実粒子は、位置と属性は（「ｑ」_ｉ，「Ａ」_ｉ）（ｉ＝１，２，・・・，Ｎ）により定まることとなる。また、各実粒子の位置と属性とは経過時間と共に変化するので、これらは時間の関数であり、「ｑ」_ｉ＝「ｑ」_ｉ（ｔ）、「Ａ」_ｉ＝「Ａ」_ｉ（ｔ）、ｉ＝１，２，・・・，Ｎと表すことができる。

そして、空間「ｒ」は流体で満たされているとし、この流体はｍ_ｆ個の流体場変数「ａ」で特徴付けられており、この流体場変数「ａ」を「ａ」＝（「ａ」^（１），「ａ」^（２），・・・，「ａ」^（ｍｆ））とする。この流体場変数「ａ」は時間と空間の関数となるので、「ａ」＝「ａ」（「ｒ」，ｔ）と表すことができる。

また、実粒子の属性「Ａ」には、次に示す数式（１８）のような時間発展方程式が与えられているとする。なお、この数式（１８）が実粒子および超粒子の属性時間発展方程式に該当する。

ただし、この数式（１８）において、「Ｆ」（「Ａ」，「ａ」（「ｑ」，ｔ））＝（Ｆ^（１）（「Ａ」，「ａ」（「ｑ」，ｔ）），Ｆ^（２）（「Ａ」，「ａ」（「ｑ」，ｔ）），・・・，Ｆ^（ｍｐ）（「Ａ」，「ａ」（「ｑ」，ｔ）））である。この「Ｆ」は、一般に実粒子の存在する位置「ｑ」における流体場「ａ」（「ｑ」，ｔ）にも依存している。また、この時間発展方程式群は、単体の実粒子が流体場中でどの様に振る舞うか（どの様な速度で移動する）を規定している。それゆえ、実粒子の位置「ｑ」は、当該速度に従って移動し、次に示す数式（１９）のようになる。なお、この数式（１９）が実粒子および超粒子の位置座標時間発展方程式に該当している。

そしてまた、流体場「ａ」の時間発展方程式群が、微積分方程式の形で与えられていて、この微積分方程式は次に示す数式（２０）のようになる。なお、この数式（２０）が実粒子の流体場時間発展方程式に該当している。

この数式（２０）において、「Ｓ」（「ｒ」，ｔ）＝（Ｓ^（１），Ｓ^（２），・・・，Ｓ^（ｍｆ））は、実粒子から流体への相互作用を記述する項であり、すべての実粒子の位置と属性とにより決定され、「Ｓ」（「ｒ」，ｔ）＝「Ｓ」（「Ａ」_１，「Ａ」_２，・・・，「Ａ」_Ｎ，「ｑ」_１，「ｑ」_２，・・・，「ｑ」_Ｎ，「ｒ」，ｔ）＝「Ｓ」（｛「Ａ」_ｉ，「ｑ」_ｉ｝，「ｒ」，ｔ）となる。

そして、実粒子同士は、確率的に衝突併合する。つまり、全体空間中の、十分小さな体積ΔＶの分割空間内に存在する２個の実粒子ｊとｋとは、時刻ｔより十分短い時間Δｔの間に、次に示す数式（２１）で与えられる確率により、衝突併合する。

この数式（２１）において、Ｃ（「Ａ」_ｊ，「Ａ」_ｋ，「ａ」（「ｒ^〜」，ｔ））は実効衝突併合断面積であり、２個の実粒子の属性「Ａ」_ｊと「Ａ」_ｋ、分割空間における流体場「ａ」（「ｒ^〜」，ｔ）の関数である。ここで、「ｒ^〜」は、ΔＶを十分に小さな体積としているため、分割空間内のどの位置で流体場を評価したとしても、誤差は高次の微少量となり、当該分割空間内のどの点で取得してもよい。「ｒ^〜」：＝（「ｑ」_ｊ＋「ｑ」_ｋ）／２と表すことも可能である。

そして、実粒子同士の衝突併合により、２個の実粒子は新たな１個の実粒子となることを、「実粒子が衝突併合する」と定義する。

この新たな実粒子も属性「Ａ」で特徴付けられる実粒子である。また、実粒子同士が衝突併合した結果、どの様な属性「Ａ」´の新たな実粒子が生成されるのかということも一般的には確率的に決定される。つまり、属性Ａ_ｊとＡ_ｋとを持つ実粒子が衝突併合するという条件の元、属性「Ａ」´の実粒子が生成される確率分布ｐ_ｃ（「Ａ」_ｊ，「Ａ」_ｋ；「Ａ」´）が与えられているとする。

なお、衝突併合により新しく生成された実粒子は、衝突が起こった分割空間内のΔＶのどこかに存在するものとして、当該実粒子の位置「ｑ」´は確率的に、或いは、決定論的に決定されるとする。一般には、実粒子（「Ａ」_ｊ，「ｑ」_ｊ）と実粒子（「Ａ」_ｋ，「ｑ」_ｋ）とが衝突併合し、属性「Ａ」´の実粒子が生成されたという条件の下、新しく生成された実粒子の位置が「ｑ」´であるという確率分布は、θ_ｃ（「Ａ」_ｊ，「Ａ」_ｋ，「ｑ」_ｊ，「ｑ」_ｋ，「Ａ」´，ΔＶ；「ｑ」´）に従うものとする。

ここで、実粒子同士が衝突した場合に、分裂する過程（衝突分裂）も想定することができる。この場合、全体空間の中の十分小さな体積ΔＶの分割空間内に存在する２個の実粒子ｊとｋとが時刻ｔより十分短い時間Δｔの間に衝突分裂する確率は、次に示す数式（２２）のように与えられる。

ただし、この数式（２２）において、Ｓ（「Ａ」_ｊ，「Ａ」_ｋ，「ａ」（「ｒ^〜」，ｔ））は実効衝突分裂断面積である。なお、衝突分裂とは、２個の実粒子が衝突後、別の属性を持つ２個以上の新たな実粒子となることである。これら新たな実粒子も複数の属性からなる属性群の中のいずれかの属性で特徴付けられる実粒子であるとする。この衝突分裂した結果、どのような属性の実粒子が何個生成されるかということも一般に確率的に決定されることになる。

ここで、属性「Ａ」_ｊを持つ実粒子と属性「Ａ」_ｋを持つ実粒子とが衝突分裂するという条件の元、属性（「Ａ」´^［１］，「Ａ」´^［２］，・・・，「Ａ」´^［ｎ］）の全部でｎ個の実粒子が生成される確率分布ｐ_ｓ（「Ａ」_ｊ，「Ａ」_ｋ；Ａ´^［１］，Ａ´^［２］，・・・，Ａ´^［ｎ］）が与えられているとする。

さらに、実粒子は、単体分裂してもよい。単体分裂とは、実粒子同士が衝突することなく、属性「Ａ」を持つ１個の実粒子から２個以上の新たな実粒子となることであり、新たな超粒子は属性（「Ａ」´^［１］，「Ａ」´^［２］，・・・，「Ａ」´^［ｎ］）を有している。なお、単体分裂の分裂過程は、確率的な或いは決定論的な法則に従うものである。

ここまで「流体中で多数の実粒子が衝突併合を繰り返す系」について一通りの説明を行ってきたが、これを超粒子法に適用する場合について説明する。

ただし、超粒子法において、衝突併合の過程よりも衝突分裂の過程が支配的な場合、当該超粒子法による数値シミュレーションの際に、全超粒子数が経過時間と共に急激に増加していくことが想定され、計算効率が悪くなる可能性がある。
或いは、超粒子法において、単体分裂の分裂過程が、衝突併合の過程や衝突分裂の過程よりも支配的な場合、当該超粒子法による数値シミュレーションの際に、全超粒子数が経過時間と共に急激に増加していくことが想定され、計算効率が悪くなる可能性がある。

ここで、図１に示したシミュレーション装置１を、超粒子モデルに適用した場合について、図１６を参照して説明する。この図１６に示したように、シミュレーション装置１Ａは、これまで説明した超粒子および超粒子モデル（超粒子法）を用いたシミュレーションを実行するもので、入力手段（変数入力手段）３ａと、データ記憶手段５ａと、演算手段１３と、出力手段１１ａとを備えている。なお、以下、シミュレーション装置１Ａの説明は、ここまでの超粒子法の説明と重複する箇所があるが、これよりは、当該超粒子法を、当該装置１Ａとして実現した場合について説明するものである。

シミュレーション装置１Ａの説明に先立ち、改めて、超粒子法について、これまで説明した「系」との関係に触れつつ説明する。超粒子法とは、「系」が有している膨大な自由度のうち、不必要に細かい部分を捨象する（或いは、粗視化する、概略化する）方法の一種と位置付けることができる。また、各実粒子は、流体中に離散する多数の離散粒子であると言え、これら「多数の離散粒子が流体中で衝突併合を繰り返す系」を「超粒子化」して（超粒子化された系として）、シミュレーション装置１Ａによって、シミュレーションする場合について説明する。

まず、超粒子化された系における時空間の定義は、元の「多数の離散粒子が流体中で衝突併合を繰り返す系」と同じである。この時空間にＮ_ｓ個の超粒子が存在しているとする。個々の超粒子は、実粒子と同様に空間上の位置と、ｍ_ｐ個の属性（「ｑ」_ｉ，「Ａ」_ｉ）を有している。そして、個々の超粒子は、多重度ｎ、ｎ＝１，２，３，・・・を有している。この多重度ｎは超粒子がｎ個の実粒子を含んでいるということを表現した量である。そうすると、全ての超粒子の状態は、位置と属性と多重度とで表され、（「ｑ」_ｉ，「Ａ」_ｉ，ｎ_ｉ）、ｎ＝１，２，３，・・・，Ｎ_ｓにより決定される。そして、これらの超粒子は次に示す数式（２３）のように表現される。

そして、超粒子の属性「Ａ」と位置「ｑ」は、実粒子と同様に数式（１８）と数式（１９）とに従い時間発展する。流体場「ａ」も元の系と同様の時間発展方程式に従うことになるが、超粒子から流体への相互作用の項は、超粒子の表現する実粒子によって評価され、次に示す数式（２４）のように表現される。なお、この数式（２４）が超粒子の流体場時間発展方程式に該当する。

この数式（２４）において、「Ｓ^〜」（「ｒ」，ｔ）は超粒子から流体への相互作用を記述する項（相互作用量）である。この項は、Ｎ_ｓ個の全超粒子が表現するＮ個の全粒子を｛「Ａ」_ｉ，「ｑ」_ｉ｝_ｒと記述して、元の系の相互作用項「Ｓ」によって、次に示す数式（２５）のように定義される。

また、超粒子同士の衝突併合は次のように定義する。まず、ｎ_ｊ個とｎ_ｋ（＞ｎ_ｊ）個の多重度を持つ超粒子が衝突併合する状況を想定すると、多重度ｎ_ｋ−ｎ_ｊの属性「Ａ」_ｋの超粒子と多重度ｎ_ｊの属性「Ａ」´を持つ超粒子が生成されるとする。次に、ｎ_ｋ＝ｎ_ｊの場合を想定すると、この場合、多重度“「ｎ_ｊ／２」”と多重度ｎ_ｊ−“「ｎ_ｊ／２」”の両方とも「Ａ」´を持つ超粒子が生成されるとする。よって、いずれの場合も「Ａ」_ｊと「Ａ」_ｋとを持つ実粒子がｍｉｎ（ｎ_ｊ，ｎ_ｋ）組だけ衝突併合した結果、属性「Ａ」´の実粒子がｍｉｎ（ｎ_ｊ，ｎ_ｋ）個生成されたということを意味する。

このことを前提とすると、超粒子同士の衝突併合確率は次に示す数式（２６）のように与えられる。

この数式（２６）に示した衝突併合確率で、全体空間中の流体場「ａ」の空間スケールより十分小さな体積ΔＶの空間内に存在する２個の超粒子ｊと超粒子ｋとは、時刻ｔより十分短い時間Δｔの間に衝突併合する。

この衝突併合した結果、新しく生成された超粒子の属性「Ａ」´は一般には確率的に決定される量であり、元の系と同様に確率分布ｐ_ｃ（「Ａ」_ｊ，「Ａ」_ｋ；「Ａ」´）に従うと定義する。この定義により、実粒子の衝突併合数の期待値を再現することとなる。今、ｎ_ｊ個とｎ_ｋ個の実粒子があるとすると、全部でｎ_ｊｎ_ｋ個のペアが作成できる。それぞれのペアはＰ_ｃの確率で衝突併合するので、衝突併合数の期待値はｎ_ｊｎ_ｋＰ_ｃとなる。

一方、超粒子の衝突併合は、ｍｉｎ（ｎ_ｊ，ｎ_ｋ）個の実粒子の衝突併合を表すため、当該超粒子の衝突併合を表現する実粒子の衝突併合の期待値は、ｍｉｎ（ｎ_ｊ，ｎ_ｋ）Ｐ_ｃ ^（ｓ）＝ｍｉｎ（ｎ_ｊ，ｎ_ｋ）ｍａｘ（ｎ_ｊ，ｎ_ｋ）Ｐ_ｃ＝ｎ_ｊｎ_ｋＰ_ｃである。よって、超粒子の衝突併合により、元の系の期待値が再現されることが確認される。

これによって、「超粒子化された系」が定義された。この系は元の系を近似的に表現したものである。ここで、全超粒子数Ｎｓを十分に多くとれば、つまり、多重度を十分に小さくとれば、シミュレーションしたい現象が最適な近似で再現されると期待できる。つまり、全ての超粒子の多重度を１にすると、「超粒子化された系」は元の系と完全に一致する。

なお、「超粒子化された系」では、元の系である実粒子の個数が衝突併合過程により減少するのに対し、超粒子の個数がほとんど減少しないという特徴を有している。この衝突併合過程により超粒子の個数がほとんど減少しない代わりに、多重度が減少するという特徴を有している。特別な場合として、ｎ_ｊ＝ｎ_ｋ＝１の多重度を持つ超粒子同士が衝突併合すると、多重度０の超粒子が生成され、超粒子の個数も減少する。

全超粒子の個数が大きく変動しないということは、「超粒子化された系」と「元の系」とが常に変わらずよく近似しているということであり、この意味において、超粒子法は、衝突併合を繰り返す系の数値計算モデルとして適していることとなる。なお、超粒子の個数を多くとれば、よりよい近似であるといえる。

そして、超粒子の衝突併合以外の時間発展（経過時間に伴う変化）は、数式（１８）および数式（１９）に基づいて、流体場「ａ」（「ｒ」，ｔ）の時間発展は数式（２４）に基づいて計算する。超粒子の衝突併合過程の計算は、以下に説明するモンテカルロ法によって行うこととする。

まず、超粒子が存在する空間を流体場の空間スケールよりも十分小さな体積を持つ格子（超水滴法における分割空間に相当）に分割する。すると、各格子の内部の超粒子同士は、数式（２６）で与えられる確率Ｐ_ｃ ^（ｓ）に従って衝突併合することとなる。このことを利用して、超粒子の衝突併合過程をモンテカルロ的に数値シミュレーションする。つまり、同一格子内の超粒子を、実際に乱数を発生させて衝突併合させる。時刻ｔにおける衝突併合操作を以下に記載するように行う。

まず、時刻ｔ（例えば、初期時刻ｔ＝０）における各超粒子の位置を特定し、格子ごとに分類しておく。例えば、ある格子の内部に、超粒子がＮ_ｓ ⁻個あるとする。各超粒子ｌ＝１，２，・・・，Ｎ_ｓ ⁻個はそれぞれｎ_ｌ個の多重度を持つとする。この全ての超粒子について、２つの対（ｊとｋ）を、乱数を用いて選択する。なお、この操作の詳細は超水滴の場合と同様である。この対を衝突併合候補対とする。格子内に奇数の超粒子がある場合は、衝突併合候補対に含まれない超粒子が１個残ることとなり、Ｎ_ｓ ⁻／２を越えない最大整数“［Ｎ_ｓ ⁻／２］”個のペアができる。

本来、格子内の超粒子同士は全ての組み合わせについて衝突する可能性があるため、Ｎ_ｓ ⁻（Ｎ_ｓ ⁻−１）／２個の全てのペアの組み合わせについて、衝突併合する可能性を判定する必要がある。

しかし、これでは、ペアの数が多くなりすぎてしまい計算効率が悪く、“［Ｎ_ｓ ⁻／２］”個のペアで代表させて計算コストを抑えることとし、当該計算コストが〜Ｎ_ｓ ^２に応じて増加するのではなく、〜Ｎ_ｓに応じて増加することとなる。衝突併合するペアの数を抑えた代わりに、衝突確率を大きくすることで、格子内の総衝突数の期待値が再現するように調整している。つまり、ｉ番目のペア（ｊ，ｋ）に対し、衝突併合確率ｐ_ｉ：＝（Ｎ_ｓ ⁻−１）Ｐ_ｃ ^（ｓ）（Ａ_ｊ，Ａ_ｋ）＝ｍａｘ（ｎ_ｊ，ｎ_ｋ）（Ｎ_ｓ ⁻−１）Ｐ_ｃ（Ａ_ｊ，Ａ_ｋ）とする。そして、次に示す数式（２７）により、格子内の総衝突併合数の期待値が再現する。

そして、さらに、全ての衝突併合候補対ｉ＝１，・・・，“［Ｎ_ｓ ⁻／２］”について以下の計算を行う。乱数Ｒａｎ∈（０，１）とすると、次に示す数式（２８）のようにγが与えられる。

このγは超粒子のペア（ｊ，ｋ）が時刻ｔからｔ＋Δｔの間に実際に衝突する回数である。本来、γは０（衝突併合しない）か１（衝突併合する）かのどちらかであるべきであるが、数値計算の効率向上を図るため、複数回衝突する場合も想定されている。γ＝０の場合には、ｉ番目のペアには何の操作もせずに、γ≠０の場合には、次に示す数式（２９）のように操作する。なお、この数式（２９）が超粒子におけるモンテカルロ演算を示したものである。

ただし、´の付いた量は、衝突併合後に更新される値である。「Ａ」_η ^〜、η＝０，１，・・・，ｍは、確率分布ｐ_ｃ（「Ａ」_ｊ，「Ａ」_η−１ ^〜；「Ａ」_η ^〜）に従う確率変数であり、「Ａ」_０ ^〜＝「Ａ」_ｋとする。この｛「Ａ」_η ^〜｝＝｛Ａ_ｋ，「Ａ」_η ^〜，・・・，「Ａ」_ｍ−１ ^〜,「Ａ」_ｍ ^〜＝「Ａ」_ｋ´｝を乱数として発生させ、更新される値を確率的に決定する。これにより、超粒子ｊとｍ回衝突併合後の超粒子ｋの属性「Ａ」_ｋ´が決定されることとなる。

同様に、「ｑ」_η ^〜、η＝０，１，・・・，ｍは確率分布θ_ｃ（「Ａ」_ｊ，「Ａ」_η―１ ^〜，「ｑ」_ｊ，「ｑ」_η−１ ^〜，「Ａ」_η ^〜，ΔＶ；「ｑ」_η ^〜）に従う確率変数であり、「ｑ」_０ ^〜＝「ｑ」_ｋとする。この｛「ｑ」_η ^〜｝＝｛「ｑ」_ｋ，「ｑ」_１ ^〜，・・・，「ｑ」_ｍ−１ ^〜，「ｑ」_ｍ ^〜＝「ｑ」_ｋ´｝を、乱数を元に発生させ、更新される値を確率的に決定する。これにより、超粒子ｊとｍ回衝突併合後の超粒子ｋの位置「ｑ」_ｋ´が決定されることとなる。

なお、ｎ_ｊ≦ｎ_ｋの場合、ｊとｋとを逆にして同様の操作を行う。そして、これらの操作を、全ての格子について行う。多重度が０になった超粒子は系から除外する。これによって、時刻ｔにおける衝突併合の操作は終了し、時刻ｔ＋Δｔの状態が得られる。

なお、この操作は、雲形成・降雨現象に適用した超水滴法における微物理モデル演算手段７および流体力学モデル演算手段９における演算を一般化したものである。雲形成・降雨現象では、ｐ_ｃ＝１、つまり１００％の確率で次に示す数式（３０）が実現されるため、乱数を発生させる必要が無く、決定論的に数式（３１）に示すものが得られる。

さらに、超水滴法の場合、数式（２９）において、衝突併合後の超水滴の位置の更新が以下に示すように、更新前の位置を引き継ぐという形で実現されている。
「ｑ」_ｊ´＝「ｑ」_ｊ、「ｑ」_ｋ´＝「ｑ」_ｋ
ちなみに、超水滴の位置「ｒ」とは、超粒子の「ｑ」に該当している。

この衝突併合前と衝突併合後とで超水滴および超粒子の位置が引き継がれるということは、次のように正当化する（保証する）ことができる。
まず、雲形成・降雨現象で適用した超水滴においては、衝突併合の結果生成された超水滴の位置「ｑ」´を決定する条件付き確率分布θ_ｃ（「Ａ」_ｊ，「Ａ」_ｋ，「ｑ」_ｊ，「ｑ」_ｋ，「Ａ」´，ΔＶ；「ｑ」´）の具体的な関数形に重要な物理的意味合いは存在していない。つまり、数値シミュレーションにより、再現される雲形成・降雨現象は、θ_ｃの与え方には大きく因らない（超水滴の前提となる、個々の実水滴の位置のΔＶ程度ずれと雲形成・降雨現象との相関関係はほとんどない）と推測される。

物理現象のモデル化（雲の微物理のモデル化）という観点からすると、流体場の空間変化のスケールよりも十分小さなΔＶの中で実水滴の位置が明確に定まらないという現象の性質が、衝突併合が確率的に起こるという解釈を与える根拠となっている。ちなみに、空間変化のスケールとは、想定している流体場の変化が一定とみなせる様な空間の長さのことをいい、例えば、流体場の気温が一定とみなせる様な長さ、風速が一定とみなせる様な長さのことである。

それゆえ、もし、θ_ｃを変えること、つまり、実水滴の位置をΔＶの中で変えることで、再現される物理現象が大きく変わるならば、実水滴同士が確率的に衝突併合するといそもそもの前提が不適切であったということになる。そこで、このような事柄と数値計算の容易さという利点から衝突併合前と衝突併合後とで超水滴および超粒子の位置が引き継がれるということを前提とした。

これより、シミュレーション装置１Ａの各構成について説明する。なお、超粒子法をシミュレーションに用いる場合の実装方法には任意性があり、ここでは、より汎用性のある実装方法の一例を説明する。
入力手段（変数入力手段）３ａは、実粒子の種別、属性等の実粒子に関するデータと、当該実粒子が存在している流体に関するデータ（周辺環境データ、流体変数）と、計算時間に関するデータ（初期時刻、時間ステップ（刻み時刻）、総計算時間）とを初期変数として入力するものである。

各実粒子は、ｍ個の属性（物理量）を有しており、この属性を「Ａ」＝（Ａ１，Ａ２，Ａ３，・・・，Ａｍ）とする。また、各属性は、流体の影響を受けながら時間変化し、この変化が前記した時間発展方程式に従っている、ここでは、ｄＡ１／ｄｔ＝Ｆ１（Ａ１，Ａ２，Ａ３，・・・，Ａｍ，流体変数）、ｄＡ２／ｄｔ＝Ｆ２（Ａ１，Ａ２，Ａ３，・・・，Ａｍ，流体変数）、ｄＡｍ／ｄｔ＝Ｆｍ（Ａ１，Ａ２，Ａ３，・・・，Ａｍ，流体変数）とし、これらをまとめてｄ「Ａ」／ｄｔ＝「Ｆ」（「Ａ」）、「Ｆ」：＝（Ｆ１，Ｆ２，・・・，Ｆｍ）とする。

そして、同一属性を持つ複数個の実粒子を超粒子として扱い、１個の超粒子に含まれる実粒子の数を多重度ｎで表す。そして、超粒子同士は、何らかの確率で衝突併合、衝突分裂または単体分裂するとする。つまり、衝突併合の場合、超粒子同士がぶつかって結合することが確率的に起こり、生成される新たな超粒子もやはり「Ａ」という属性群で特徴付けられる。

この多重度ｎを少なくしていくと、シミュレーション装置１Ａによるシミュレーションの精度は向上することになる。ちなみに、ｎ＝１にすれば、１個の超粒子は１個の実粒子と同じになる。なお、多重度ｎを小さくすると言うことは、超粒子の個数を多くすることを意味し、計算コストが高くすることになる。このシミュレーション装置１Ａにおいて、ある程度のシミュレーションの精度を維持しつつ、計算コストを削減するために、どのくらい多重度ｎを小さくする、つまり、超粒子の個数を多くとるべきかという基準は、超粒子の分布関数がε程度で再現するようにすることで評価することができる。

ただし、実際にシミュレーション装置１Ａでシミュレーションを実行する際には、シミュレーションしたい現象（例えば、粉塵の振る舞い、マイクロバブルの振る舞い等）が再現可能な程度にｎを十分小さくとればよい。例えば、微小な粉塵が空中で結合して大きな粉塵となり、工場内の床面等に落下する現象であり、シミュレーション装置１Ａでシミュレーションを実行後、再びｎの値を少し変えて実行したとしても、シミュレーションの結果がほとんど変わらないのであれば、そのｎの値で十分であると判断できる。

データ記憶手段５ａは、入力手段３ａで入力された実粒子の種別、属性等の実粒子に関するデータと、当該実粒子が存在している流体に関するデータ（周辺環境データ、流体変数）と、演算手段１３で計算途中であるデータとを記憶し、必要に応じて適宜出力するものである。このデータ記憶手段５ａは、一般的なメモリやハードディスク等によって構成されている。

演算手段１３は、超粒子の各属性の時間発展（経過時間に伴う変化）について演算を行うもので、超粒子運動演算手段１３ａと、超粒子属性更新手段１３ｂと、超粒子衝突併合演算手段１３ｃと、流体場更新手段１３ｄと、超粒子衝突分裂演算手段１３ｅと、超粒子単体分裂演算手段１３ｆとを備えている。

超粒子運動演算手段１３ａは、前記した数式（１８）〜数式（２０）に従った演算を実行し、超粒子の速度について演算するものである。
超粒子属性更新手段１３ｂは、前記した数式（２９）に従った演算を実行し、超粒子の属性について演算するものである。

超粒子衝突併合演算手段１３ｃは、前記した数式（２６）に従った演算を実行し、超粒子の衝突併合について演算するものである。
流体場更新手段１３ｄは、前記した数式（２４）および数式（２５）に従った演算を実行し、超粒子が存在している流体について演算するものである。

超粒子衝突分裂演算手段１３ｅは、前記した数式（２２）に従ったモンテカルロ演算を実行し、超粒子の衝突分裂について演算するものである。
超粒子単体分裂演算手段１３ｆは、超粒子がある確率、或いは決定論に従って２個以上に分裂することについて演算するものである。

なお、シミュレーション装置１Ａでは、シミュレーションの対象となる実粒子が物質的性質・状況において単体分裂を起こさないならば、当初から超粒子単体分裂演算手段１３ｆの機能を停止しておくことも可能である。同様に、シミュレーション装置１Ａでは、超粒子衝突併合演算手段１３ｃおよび超粒子衝突分裂演算手段１３ｅも必要に応じて機能を停止することができる。

この演算手段１３による演算について、一連の流れについて説明する。各属性の時間発展は、ｄ「Ａ」／ｄｔにより行われ、超粒子は自らの速度に従って全体空間内を自由に移動する。この速度を、超粒子運動演算手段１３ａによって演算する。

続いて、超粒子同士の衝突併合は、与えられた衝突併合確率を元に、前記したモンテカルロ法を使って行う。そして、衝突併合確率を与えるには一般に、数式（２１）における衝突併合断面積Ｃ（「Ａ」_ｊ，「Ａ」_ｋ，「ａ」（「ｒ^〜」，ｔ））、或いは超粒子の衝突効率Ｅ_ｊｋのどちらかを与えておけば十分である。なお、数式（３２）において、実効衝突併合断面積Ｃと超粒子の衝突効率Ｅ_ｊｋとの関係を規定している。

そして、超粒子が衝突して分裂する場合（衝突分裂）および超粒子が単体で分裂する場合（単体分裂）について、超粒子衝突分裂演算手段１３ｅおよび超粒子単体分裂演算手段１３ｆによって演算を行う。ここでは、数式（２２）に従うか、何らかの別の決定論的な方程式が与えられていて、これにより分裂の仕方を決定するか、あるいは、分裂確率が与えられていて、それに基づいて、乱数を使ってモンテカルロ的に分裂過程を扱って演算すればよい。

この演算手段１３で演算する対象について例示する。
前記した雲・降雨の超水滴モデルでは、実水滴が実粒子に当たり、湿潤大気が流体に当たる。そして、実水滴の属性は半径と凝結核の量である。速度は終端速度を与える事にしたので、これは独立な属性でない。また、半径の時間発展は数式（２）で与えられる。凝結核の量は衝突併合の際を除いて時間変化しないので、明示的に書けばｄＭ／ｄｔ＝０である。また、衝突効率はＥ_ｊｋ（Ｒ_ｊ，Ｒ_ｋ）として与えた。

また、超粒子の属性として、温度Ｔや帯電荷Ｃを加えることができる。実粒子の温度の時間発展方程式ｄＴ／ｄｔ＝・・・を取り扱うことや帯電荷の時間発展方程式ｄＣ／ｄｔ＝・・・を取り扱うことが可能である。なお、これらの属性により、衝突併合確率や衝突分裂確率に対する依存性（帯電すると衝突しやすくなったり衝突しにくくなったりすること）、つまり、衝突効率Ｅ_ｊｋ（Ｒ_ｊ，Ｍ_ｊ，Ｔ_ｊ，Ｃ_ｊ，Ｒ_ｋ，Ｍ_ｋ，Ｔ_ｋ，Ｃ_ｋ，流体変数）を明確にしておく必要がある。

また、雲・降雨の場合の実水滴のみならず、取り扱う対象を、雪、ひょう、あられに拡張し、これらをまとめて同種の実粒子と見なして、超粒子モデルを用いる。そのために、こられの状態を区別する属性Ｇを導入する。例えば、Ｇ＝１は水、Ｇ＝２は樹状結晶の雪、Ｇ＝３は粉雪等に属性Ｇによって区別することとする。この場合、属性Ｇの時間発展（つまり、どのように水が樹状結晶になったり、粉雪になったりするか等）を含めた当該属性の時間発展方程式群ｄ「Ａ」／ｄｔ＝「Ｆ」（「Ａ」）、衝突効率Ｅ_ｊｋ（「Ａ」_ｊ，「Ａ」_ｋ，流体変数）、さらに、衝突併合前後または衝突分裂前後での属性Ｇの変化を明確にしておく必要がある。

また、超粒子としてマイクロバブルを取り扱って、当該マイクロバブルの振る舞いを演算する場合について説明する。マイクロバブルとは微小な泡であり、この泡を実粒子として見なして、水を流体として想定し、超粒子モデルを用いる。泡の半径Ｒを属性の１つとし、当該泡の速度は浮力と空気抵抗とが釣り合っているとし終端速度で与えられる半径Ｒの関数とする。この場合、属性Ｒを時間発展方程式（泡が大きくなったり、小さくなったり、破裂したりすることについて記述した方程式）とし、衝突効率Ｅ_ｊｋ（Ｒ_ｊ，Ｒ_ｋ）とを明確にしておく必要がある。

また、超粒子として粉塵を取り扱って、当該粉塵の振る舞いを演算する場合について説明する。粉塵とは空気中を飛散している微小な粉末や塵（固体）であり、この粉塵を実粒子と見なして、空気を流体として想定し、超粒子モデルを用いる。小さな粉塵同士が空気中で衝突することで大きな粉塵になると取り扱う。この場合、粉塵の運動方程式（例えば、粉塵の大きさがある程度以上になると、自重により浮力と空気抵抗との釣り合いが崩れることについて記述した方程式）、衝突併合確率等を明確にしておく必要がある。

また、超粒子として液滴分散系(エマルジョン系)を取り扱って、当該液滴分散系の振る舞いを演算する場合について説明する。液滴分散系とは、親和結合しない液敵同士が分散していることを指し、例えば、水の中に油滴が浮遊している状況、簡便にいえば、良く撹拌されたドレッシング等や、重工業に用いられるエマルジョン燃料に該当する。油滴を実粒子と見なして、水を流体として想定し、超粒子モデルを用いる。この場合、油滴の運動方程式、油滴同士の衝突併合確率等を明確にしておく必要がある。

さらにまた、超粒子として燃料液滴を取り扱って、当該燃料液滴の振る舞い、いわゆる噴霧燃焼を演算する場合について説明する。この燃料液滴を実粒子と見なして、空気と燃料の混合液体を流体として想定し、超粒子モデルを用いる。例えば、属性群「Ａ」＝（「ｖ」，ｒ，Ｔ，ｑ，ｑ´）とし、液滴の速度「ｖ」、実効半径（換算半径）ｒ、温度Ｔ、変形率（球形からのズレ）ｑ、変形率の時間微分ｑ´を属性と取り扱う。この場合、これら各属性の時間発展方程式ｄ「Ａ」／ｄｔ＝「Ｆ」（「Ａ」）と、液滴同士の衝突併合確率等を明確にしておく必要がある。

出力手段１１ａは、入力手段３ａに入力された計算時間に関するデータ、時間ステップ（刻み時刻）と、総計算時間とに従って、演算手段１３で演算された結果を出力するものである。

このシミュレーション装置１Ａによれば、演算手段１３において、分割空間内で同じ属性を有している複数の実粒子の超粒子として取り扱って、超粒子モデルを、マイクロバブルの振る舞い、粉塵の振る舞い、液滴分散系の振る舞い、燃料液滴の振る舞いの演算に用いることで、実粒子を分布関数として取り扱うことなく、当該実粒子の属性数が増加しても、時間発展（経過時間に伴う変化）の計算時間の増加を抑制することができ、高い精度で予測することができる。

（超粒子法による場合のシミュレーションの概略動作）
次に、図１７に示すフローチャートを参照して、シミュレーション装置１Ａの動作について説明する（適宜、図１６参照）。なお、この動作において、超粒子属性更新手段１３ｂおよび流体場更新手段１３ｄによって、超粒子の属性の更新および流体場の更新は適宜行われていくこととして、説明を省略している。

まず、シミュレーション装置１Ａは、入力手段３ａに、初期変数が入力され、これにより、当該装置１で演算に必要とされる種々の変数が設定される（ステップＳ２１）。続いて、シミュレーション装置１Ａは、入力手段３ａによって入力された初期変数を、データ記憶手段５ａに記憶する。

そうすると、シミュレーション装置１Ａは、演算手段１３の超粒子運動演算手段１３ａによって、超粒子の速度を演算する（ステップＳ２２）。そして、シミュレーション装置１Ａは、超粒子の衝突併合が生じている時間ステップであるか否かを判定し（ステップＳ２３）、衝突併合が生じている時間ステップであると判定した場合（ステップＳ２３、Ｙｅｓ）には、超粒子衝突併合演算手段１３ｃによって、超粒子の衝突併合の演算を実行する（ステップＳ２４）。衝突併合が生じている時間ステップであると判定しなかった場合（ステップＳ２３、Ｎｏ）には、ステップＳ２５に移行する。

続けて、シミュレーション装置１Ａは、超粒子の衝突分裂が生じている時間ステップであるか否かを判定し（ステップＳ２５）、衝突分裂が生じている時間ステップであると判定した場合（ステップＳ２５、Ｙｅｓ）には、超粒子衝突分裂演算手段１３ｅによって、超粒子の衝突分裂の演算を実行する（ステップＳ２６）。衝突分裂が生じている時間ステップであると判定しなかった場合（ステップＳ２５、Ｎｏ）には、ステップＳ２７に移行する。

さらに続けて、シミュレーション装置１Ａは、超粒子の単体分裂が生じている時間ステップであるか否かを判定し（ステップＳ２７）、単体分裂が生じている時間ステップであると判定した場合（ステップＳ２７、Ｙｅｓ）には、超粒子単体分裂演算手段１３ｆによって、超粒子の単体分裂の演算を実行する（ステップＳ２８）。単体分裂が生じている時間ステップであると判定しなかった場合（ステップＳ２７、Ｎｏ）には、ステップＳ２９に移行する。

そして、シミュレーション装置１Ａは、時間ステップ（Δｔ）を加算して（ステップＳ２９）、加算した時間の合計が、総計算時間に達したか否かを判定する（ステップＳ３０）。そして、シミュレーション装置１Ａは、総計算時間に達したと判定しなかった場合（ステップＳ３０、Ｎｏ）、ステップＳ２２に戻り、動作を継続し、総計算時間に達したと判定した場合（ステップＳ３０、Ｙｅｓ）、動作を終了する。

（従来技術との顕著な差について）
次に、直近の従来技術であるＫＩＶＡ＋拡張ＮＴＣ法（D.P.Schmidt and C.J.Rutland,”A New Droplet Collision Algorithm,”J.Comput.Phys.,164,62-80(2000)）と本願との顕著な差について説明する。
このＫＩＶＡ＋拡張ＮＴＣ法と超粒子法（超水滴法）との相違点は、衝突併合操作の際のペアの作り方である。超粒子法では、ある分割空間内のＮ個の超粒子のうち、各超粒子の相手が一つに定まった“[Ｎ／２個]”のペアをつくるが、ＫＩＶＡ＋拡張ＮＴＣ法はそういった操作を行わない。

この違いにより、超粒子法では各分割空間内での衝突併合操作を、数値計算の際にベクトル化できる事である。ベクトル化とは、数値計算の際に演算を並列化する方法の一種であり、計算効率を大幅に向上させる(使用するコンピュータの機種に依存するが、例えば、２５６倍)ことが可能である。超粒子法では、各超粒子の衝突相手が一つに定まった“[Ｎ／２個]”のペアとなり、これらペア群の中に同一の超粒子は、一度ずつしか出てこず重複しない。それゆえ、超粒子法では、ベクトル化が可能である。これはＫＩＶＡ＋拡張ＮＴＣ法ではできない事である。

以上、本発明の実施形態について説明したが、本発明は前記実施形態には限定されない。例えば、本実施形態では、シミュレーション装置１として説明したが、当該装置１の各構成の処理を、一般的または特殊なコンピュータ言語で記述したシミュレーションプログラムとして捉えることも可能であり、さらに、当該装置１で実行される処理を１つずつの過程としたシミュレーション方法として捉えることも可能である。そして、これらはシミュレーション装置１と同様の効果を奏する。

本発明の実施形態に係るシミュレーション装置のブロック図である。 超水滴と実水滴との関係を表した模式図である。 超水滴の衝突併合と実水滴の衝突併合との対応関係を示した模式図である。 実水滴の半径と終端速度との関係を示した図である。 ある時刻において、実水滴同士がスウィープする体積について説明した図である。 図１に示したシミュレーション装置の概略の動作を示したフローチャートである。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 シミュレーション装置による演算結果を示した図である。 本発明の実施形態に係るシミュレーション装置（超粒子法に適用した場合）のブロック図である。 図１６に示したシミュレーション装置の動作を示したフローチャートである。

符号の説明

１、１Ａ シミュレーション装置
３、３ａ 入力手段
５、５ａ データ記憶手段
７ 微物理モデル演算手段
７ａ 超水滴運動演算手段
７ｂ 超水滴凝結成長演算手段
７ｃ 超水滴衝突併合演算手段
９ 流体力学モデル演算手段
１１、１１ａ 出力手段
１３ 演算手段
１３ａ 超粒子運動演算手段
１３ｂ 超粒子属性更新手段
１３ｃ 超粒子衝突併合演算手段
１３ｄ 流体場更新手段
１３ｅ 超粒子衝突分裂演算手段
１３ｆ 超粒子単体分裂演算手段

Claims (10)

1. 観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、
   前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、
   前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、
   予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、
   前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するシミュレーション方法であって、
   前記初期時刻と前記超粒子の属性と前記超粒子の総個数と前記体積と前記超粒子の速度と前記超粒子の位置座標と前記流体場変数とを初期変数として入力する入力ステップと、
   この入力ステップにて入力された初期変数に基づいて、前記経過時間に伴った実粒子の運動が前記体積と前記速度と前記位置座標と前記流体場変数とに従って前記属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の前記超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する演算ステップと、
   この演算ステップにより前記超粒子の属性、速度、位置座標、個数および多重度の演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する出力ステップと、
   を含むことを特徴とするシミュレーション方法。
2. 観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、
   前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、
   前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、
   予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、
   前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するために、コンピュータを、
   前記初期時刻と前記超粒子の属性と前記超粒子の総個数と前記体積と前記超粒子の速度と前記超粒子の位置座標と前記流体場変数とを初期変数として入力する入力手段、
   この入力手段で入力された初期変数に基づいて、前記経過時間に伴った実粒子の運動が前記体積と前記速度と前記位置座標と前記流体場変数とに従って前記属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の前記超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する演算手段、
   この演算手段により前記超粒子の属性、速度、位置座標、個数および多重度の演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する出力手段、
   として機能させることを特徴とするシミュレーションプログラム。
3. 観測される空間内に存在する実粒子において当該実粒子同士が当該空間内における所定区画とする体積内で所定時間内に所定確率で衝突する際に、
   前記実粒子が、任意個数の属性と初期時刻における当該属性の１つである速度と初期時刻における前記空間上の位置座標とで表され、
   前記空間内を満たす流体が前記初期時刻からの経過時間と前記空間との関数で表される任意個数の流体場変数によって特徴付けられ、
   予め設定した所定の同属性を有する前記実粒子の任意数ごとの集合を超粒子とし、当該任意数を当該超粒子の多重度とし、当該超粒子が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突した際に前記多重度が変化する場合に、
   前記超粒子に関するデータを演算することで、任意時間後の前記実粒子に関するデータを出力するシミュレーション装置であって、
   前記初期時刻と前記超粒子の属性と前記超粒子の総個数と前記体積と前記超粒子の速度と前記超粒子の位置座標と前記流体場変数とを初期変数として入力する入力手段と、
   この入力手段で入力された初期変数に基づいて、前記経過時間に伴った実粒子の運動が前記体積と前記速度と前記位置座標と前記流体場変数とに従って前記属性ごとに規定される属性時間発展方程式、前記実粒子の速度と位置座標との関係を規定した位置座標時間発展方程式および所定時間内に前記超粒子同士が前記確率で衝突するとしたモンテカルロ演算により、衝突した後の前記超粒子の属性、速度、位置座標、個数および多重度を演算すると共に、前記流体の変化が前記流体場変数および前記超粒子の属性、速度、位置座標、個数および多重度に基づいて規定される流体場時間発展方程式により、前記流体場変数を演算する演算手段と、
   この演算手段により前記超粒子の属性、速度、位置座標、個数および多重度の演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記実粒子の任意時間後の属性、速度、位置座標および個数に関するデータとして出力すると共に、前記流体に関するデータを出力する出力手段と、
   を備えることを特徴とするシミュレーション装置。
4. 観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、
   前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、
   予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、
   前記超水滴に関するデータを演算することで、任意時間後の前記実水滴に関するデータを出力するシミュレーション方法であって、
   前記初期時刻と前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積と前記超水滴の位置座標と当該分割空間ごとの前記実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する変数入力ステップと、
   前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積とに基づいて、前記全体空間内における前記超水滴の運動による位置変化、前記超水滴の凝結成長による水量変化および前記超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る微物理モデル演算ステップと、
   この微物理モデル演算ステップにて得られた前記相互作用量および前記周辺環境データに基づいて、前記実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を前記微物理モデル演算ステップにフィードバックする流体力学モデル演算ステップと、
   この流体力学モデル演算ステップと前記微物理モデル演算ステップとによる演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記任意時間後の前記実水滴に関するデータとして出力すると共に、前記任意時間後の前記周辺環境データを出力する出力ステップと、
   を含むことを特徴とするシミュレーション方法。
5. 観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、
   前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、
   予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、
   前記超水滴に関するデータを演算することで、任意時間後の前記実水滴に関するデータを出力するために、コンピュータを、
   前記初期時刻と前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積と前記超水滴の位置座標と当該分割空間ごとの前記実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する変数入力手段、
   前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積とに基づいて、前記全体空間内における前記超水滴の運動による位置変化、前記超水滴の凝結成長による水量変化および前記超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る微物理モデル演算手段、
   この微物理モデル演算手段で得られた前記相互作用量および前記周辺環境データに基づいて、前記実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を前記微物理モデル演算手段にフィードバックする流体力学モデル演算手段、
   この流体力学モデル演算手段と前記微物理モデル演算手段とによる演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記任意時間後の前記実水滴に関するデータとして出力すると共に、前記任意時間後の前記周辺環境データを出力する出力手段、
   として機能させることを特徴とするシミュレーションプログラム。
6. 前記微物理モデル演算手段は、
   前記超水滴にかかる重力および空気抵抗が釣り合った状態で、前記超水滴の運動が風速によって変化し、当該風速に対して一定の相対速度である終端速度で運動するものとして、当該終端速度を演算する超水滴運動演算手段と、
   前記超水滴に含まれる水量が前記周辺環境データに含まれる湿度によって変化するものとして、当該水量を演算する超水滴凝結成長演算手段と、
   前記超水滴同士が衝突する任意のペアを設定し、当該ペアが衝突併合する確率を衝突併合確率とし、前記ペアの数を所定数に削減すると共に前記衝突併合確率を所定幅上昇させて、前記超水滴同士の衝突併合による超水滴の属性、多重度および個数を演算する超水滴衝突併合演算手段と、
   を有することを特徴とする請求項５に記載のシミュレーションプログラム。
7. 超水滴衝突併合演算手段は、前記衝突併合の過程をモンテカルロ法による数値シミュレーションに従って演算することを特徴とする請求項６に記載のシミュレーションプログラム。
8. 観測される全体空間内に存在する実水滴において当該実水滴同士が所定体積内で所定時間内に所定確率で衝突する際に、
   前記実水滴が、任意個数の属性と初期時刻における前記全体空間を分割した分割空間内の位置座標とで表され、
   予め設定した所定の同属性を有する前記実水滴の任意数ごとの集合を超水滴とし、当該任意数を当該超水滴の多重度とし、当該超水滴が前記所定確率を基準に前記多重度に応じた確率で衝突し、衝突して併合した際に前記多重度が変化する場合に、
   前記超水滴に関するデータを演算することで、任意時間後の前記実水滴に関するデータを出力するシミュレーション装置であって、
   前記初期時刻と前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積と前記超水滴の位置座標と当該分割空間ごとの前記実水滴の周辺環境に関するデータである周辺環境データとを、初期変数として入力する変数入力手段と、
   前記超水滴の属性と前記超水滴の総個数と前記全体空間の体積と前記分割空間の体積とに基づいて、前記全体空間内における前記超水滴の運動による位置変化、前記超水滴の凝結成長による水量変化および前記超水滴同士の衝突併合による超水滴の属性、多重度および個数変化を演算し、当該超水滴に基づいて規定される実水滴の質量を求め、当該実水滴の質量に基づき前記超水滴から大気への相互作用量を得る微物理モデル演算手段と、
   この微物理モデル演算手段で得られた前記相互作用量および前記周辺環境データに基づいて、前記実水滴の存在する大気の流体力学過程の演算を行うと共に、この演算結果を前記微物理モデル演算手段にフィードバックする流体力学モデル演算手段と、
   この流体力学モデル演算手段と前記微物理モデル演算手段とによる演算を、前記任意時間に達するまで繰り返した後に、繰り返した後の結果を前記任意時間後の前記実水滴に関するデータとして出力すると共に、前記任意時間後の前記周辺環境データを出力する出力手段と、
   を備えることを特徴とするシミュレーション装置。
9. 前記微物理モデル演算手段は、
   前記超水滴にかかる重力および空気抵抗が釣り合った状態で、前記超水滴の運動が風速によって変化し、当該風速に対して一定の相対速度である終端速度で運動するものとして、当該終端速度を演算する超水滴運動演算手段と、
   前記超水滴に含まれる水量が前記周辺環境データに含まれる湿度によって変化するものとして、当該水量を演算する超水滴凝結成長演算手段と、
   前記超水滴同士が衝突する任意のペアを設定し、当該ペアが衝突併合する確率を衝突併合確率とし、前記ペアの数を所定数に削減すると共に前記衝突併合確率を所定幅上昇させて、前記超水滴同士の衝突併合よる超水滴の属性、多重度および個数を演算する超水滴衝突併合演算手段と、
   を有することを特徴とする請求項８に記載のシミュレーション装置。
10. 前記超水滴衝突併合演算手段は、前記衝突併合の過程をモンテカルロ法による数値シミュレーションに従って演算することを特徴とする請求項９に記載のシミュレーション装置。

Images