(57)【要約】【課題】  界面形状を表現する線分列の各線分の頂点の順番が入れ替わるためにループを作ってしまう現象が現れるという不都合を生じない界面軌跡シミュレーション方法を提供する。【解決手段】  界面追跡シミュレーション方法にオイラー描像的設計思想を取り入れ、界面に沿った距離を、２次元ユークリッド空間から自然に誘導されるものとして定義し、界面曲線の微小線素の長さを界面変形に対して不変とし、界面変形を表現する速度場の算出にあたって等長的な界面上での関数であるという性質を利用した。これにより、時間発展中の任意の時刻の界面形状において、左右非対称性を生ずることなく界面形状を追跡するためのシミュレーションが可能となる。【選択図】  図６
2005183411.tif 000002
【特許請求の範囲】
【請求項１】
２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーション方法であって、
界面を表現する曲線を、等長の線分に分割してなる線分列により近似し、
前記線分の長さを維持しながら、次の時間ステップの界面形状を決定する、
ことを特徴とするコンピュータ・シミュレーション方法。
【請求項２】
前記線分列の各頂点を、整数ｎ＝１，…，  Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数１】
2005183411.tif 000066
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数２】
2005183411.tif 000067
により定義し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応した前記線分列のｎ＝０及びｎ＝Ｎの周りでの振る舞いが、解析的に知られた漸近的振る舞いを近似できるように、Ｖ０の値を決定し、該法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応して、
【数３】
2005183411.tif 000068
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【数４】
2005183411.tif 000069
により算出すること、
を特徴とする請求項１に記載のコンピュータ・シミュレーション方法。
【請求項３】
前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数５】
2005183411.tif 000070
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数６】
2005183411.tif 000071
により定義し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【数７】
2005183411.tif 000072
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【数８】
2005183411.tif 000073
により算出し、
前記中間的な接角を用いて界面の更新を行い、
前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、
前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【数８】
2005183411.tif 000074
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【数９】
2005183411.tif 000075
として求め、
接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新する、
ことを特徴とする請求項１に記載のコンピュータ・シミュレーション方法。
【請求項４】
２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーション方法であって、
界面を表現する曲線を向きづけ、
前記界面曲線を等長の線分に分割近似し、
前記線分の長さを維持しながら、次の時間ステップの界面形状を決定し、
各時間ステップ毎に前記曲線の向きの入れ替えを行う、
ことを特徴とするコンピュータ・シミュレーション方法。
【請求項５】
前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数１０】
2005183411.tif 000076
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数１１】
2005183411.tif 000077
により定義し、
Ｖ０の値をゼロとして、
【数１２】
2005183411.tif 000078
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【数１３】
2005183411.tif 000079
により算出し、
各時間ステップ毎に、前記接角に割り振られた番号の向きが逆になるように前記線分列の各頂点の番号を入れ替えて界面形状を求める、
ことを特徴とする請求項４に記載のコンピュータ・シミュレーション方法。
【請求項６】
前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数１４】
2005183411.tif 000080
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数１５】
2005183411.tif 000081
により定義し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の小さい順に
【数１６】
2005183411.tif 000082
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【数１７】
2005183411.tif 000083
により算出し、
前記中間的な接角を用いて界面の更新を行い、
前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、
前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、
ｎ＝Ｎに対応する線分の接線方向の速度ＶＮを適当に決定し、
接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の大きな順に
【数１８】
2005183411.tif 000084
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【数１９】
2005183411.tif 000085
として求め、
接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新する、
ことを特徴とする請求項４に記載のコンピュータ・シミュレーション方法。
【請求項７】
コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、
界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数２０】
2005183411.tif 000086
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数２１】
2005183411.tif 000087
により定義し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応した前記線分列のｎ＝０及びｎ＝Ｎの周りでの振る舞いが、解析的に知られた漸近的振る舞いを近似できるように、Ｖ０の値を決定し、該法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応して、
、
【数２２】
2005183411.tif 000088
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【数２３】
2005183411.tif 000089
とする、
ことを特徴とするプログラム。
【請求項８】
コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、
界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数２４】
2005183411.tif 000090
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数２５】
2005183411.tif 000091
により定義し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【数２６】
2005183411.tif 000092
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【数２７】
2005183411.tif 000093
により算出し、
前記中間的な接角を用いて界面の更新を行い、
前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、
前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【数２８】
2005183411.tif 000094
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【数２９】
2005183411.tif 000095
として求め、
接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新する、
ことを特徴とするプログラム。
【請求項９】
コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、
界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数３０】
2005183411.tif 000096
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数３１】
2005183411.tif 000097
により定義し、
Ｖ０の値をゼロとして、
【数３２】
2005183411.tif 000098
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【数３３】
2005183411.tif 000099
により算出し、
各時間ステップ毎に、前記接角に割り振られた番号の向きが逆になるように前記線分列の各頂点の番号を入れ替えて界面形状を求める、
ことを特徴とするプログラム。
【請求項１０】
コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、
界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、
前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、
前記接角を用いて
【数３４】
2005183411.tif 000100
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【数３５】
2005183411.tif 000101
により定義し、
ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、
各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の小さい順に
【数３６】
2005183411.tif 000102
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【数３７】
2005183411.tif 000103
により算出し、
前記中間的な接角を用いて界面の更新を行い、
前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、
前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、
ｎ＝Ｎに対応する線分の接線方向の速度ＶＮを適当に決定し、
接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の大きな順に
【数３８】
2005183411.tif 000104
を満足するように決定し、
前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【数３９】
2005183411.tif 000105
として求め、
接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新する、
ことを特徴とするプログラム。
【請求項１１】
請求項７〜１０いずれかに記載のプログラムを記憶させ、コンピュータ読取可能な記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、界面追跡シミュレーション方法、および、そのシミュレーション実行用プログラムに関し、より詳細には、半導体プロセス、結晶成長、２流体間界面、酸化還元反応を伴う界面等において、これらの界面形状の時間的発展を、コンピュータを用いてシミュレーションする方法、および、その実行用プログラムに関する。
【背景技術】
【０００２】
物理現象を理解するうえで、特定の２次元平面内に形成された界面が時間的に発展する様子を、コンピュータ・シミュレーションにより解析することが必要とされる場合がある。その典型例が、半導体プロセスシミュレーション（例えば、特許文献１、および、特許文献２参照）であり、産業的観点からも極めて重要な意義をもつ技術である。
【０００３】
例えば、特許文献１には、半導体プロセスの１過程である酸化膜成長過程において、酸化膜の成長に伴って微妙に高さの異なる細かな凹凸のある酸化膜同士に重なりが生じた場合でも、酸化膜成長の時間刻みが不必要に細分化されることがなく、このため時間的進行に停滞が生じないこととした酸化シミュレーション方法が記載されている。しかし、このような従来のシミュレーション方法には、以下のような問題がある。
【０００４】
図１１（ａ）、（ｂ）は、長さの異なる複数の線分からなる線分列によって界面形状を表現したものである。各図中の黒丸は、界面の形状を複数の線分に分割するためのポイント（頂点）であり、相互に隣接するポイント（頂点）同士を直線で結んで得られる線分を線分列として形成し、界面の形状を近似している。
【０００５】
界面の形状は時々刻々変化するが、この時間変化は、各頂点に与えられる移動ベクトルによって特徴づけられる。すなわち、各頂点において２つの線分がなす角度の２等分角方向に、単位時間刻み当たりの移動距離に対応するノルムをもったベクトルを定義し、これをこの頂点の移動ベクトルとするのである。その結果、単位時間刻み経過後の界面は、図中に細線で示した形状をとることとなる。
【０００６】
図１２（ａ）〜（ｆ）は、界面形状の経時変化の様子を説明するための図であり、各時刻ｔ１、ｔ２´、およびｔ２における界面形状を示している。時刻ｔ１において図１２（ａ）に示した形状を有する界面は、時刻ｔ２´において図１２（ｂ）、時刻ｔ２において図１２（ｃ）の形状へと変化する。同様に、時刻ｔ１において図１２（ｄ）に示した形状を有する界面は、時刻ｔ２´において図１２（ｅ）、時刻において図１２（ｆ）の形状へと変化する。
【０００７】
この場合、図１２（ｃ）に示した界面形状では、近傍の線分列の頂点の順番が入れ替わるためにループを作ってしまう現象（現象１）が起こっており、一方、図１２（ｆ）に示した界面形状では、異なる２つの線分列を構成する線分が互いに交差するという現象（現象２）が起こっている。
【０００８】
また、上述した特許文献２に記載のシミュレーション方法によれば、線分の交差現象が時刻ｔ２で生じた場合には、先ず、シミュレーションパラメータである時間刻みを小さく設定し直し、図１２（ａ）又は（ｄ）に示した時刻ｔ１での界面形状を基に、時刻ｔ２に先行する時刻であるｔ２´までシミュレーションを実行し、線分の交差が生じていない図１２（ｂ）又は（ｅ）に示した界面形状を得る。そして、図１２（ｂ）又は（ｅ）に示した界面形状を表現する各頂点の再配置等による初期化処理を行った後に、再度シミュレーションを実行し直すという取り扱いがなされる。
【０００９】
【特許文献１】特開平５−８５８９１号公報
【特許文献２】特開平９−１９９４９２号公報
【非特許文献１】Ａｒｎｏｌｄ  Ｖ．  Ｉ．，“Ｍａｔｈｅｍａｔｉｃａｌ  Ｍｅｔｈｏｄ  ｏｆ  Ｃｌａｓｓｉｃａｌ  Ｍｅｃｈａｎｉｃｓ”，  Ｓｅｃｏｎｄ  Ｅｄｉｔｉｏｎ，  Ｓｐｒｉｎｇｅｒ  （１９８９）.
【非特許文献２】Ｒ．Ｃ．Ｂｒｏｗｅ  ｅｔ  ａｌ．，  “Ｇｅｏｍｅｔｒｉｃａｌ  ｍｏｄｅｌｓ  ｏｆ  ｉｎｔｅｒｆａｃｅ  ｅｖｏｌｕｔｉｏｎ”  Ｐｈｙｓｉｃａｌ  Ｒｅｖｉｅｗ  Ａ，  ｖｏｌ．２９,  ｐ．１３３５−１３４２  （１９８４）．
【非特許文献３】Ｓ．Ｓａｉｔｏ  ｅｔ  ａｌ．，Ｊｏｕｒｎａｌ  ｏｆ  Ｐｈｙｓｉｃｓ  Ａ：  Ｍａｔｈｅｍａｔｉｃａｌ  ａｎｄ  Ｇｅｎｅｒａｌ，  ｖｏｌ．３０，  ｐｐ．６９５１  （１９９７）．
【発明の開示】
【発明が解決しようとする課題】
【００１０】
しかしながら、上述したようなシミュレーション方法によれば、界面形状を表現する線分列を構成する線分の交差が生じるたびに、シミュレーションパラメータである時間刻みを設定し直して、再度シミュレーションをやり直す必要があり、シミュレーションに長時間を要してしまうという問題がある。
【００１１】
線分の交差が生じてしまう現象の本質は、従来のシミュレーション方法では、流体の流れ解析で広く採用されるラグランジュの流れ関数を基礎とする計算体系（ラグランジュ描像）を採用しているという点にある。
【００１２】
流体の流れを解析するためのアプローチには、ラグランジュ描像とオイラー描像という２つの見方がある。具体的には、ラグランジュ描像では流体中を運動する粒子の軌跡を記述する運動方程式から出発するのに対し、オイラー描像は流体中に形成される「速度場」を仮定し、この速度場の時間的変化を記述する運動方程式から出発するのである。
【００１３】
これを数学的に解釈すると、ラグランジュ描像では、位相空間内において定義された、ある群の作用を受ける被作用子の軌跡を時間的に追跡するのに対して、オイラー描像では、群の作用そのものの時間的な変化を追跡するということになる。
【００１４】
例えば、２次元完全流体に対するオイラー方程式は、速度場をｖ、圧力をｐ、時間をｔとして、
【００１５】
【数１】
2005183411.tif 000003
【００１６】
で与えられる。
【００１７】
このオイラー方程式を、位相空間内での群の作用という観点から解釈しなおすと、２次元空間内に形成された、ある領域の全体にわたって定義される滑らかな関数の集合において、その集合を構成するある要素から他の要素へ作用する、微分同相群の作用について表現した方程式であると見なすことが可能である（例えば、非特許文献１）。
【００１８】
このような数学的解釈の基に、基底空間として、変形のパラメータである時間ｔの軸を仮定し、ファイバーとして、２次元ユークリッド空間内の所望の領域でのファイバー束を仮定して、これらの基底空間とファイバーから構成される多様体を考える。
【００１９】
各ファイバーへの作用が微分同相であれば、ファイバー束を微分同相群のフレーム束として取り扱うことが可能であり、微分同相群の主ファイバー束として取り扱うことも可能である。これをフレーム束として取り扱うこととするのがラグランジュ描像に対応し、主ファイバー束として取り扱うこととするのがオイラー描像に対応する。
【００２０】
このような取り扱いにおいては、上述したオイラー方程式は、２次元ユークリッド空間上の微分同相群を、調和写像理論（あるいは非線形シグマ模型）におけるエネルギー汎関数と仮定し、このエネルギー汎関数に最小原理を適用して導かれた運動方程式として理解されることになる。また、上述した２次元ユークリッド空間内に形成される速度場は、無限次元リー群である微分同相群の、リー環の生成元として理解される。更に、ラグランジュ微分は、２次元ユークリッド空間での幾何学的接続として理解される。
【００２１】
ラグランジュ描像を採用すると、生起する事象を直感的に理解し易くなるという利点がある一方、位相を定義することが困難となるという問題がある。例えば、速度場の配位の状態によっては、着目した流体粒子の近傍点が、ある時間の経過後に、その近傍系に移動しないという現象が起こり得る。
【００２２】
図１３は、上述した問題点を説明するための図である。理解を容易にするために、流体粒子の運動場である２次元平面を等間隔の格子で分割している。
【００２３】
図中に示した太い矢印は、各格子点での速度場をベクトル表現したものである。また、細線の曲線は、各流体粒子の軌跡を意味している。すなわち、時間の経過に伴って、ａの流体粒子はａ´へ、ｂの流体粒子はｂ´へ、ｃの流体粒子はｃ´へと移動している。各流体粒子間の相対的な距離関係の時間変化に着目すると、初期状態では、ａはｂの近傍にあり、ｃはｂから遠くに存在するが、ある時間の経過後には、ａ´はｂ´の近傍になく、ｃ´はｂ´の近傍に存在しており、各流体粒子の相対的距離関係が、時間的に変化している。
【００２４】
しかし、オイラー描像を採用すれば、かかる不自然さは解消される。すなわち、ファイバーの位相は、２次元ユークリッド空間から誘導される距離空間となり、また、速度場は、２次元ユークリッド空間上に形成された場として取り扱うことが可能である（ここで、速度場の配位空間に如何なる位相を与えるかは別問題である）。
【００２５】
従来のシミュレーション方法で発生する上述した２つの現象（現象１及び現象２）は、ラグランジュ描像が内包する問題に起因するものである。
【００２６】
しかし、これらの現象のうち、現象２は、従来の界面追跡シミュレーションの設計思想が、広い意味でのラグランジュ的描像に基礎を置いていることに起因して生じる現象である。即ち、現象２は、従来の界面追跡シミュレーションの設計思想自身がもつ特徴であると考えるべきである。また、ある種の界面現象では、現象２のような物理現象が起こりにくい場合があることが知られている。
【００２７】
他方、現象１は、シミュレーション方法の未熟さに起因するものと考えることができる。すなわち、現象１は、従来のシミュレーション方法が、時間的に変化する界面の形状を数学的に表現できていないことを意味しており、致命的な欠陥である。
【００２８】
本発明は、このような問題に鑑みてなされたもので、その目的とするところは、シミュレーション方法に、オイラー描像的設計思想を取り入れることにより、上述した現象１を克服する界面追跡シミュレーション方法を提供することにある。従って、オイラー描像に従って、界面に沿った距離を、２次元ユークリッド空間から自然に誘導されるものとして定義する。この場合、界面曲線の微小線素の長さは、界面変形に対して不変とし、界面変形を表現する速度場は、等長的な界面上での関数であることを利用する。
【課題を解決するための手段】
【００２９】
本発明は、このような目的を達成するために、請求項１に記載の発明は、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーション方法であって、界面を表現する曲線を、等長の線分に分割してなる線分列により近似し、前記線分の長さを維持しながら、次の時間ステップの界面形状を決定することを特徴とするコンピュータ・シミュレーション方法である。
【００３０】
また、請求項２に記載の発明は、請求項１に記載のコンピュータ・シミュレーション方法において、前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００３１】
【数２】
2005183411.tif 000004
【００３２】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００３３】
【数３】
2005183411.tif 000005
【００３４】
により定義し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応した前記線分列のｎ＝０及びｎ＝Ｎの周りでの振る舞いが、解析的に知られた漸近的振る舞いに近似できるように、Ｖ０の値を決定し、該法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応して、
【００３５】
【数４】
2005183411.tif 000006
【００３６】
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【００３７】
【数５】
2005183411.tif 000007
【００３８】
とすることを特徴とするコンピュータ・シミュレーション方法である。
【００３９】
また、請求項３に記載の発明は、請求項１に記載のコンピュータ・シミュレーション方法において、前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００４０】
【数６】
2005183411.tif 000008
【００４１】
を求め、
前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００４２】
【数７】
2005183411.tif 000009
【００４３】
により定義し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【００４４】
【数８】
2005183411.tif 000010
【００４５】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【００４６】
【数９】
2005183411.tif 000011
【００４７】
により算出し、前記中間的な接角を用いて界面の更新を行い、前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【００４８】
【数１０】
2005183411.tif 000012
【００４９】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【００５０】
【数１１】
2005183411.tif 000013
【００５１】
として求め、接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新することを特徴とするコンピュータ・シミュレーション方法である。
【００５２】
また、請求項４に記載の発明は、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーション方法であって、界面を表現する曲線を向きづけ、前記界面曲線を等長の線分に分割近似し、前記線分の長さを維持しながら、次の時間ステップの界面形状を決定し、各時間ステップ毎に前記曲線の向きの入れ替えを行うことを特徴とする。
【００５３】
また、請求項５に記載の発明は、請求項４に記載のコンピュータ・シミュレーション方法において、前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００５４】
【数１２】
2005183411.tif 000014
【００５５】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００５６】
【数１３】
2005183411.tif 000015
【００５７】
により定義し、Ｖ０の値をゼロとして、
【００５８】
【数１４】
2005183411.tif 000016
【００５９】
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【００６０】
【数１５】
2005183411.tif 000017
【００６１】
により算出し、各時間ステップ毎に、前記接角に割り振られた番号の向きが逆になるように前記線分列の各頂点の番号を入れ替えて界面形状を求めることを特徴とする。
【００６２】
また、請求項６に記載の発明は、請求項４に記載のコンピュータ・シミュレーション方法において、前記線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００６３】
【数１６】
2005183411.tif 000018
【００６４】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００６５】
【数１７】
2005183411.tif 000019
【００６６】
により定義し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の小さい順に
【００６７】
【数１８】
2005183411.tif 000020
【００６８】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【００６９】
【数１９】
2005183411.tif 000021
【００７０】
により算出し、前記中間的な接角を用いて界面の更新を行い、前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、ｎ＝Ｎに対応する線分の接線方向の速度ＶＮを適当に決定し、接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の大きな順に
【００７１】
【数２０】
2005183411.tif 000022
【００７２】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【００７３】
【数２１】
2005183411.tif 000023
【００７４】
として求め、接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新することを特徴とする。
【００７５】
また、請求項７に記載の発明は、コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００７６】
【数２２】
2005183411.tif 000024
【００７７】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００７８】
【数２３】
2005183411.tif 000025
【００７９】
により定義し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応した前記線分列の漸近的な振る舞いが、ｎ−１とｎ＋１の入替えに対して対称となるようにＶ０の値を決定し、
【００８０】
【数２４】
2005183411.tif 000026
【００８１】
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【００８２】
【数２５】
2005183411.tif 000027
【００８３】
とすることを特徴とするプログラムである。
【００８４】
また、請求項８に記載の発明は、コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００８５】
【数２６】
2005183411.tif 000028
【００８６】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【００８７】
【数２７】
2005183411.tif 000029
【００８８】
により定義し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【００８９】
【数２８】
2005183411.tif 000030
【００９０】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【００９１】
【数２９】
2005183411.tif 000031
【００９２】
により算出し、前記中間的な接角を用いて界面の更新を行い、前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、
【００９３】
【数３０】
2005183411.tif 000032
【００９４】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【００９５】
【数３１】
2005183411.tif 000033
【００９６】
として求め、接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新することを特徴とする。
【００９７】
また、請求項９に記載の発明は、コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【００９８】
【数３２】
2005183411.tif 000034
【００９９】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【０１００】
【数３３】
2005183411.tif 000035
【０１０１】
により定義し、Ｖ０の値をゼロとして、
【０１０２】
【数３４】
2005183411.tif 000036
【０１０３】
を満足するように接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を
【０１０４】
【数３５】
2005183411.tif 000037
【０１０５】
により算出し、各時間ステップ毎に、前記接角に割り振られた番号の向きが逆になるように前記線分列の各頂点の番号を入れ替えて界面形状を求めることを特徴とする。
【０１０６】
また、請求項１０に記載の発明は、コンピュータを用いて、２次元平面内に形成された界面の時間的発展を追跡するためのコンピュータ・シミュレーションを行なうためのプログラムであって、界面を近似する線分列の各頂点を、整数ｎ＝１，…，Ｎによって指定し、前記各線分の基準角からの接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝を算出し、前記接角を用いて
【０１０７】
【数３６】
2005183411.tif 000038
【０１０８】
を求め、前記各線分の接線ベクトルｔｎおよび法線ベクトルｎｎを、
【０１０９】
【数３７】
2005183411.tif 000039
【０１１０】
により定義し、ｎ＝０に対応する線分の接線方向の速度Ｖ０を適当に決定し、各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝に対応する接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の小さい順に
【０１１１】
【数３８】
2005183411.tif 000040
【０１１２】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の中間的な接角を、
【０１１３】
【数３９】
2005183411.tif 000041
【０１１４】
により算出し、前記中間的な接角を用いて界面の更新を行い、前記中間的な接角を、新たな接角｛φｎ｜ｎ＝０，…，Ｎ，φ０＝φ１，φＮ＝φＮ−１｝として定義し直し、前記更新後の界面を近似する各線分の法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝を算出し、ｎ＝Ｎに対応する線分の接線方向の速度ＶＮを適当に決定し、接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を、ｎ値の大きな順に
【０１１５】
【数４０】
2005183411.tif 000042
【０１１６】
を満足するように決定し、前記法線方向の速度｛Ｕｎ｜ｎ＝０，…，Ｎ｝と、前記接線方向の速度｛Ｖｎ｜ｎ＝０，…，Ｎ｝を基に、微少時間間隔εに対して、ｎ番目の新たな接角を、
【０１１７】
【数４１】
2005183411.tif 000043
【０１１８】
として求め、接線ベクトルｔｎ、法線ベクトルｎｎ、および接角を基に界面形状を更新することを特徴とする。
【０１１９】
更に、請求項１１に記載の発明は、請求項７〜１０に記載のプログラムを記憶し、コンピュータ読取可能な記憶媒体である。
【発明の効果】
【０１２０】
本発明の界面追跡シミュレーション方法では、オイラー描像的設計思想を取り入れることにより、界面に沿った距離を、２次元ユークリッド空間から自然に誘導されるものとして定義し、界面曲線の微小線素の長さを界面変形に対して不変とし、界面変形を表現する速度場の算出にあたって等長的な界面上での関数であるという性質を利用したので、界面形状を表現する線分列の各線分の頂点の順番が入れ替わるためにループを作ってしまう現象が現れるという不都合を生じないシミュレーション方法の提供が可能となる。
【発明を実施するための最良の形態】
【０１２１】
まず、簡単のために、取り扱うべき界面の形状を連続する１本の曲線として取り扱う場合について説明する。この場合は、界面の解析的取り扱いの観点からの研究がなされており、既に非特許文献２において報告されている。
【０１２２】
図１は、滑らかな１本の曲線で界面形状を表現した図である。この図において、曲線の各点の位置ベクトルをＸとすると、この曲線は、各点毎に定められる、接ベクトルｔと、法線ベクトルｎとによって特徴づけられる。ここで、ｓを弧長とすると、法線ベクトルは、ｎ＝ｄＸ／ｄｓで与えられる。
【０１２３】
これらのベクトルは、いわゆる、フレネ・セレの関係式、
【０１２４】
【数４２】
2005183411.tif 000044
【０１２５】
を通して関係づけられる。ここでκは曲線の曲率であり、曲率は曲率半径の逆数として定義されたものである。従って、ある基準角からの接角の角度をφとすると、接線ベクトルｔ、および法線ベクトルｎは、それぞれ、ｔ＝（ｃｏｓφ、ｓｉｎφ）、及び、ｎ＝（−ｓｉｎφ、ｃｏｓφ）、で与えられ、曲率は、κ＝ｄφ／ｄｓ、で与えられることとなる。
【０１２６】
次に、この曲線の等長変形を考える。変形のパラメータをｔとして、
【０１２７】
【数４３】
2005183411.tif 000045
【０１２８】
で与えられるものとする。
このとき、
【０１２９】
【数４４】
2005183411.tif 000046
【０１３０】
という条件を課すと、曲線の局所的な長さが不変な状態での変形を与えることとなる。ここで、
【０１３１】
【数４５】
2005183411.tif 000047
【０１３２】
という略記法を使用した。
式（３）に対して、微分処理を施すと
【０１３３】
【数４６】
2005183411.tif 000048
【０１３４】
となる。
従って、法線方向の「速度」Ｕを与え、
【０１３５】
【数４７】
2005183411.tif 000049
【０１３６】
で定義される「速度」Ｖを用いて曲線の変形を考えることとすれば、この変形は等長変形を意味する。
【０１３７】
次に、取り扱うべき対象の離散化（あるいは差分化とも呼ぶ）を施す。つまり、取り扱うべき対象を有限（あるいは加算）個の集合で近似表現する。この離散化には様々な方法が知られており、連続極限をとることにより微分方程式となる差分方程式は、連続濃度の数で存在することが知られている。
【０１３８】
例えば、カオスを形成することで有名なロジスティック・マップという差分方程式（漸化式）は、連続極限をとるとカオスを形成しない微分方程式であるロジスティック方程式となることが知られており、また、連続極限の後にこの微分方程式となる差分方程式であって、カオス信号を発生しない差分方程式があることも知られている（例えば、非特許文献３を参照）。すなわち、離散化という過程は、非常に豊富な情報を生成する過程であり、振る舞いの「悪い」離散化を行うことによって系が不安定化することがあり得ることを意味している。（もちろん、離散化によって連続の際、不安定であった系が安定化する場合もある。）
【０１３９】
従って、本発明の目的を達するためには、式（５）を離散化する必要があるが、単に、中心差分あるいは前進差分等によって離散化すればよい訳ではない。以下で述べるように、連続的な描像で得られた個々の幾何学量は、それぞれが幾何学的背景に基づいて定義可能な量であり、その背景を無視して離散化あるいは離散化を行うと、ある有限の時間発展後にシミュレーション結果が不安定化することが予想される。
【０１４０】
そこで、本発明のシミュレーション方法では、上述した連続幾何学的概念を、シミュレーション結果に矛盾を生ぜしめることがないように、離散幾何学的な立場に立脚して離散化し、曲線の変形を記述する差分方程式のスキームを与える。
【０１４１】
図２は、１本の曲線に離散化を施して得られた線分列を示す。この図で示すように、各線分は長さＬを持っており、並びの番号ｎ（ｎ＝０，１，２，…Ｎ）で区別される。また、各線分は、各頂点において定まる、位置ベクトルｒｎ、接線ベクトルｔｎ、法線ベクトルｎｎ、および接角φｎの各パラメータによって特徴づけられている。
【０１４２】
これらのパラメータは、以下の離散フレネ・セレの関係式を満足する。
【０１４３】
【数４８】
2005183411.tif 000050
【０１４４】
この線分列において、仮想的に以下の微小変形を考える。（Ｕｎに関しては図１８に記してある。）
【０１４５】
【数４９】
2005183411.tif 000051
【０１４６】
また、
【０１４７】
【数５０】
2005183411.tif 000052
【０１４８】
を定義すると、
【０１４９】
【数５１】
2005183411.tif 000053
【０１５０】
は、個々の線分の長さを維持するための条件となる。
式（７）と式（１０）の関係式から、
【０１５１】
【数５２】
2005183411.tif 000054
【０１５２】
が得られる。
更に、式（１１）の左辺の最後の項の（ｔｎ−１，ｎｎ−１）を、（ｔｎ，ｎｎ）に分解し、
【０１５３】
【数５３】
2005183411.tif 000055
【０１５４】
として、φｎεを導入すると、線分の長さを変化させない変形は、
【０１５５】
【数５４】
2005183411.tif 000056
【０１５６】
となる。
更に、式（１３）の左辺を、オーダーεで近似すると、
【０１５７】
【数５５】
2005183411.tif 000057
【０１５８】
となる。
これは、式（５）を離散化したものとなっており、陽解法により解を求めることができる。
【０１５９】
ある界面に対して、法線方向の変形に対応する速度分布｛Ｕｎ｝が与えられれば、Ｖ０を適当に定め、順次、
【０１６０】
【数５６】
2005183411.tif 000058
【０１６１】
を解いてゆくことにより｛Ｖｎ｝を定めることができる。
そして、｛Ｕｎ｝｛Ｖｎ｝と｛φｎ｝とによって新たな角度が、
【０１６２】
【数５７】
2005183411.tif 000059
【０１６３】
から定まり、ある時間ステップ後の形状が定まることとなる。
【０１６４】
但し、式（１４）から明らかなように、これらの式は、ｎ−１とｎ＋１とに入れ替えに対して不変ではない。すなわち、取り扱われる物理系が本来は有しないはずの左右非対称性が、離散化を施すことによって生じたこととなっている。
【０１６５】
以下の実施例では、この左右非対称性を生じさせない差分方程式の解法について説明する。
以下に説明する本発明の差分方程式の解法では、何らかの物理的条件により、界面に垂直な方向の速度Ｕｎが与えられる場合を想定している。そして、界面の移動をオイラー描像的に表現するために、与えられた速度Ｕｎに対して、式（１５）に基づく計算を逐次実行し、界面の接線方向の速度Ｖｎを得るのである。
【０１６６】
このような計算が意味をもつためには、任意の弧長ｓについて、正負充分大きい領域での界面の漸近的振舞いが均一でかつ既知である必要がある。
【０１６７】
図３は、着目した領域５２の領域から充分に離れた領域５１および５３では、Ｕｎが一定であり、かつ、曲率ｋｎがゼロであると仮定した場合の線分列の様子を示した図である。このような線分列の様子は、界面の変化を表現するものとして充分に一般的であり、半導体プロセス等において生起する界面の変化をシミュレーションする場合に、頻繁に生じる境界条件である。
【０１６８】
（実施例）
図４は、本発明の界面追跡シミュレーションを実行させるためのシミュレーション・システムの例を示す。このシミュレーション・システムは、シミュレーションの演算を司る計算機システムと、プログラムの実行ファイルを格納する記憶媒体とから構成される。計算機システムは、計算機本体であるコンピュータ１１１と、コンピュータ１１１に内蔵されたフロッピー（登録商標）ディスクドライブ１１３と、表示装置であるディスプレイ１１２と、入力装置であるマウス１１４およびキーボード１１５とから構成され、記憶媒体としては、フロッピー（登録商標）ディスク１１０が用いられる。
【０１６９】
図５は、本シミュレーション・システムを構成する、表示装置、入力装置、記憶装置、および演算装置の相互の関係を説明するブロック図である。
【０１７０】
表示装置であるディスプレイ１１２と、入力装置であるマウス１１４およびキーボード１１５は、コンピュータ１１１に内蔵された演算装置であるＣＰＵ１１６に接続されており、入力装置を介して入力された信号が入力されるとともに、演算結果をディスプレイ１１２に出力する。また、コンピュータ１１１に内蔵されたハードディスク１１７とメモリ１１８、および、フロッピー（登録商標）ディスクドライブ１１３を介して読込み・書出しがなされるフロッピー（登録商標）ディスク１１０が、ＣＰＵ１１６に接続され、演算に必要な情報がＣＰＵ１１６により読み込まれるとともに、演算結果が記憶されるようになっている。
【０１７１】
本発明の界面追跡アルゴリズムに基づいたプログラムの実行ファイルは、コンピュータ１１１に内蔵されたハードディスク１１７、あるいは外部記憶装置であるフロッピー（登録商標）ディスク１１０に予め記憶されている。界面追跡シミュレーション・プログラムの内容をＣＰＵ１１６により実行させるため、キーボード１１５あるいはマウス１１４からの指示入力によりメモリ１１８上にロードし、本発明の界面追跡シミュレーションが実行される。
【実施例１】
【０１７２】
本実施例では、本発明のシミュレーション方法における差分方程式の計算方法について説明する。
【０１７３】
従来の計算方法の問題点は、ある時間ステップ後に生じる界面形状が、左右非対称性を有することである。この問題を克服するために、本実施例の計算方法では、界面の接線方向の速度Ｖ０を、その非対称性を打ち消すように設定する。
【０１７４】
図６は、本発明の計算方法のアルゴリズムを示すフローチャートである。また、図７は、このアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子を示す図である。
【０１７５】
先ず、初期化を行い、ハードディスク１１７からメモリ１１８上に実行プログラムをロードし、必要な計算領域を確保したうえで、所望の方法によって入力された初期界面形状を形成する。また、界面の接線方向の速度の初期値Ｖ０をゼロに設定する（Ｓ２１）。
【０１７６】
次に、物理的考察に基づいて、界面の法線方向の、離散化された速度場Ｕｎを計算する（Ｓ２２）。例えば、実行されるべきシミュレーションが、酸化による界面の時間変化という現象に関するものであれば、酸化という物理現象を記述するための差分方程式により、速度場Ｕｎを計算する。一般的に、速度場は一形式の量であるので、離散化された場合も線分上で定義される。本実施例では、式（８）を考慮して、図１８に示すように位置ベクトルｒｎとｒｎ−１とで定まる線分上で定義された速度場をＵｎとした。これらの計算の結果として、｛Ｕｎ｜ｎ＝１，…，Ｎ｝が定まる。
【０１７７】
更に、Ｕ０を、Ｕ０＝Ｕ１として設定する（Ｓ２３）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０１７８】
引き続いて、与えられたＶ０の値に対して、式（１５）を順次解くことにより、｛Ｖｎ｝を定め（Ｓ２４）、更に、式（１６）に従って新しい接角分布｛φｎε｝を決定し、改めて｛φｎ｝と定義し直す（Ｓ２５）。
【０１７９】
これらのパラメータ設定を基に、界面形状の再構成をおこなう（Ｓ２６）。本実施例の場合には、ｒ０での時間発展が既知であるので、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、形状｛φｎ｝のデータから構成し直すことになる。
【０１８０】
そして、図７において、ｔ＋εの時点での漸近的な長さの差であるδｓを求め、次のステップでのＶ０の値を−δｓとしておく（Ｓ２７）。
【０１８１】
最後に、終了判定（Ｓ２８）を行い、終了と判定（Ｓ２８：Ｙｅｓ）されれば終了し（Ｓ２９）、継続と判定（Ｓ２８：Ｎｏ）されれば、Ｓ２２〜Ｓ２８のステップを繰り返す。このとき、Ｖ０＝−δｓと設定して計算を行うこととするので、図７に示すように左右の非対称性が起こらなくなる。
【実施例２】
【０１８２】
実施例１で示した差分スキームを含んだ界面追跡シミュレータ、あるいはサブルーチンを、図４に示したフロッピー（登録商標）ディスク１１０またはハードディスク１１７等の記憶手段に格納しておく。これにより、容易にシミュレーションができることとなる。
【実施例３】
【０１８３】
従来の計算方法により得られる計算結果に左右非対称性が現れるという問題点は、線分上で定義されるべき速度ベクトルを、頂点上で定義していることに起因している。この問題を克服するために、本実施例では、速度ベクトルの頂点への割り当てを変更することで、左右非対称性を打ち消すこととしている。
【０１８４】
図８は、本実施例の計算方法のアルゴリズムを示すフローチャートである。
【０１８５】
先ず、初期化を行い、ハードディスク１１７から実行プログラムをメモリ１１８上にロードし、必要な計算領域を確保した上で、所望の方法により入力された初期界面形状を形成する（Ｓ２２１）。
【０１８６】
次に、物理的考察に基づいて、界面の法線方向の、離散化された速度場Ｕｎを計算する（Ｓ２２２）。例えば、実行されるべきシミュレーションが酸化による界面の変化時間という現象についてのものであれば、酸化という物理現象を記述するための差分方程式により、速度場Ｕｎを計算する。一般的に、速度場は一形式の量であるので、離散化された場合も線分上で定義される。
【０１８７】
図１８は、式（８）を変形して、速度ベクトルの頂点への自然な割り当てを行うことにより、位置ベクトルｒｎとｒｎ−１の上で速度場Ｕｎを定義した様子を示す図である。この割り当てにより、｛Ｕｎ｜ｎ＝１，…，Ｎ｝の速度が定まる。
【０１８８】
更に、Ｕ０を、Ｕ０＝Ｕ１として設定する（Ｓ２２３）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０１８９】
引き続いて、与えられたＶ０の値に対して、式（１５）を順次解くことにより、｛Ｖｎ｝を定め（Ｓ２２４）、更に、式（１６）に従って新しい接角分布｛φｎε｝を決定し、改めて｛φｎ｝と定義し直す（Ｓ２２５）。
【０１９０】
これらのパラメータ設定を基に、界面形状の再構成をおこなう（Ｓ２２６）。本実施例の場合には、ｒ０での時間発展が既知であるので、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、形状｛φｎ｝のデータから構成し直すこととなる。
【０１９１】
また、Ｕｎの決定２として、再び、物理的考察により法線方向の離散化された速度場Ｕｎを計算する（Ｓ２２７）。決定２では図１８とは異なる図９に示す速度場Ｕｎの割り当て方を行う。つまり、位置ベクトルｒｎとｒｎ−１の上で定義された速度場Ｕｎを考えているから、界面形状は、
【０１９２】
【数５８】
2005183411.tif 000060
【０１９３】
を満足することとなる。この割り当てにより、｛Ｕｎ｜ｎ＝１，…，Ｎ−１｝の速度が定まる。
【０１９４】
次に、ＵＮを、ＵＮ＝ＵＮ−１として設定する（Ｓ２２８）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０１９５】
また、最初の速度であるＶ０をゼロとして、図９に示す線分列を構成する各線分の長さを保存する方程式を、式（１６）により計算し、その計算結果を還元した式、
【０１９６】
【数５９】
2005183411.tif 000061
【０１９７】
を解くことにより、｛Ｖｎ｝を定める（Ｓ２２９）。
更に、
【０１９８】
【数６０】
2005183411.tif 000062
【０１９９】
に従って新しい接角分布｛φｎε｝を決定し、その結果を改めて｛φｎ｝と定義し直し（Ｓ２３０）、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、｛φｎ｝のデータから界面形状を構成し直す（Ｓ２３１）。
【０２００】
最後に、終了判定（Ｓ２３２）を行い、終了と判定（Ｓ２３２：Ｙｅｓ）されれば終了し（Ｓ２３３）、継続と判定（Ｓ２３２：Ｎｏ）されれば、Ｓ２２２〜Ｓ２３２のステップを繰り返す。
【０２０１】
図１０は、上述したアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子である。この図から明らかなように、上述したアルゴリズムに従って差分方程式を解くことにより、関心のある領域における左右の非対称性が起こらなくなる。
【実施例４】
【０２０２】
実施例３で示した差分スキームを含んだ界面追跡シミュレータまたはサブルーチンを、フロッピー（登録商標）ディスク１１０やハードディスク１１７等の記憶媒体に収めることにより、手軽にシミュレーションができることとなる。
【実施例５】
【０２０３】
本実施例では、各ステップ毎に界面形状を表現する曲線の向きの入れ替えを行うことにより界面形状の左右非対称性を打ち消すようにした、本発明のシミュレーション方法における差分方程式の計算方法について説明する。具体的には、各ステップ毎に格子の順番を変更する取り扱いを行っている。
【０２０４】
図１４は、本発明の計算方法のアルゴリズムを示すフローチャートである。また、図１５は、このアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子を示す図である。
【０２０５】
先ず、初期化を行い、ハードディスク１１７からメモリ１１８上に実行プログラムをロードし、必要な計算領域を確保したうえで、所望の方法によって入力された初期界面形状を形成する（Ｓ３１）。
【０２０６】
次に、物理的考察に基づいて、界面の法線方向の、離散化された速度場Ｕｎを計算する（Ｓ３２）。例えば、実行されるべきシミュレーションが、酸化による界面の時間変化という現象に関するものであれば、酸化という物理現象を記述するための差分方程式により、速度場Ｕｎを計算する。一般的に、速度場は一形式の量であるので、離散化された場合も線分上で定義される。本実施例では、式（８）を考慮して、位置ベクトルｒｎとｒｎ−１とで定まる線分上で定義された速度場をＵｎとした。これらの計算の結果として、｛Ｕｎ｜ｎ＝１，…，Ｎ｝が定まる。
【０２０７】
更に、Ｕ０を、Ｕ０＝Ｕ１として設定する（Ｓ３３）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０２０８】
引き続いて、与えられたＶ０の値に対して、式（１５）をｎ値の小さい順に順次解くことにより、｛Ｖｎ｝を定め（Ｓ３４）、更に、式（１６）に従って新しい接角分布｛φｎε｝を決定し、改めて｛φｎ｝と定義し直す（Ｓ３５）。
【０２０９】
これらのパラメータ設定を基に、界面形状の再構成をおこなう（Ｓ３６）。本実施例の場合には、ｒ０での時間発展が既知であるので、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、形状｛φｎ｝のデータから構成し直すことになる。
【０２１０】
図１５に示すように、ｔ＋εの時点において、漸近的な長さの差δｓが存在する。本発明のシミュレーションに用いる差分方程式の計算では、図１５に示すように、頂点列の名前を反転させ、曲線の向きを逆転させる。また、この逆転操作に伴って、変数｛φｎ｝、｛Ｕｎ｝、｛Ｖｎ｝も添え字を逆転させる（Ｓ３７）。なお、プログラムの実装上は、それぞれの向き付け用のメモリ領域を取り、それらを入れ替えて、順番に使用してもよい。
【０２１１】
最後に、終了判定（Ｓ３８）を行い、終了と判定（Ｓ３８：Ｙｅｓ）されれば終了し（Ｓ３９）、継続と判定（Ｓ３８：Ｎｏ）されれば、Ｓ３２〜Ｓ３８のステップを繰り返す。
【０２１２】
上述した差分方程式の計算を行うことにより、左右の非対称性が起こらなくなる。
【実施例６】
【０２１３】
実施例５で示した差分スキームを含んだ界面追跡シミュレータまたはサブルーチンを、フロッピー（登録商標）ディスク１１０やハードディスク１１７等の記憶媒体に収めることにより、手軽にシミュレーションができることとなる。
【実施例７】
【０２１４】
従来の計算方法により得られる計算結果に左右非対称性が現れるという問題点は、線分上で定義されるべき速度ベクトルを、頂点上で定義していることに起因している。この問題を克服するために、本実施例では、速度ベクトルの頂点への割り当てを変更し、解く順序を変更することで、左右非対称性を打ち消すこととしている。
【０２１５】
図１６は、本実施例の計算方法のアルゴリズムを示すフローチャートである。
【０２１６】
先ず、初期化を行い、ハードディスク１１７から実行プログラムをメモリ１１８上にロードし、必要な計算領域を確保した上で、所望の方法により入力された初期界面形状を形成する（Ｓ３２１）。
【０２１７】
次に、物理的考察に基づいて、界面の法線方向の、離散化された速度場Ｕｎを計算する（Ｓ３２２）。例えば、実行されるべきシミュレーションが酸化による界面の変化時間という現象についてのものであれば、酸化という物理現象を記述するための差分方程式により、速度場Ｕｎを計算する。一般的に、速度場は一形式の量であるので、離散化された場合も線分上で定義される。
【０２１８】
図１８は、式（８）を変形して、速度ベクトルの頂点への自然な割り当てを行うことにより、位置ベクトルｒｎとｒｎ−１の上で速度場Ｕｎを定義した様子を示す図である。この割り当てにより、｛Ｕｎ｜ｎ＝１，…，Ｎ｝の速度が定まる。
【０２１９】
更に、Ｕ０を、Ｕ０＝Ｕ１として設定する（Ｓ３２３）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。引き続き、ＵＮをＵＮ＝ＵＮ−１として設定する。このことにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０２２０】
更に、Ｖ０の値をゼロとし、式（１５）をｎ値の小さい順に順次解くことにより、｛Ｖｎ｝を定め（Ｓ３２４）、更に、式（１６）に従って新しい接角分布｛φｎε｝を決定し、改めて｛φｎ｝と定義し直す（Ｓ３２５）。
【０２２１】
これらのパラメータ設定を基に、界面形状の再構成をおこなう（Ｓ３２６）。本実施例の場合には、ｒ０での時間発展が既知であるので、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、形状｛φｎ｝のデータから構成し直すこととなる。
【０２２２】
また、Ｕｎの決定２として、再び、物理的考察により法線方向の離散化された速度場Ｕｎを計算する（Ｓ３２７）。本実施例では、図９に示すように、位置ベクトルｒｎとｒｎ−１の上で定義された速度場Ｕｎを考えているから、界面形状は、
【０２２３】
【数６１】
2005183411.tif 000063
【０２２４】
を満足することとなる。この割り当てにより、｛Ｕｎ｜ｎ＝１，…，Ｎ−１｝の速度が定まる。
【０２２５】
次に、ＵＮを、ＵＮ＝ＵＮ−１として設定する（Ｓ３２８）。これにより、｛Ｕｎ｜ｎ＝０，…，Ｎ｝が設定される。
【０２２６】
また、ＶＮをゼロとして、図９に示す線分列を構成する各線分の長さを保存する方程式を、式（１６）により計算し、その計算結果を還元した式、
【０２２７】
【数６２】
2005183411.tif 000064
【０２２８】
をｎ値の大きいほうから順次解くことにより、｛Ｖｎ｝を定める（Ｓ３２９）。
更に、
【０２２９】
【数６３】
2005183411.tif 000065
【０２３０】
に従って新しい接角分布｛φｎε｝を決定し、その結果を改めて｛φｎ｝と定義し直し（Ｓ３３０）、ｒ０＋εＵ０ｎ０＋εＶ０ｔ０として、｛φｎ｝のデータから界面形状を構成し直す（Ｓ３３１）。これは実質的には、実施例５で示した方向の逆転と同じことを行ったこととなる。
【０２３１】
最後に、終了判定（Ｓ３３２）を行い、終了と判定（Ｓ３３２：Ｙｅｓ）されれば終了し（Ｓ３３３）、継続と判定（Ｓ３３２：Ｎｏ）されれば、Ｓ３２２〜Ｓ３３２のステップを繰り返す。
【０２３２】
図１７は、上述したアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子である。この図から明らかなように、上述したアルゴリズムに従って差分方程式を解くことにより、関心のある領域における左右の非対称性が起こらなくなる。
【０２３３】
上述した差分方程式の計算を行うことにより、左右の非対称性が起こらなくなる。
【実施例８】
【０２３４】
実施例７で示した差分スキームを含んだ界面追跡シミュレータまたはサブルーチンを、フロッピー（登録商標）ディスク１１０やハードディスク１１７等の記憶媒体に収めることにより、手軽にシミュレーションができることとなる。
【産業上の利用可能性】
【０２３５】
界面形状を表現する線分列の各線分の頂点の順番が入れ替わるためにループを作ってしまう現象が現れるという不都合を生じないシミュレーション方法の提供が可能となる。
【図面の簡単な説明】
【０２３６】
【図１】滑らかな曲線についての古典幾何学の説明図である。
【図２】１本の曲線に離散化を施して得られた線分列を説明する図である。
【図３】着目した領域から充分離れた領域ではＵｎが一定で、かつ、曲率ｋｎがゼロである界面を説明する図である。
【図４】本発明の界面追跡方法を実行させるシミュレーション手段を説明するための図である。
【図５】本発明のシミュレーション手段を説明するためのブロック図である。
【図６】本発明の計算方法のアルゴリズムを説明するためのフローチャートである。
【図７】本発明の計算方法のアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子を説明するための図である。
【図８】本発明の計算方法のアルゴリズムを説明するためのフローチャートである。
【図９】速度ベクトルの頂点への自然な割り当てを行うことにより、位置ベクトルｒｎとｒｎ−１の上で速度場Ｕｎを定義した様子を示す図である。
【図１０】本発明の計算方法のアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子である。
【図１１】長さの異なる複数の線分からなる線分列によって界面形状を表現したものである。
【図１２】界面形状の経時変化の様子を説明するための図である。
【図１３】従来のシミュレーション方法で生じる問題点を説明するための図である。
【図１４】本発明の計算方法のアルゴリズムを説明するためのフローチャートである。
【図１５】本発明の計算方法のアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子である。
【図１６】本発明の計算方法のアルゴリズムを説明するためのフローチャートである。
【図１７】本発明の計算方法のアルゴリズムに従って差分方程式を解いて得られた界面形状の時間変化の様子である。
【図１８】速度ベクトルの頂点への自然な割り当てを行うことにより、位置ベクトルｒｎとｒｎ−１の上で速度場Ｕｎを定義した様子を示す図である。
【符号の説明】
【０２３７】
５１  着目する領域から離れた領域
５２  着目する領域
５３  着目する領域から離れた領域
１１０  フロッピー（登録商標）ディスク
１１１  コンピュータ
１１２  ディスプレイ
１１３  フロッピー（登録商標）ディスクドライブ
１１４  マウス
１１５  キーボード
１１６  ＣＰＵ
１１７  ハードディスク
１１８  メモリ
【図１】
2005183411.tif 000106
【図２】
2005183411.tif 000107
【図３】
2005183411.tif 000108
【図４】
2005183411.tif 000109
【図５】
2005183411.tif 000110
【図６】
2005183411.tif 000111
【図７】
2005183411.tif 000112
【図８】
2005183411.tif 000113
【図９】
2005183411.tif 000114
【図１０】
2005183411.tif 000115
【図１１】
2005183411.tif 000116
【図１２】
2005183411.tif 000117
【図１３】
2005183411.tif 000118
【図１４】
2005183411.tif 000119
【図１５】
2005183411.tif 000120
【図１６】
2005183411.tif 000121
【図１７】
2005183411.tif 000122
【図１８】
2005183411.tif 000123