(57)【要約】【課題】  より物理的な接触角を表現する技術を提供すること。【解決手段】  空間内各点においてＦ１＋Ｆ２＋Ｆ３＝１を満たす関数Ｆ１とＦ２とＦ３をそれぞれ３種類の物質の体積占有率とし、勾配の絶対値｜∇Ｆ１｜、｜∇Ｆ２｜、｜∇Ｆ３｜と３つの実定数σ１、σ２、σ３に対して、エネルギー積分を微小体積要素ｄｖｏｌを利用して、Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ＋σ３∫｜∇Ｆ３｜ｄｖｏｌと定義し、このエネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式により３つの物質の表面張力を近似的に表現する。【選択図】  図１
2006284352.tif 000002
【特許請求の範囲】
【請求項１】
２次元以上の空間内の３種類の物質の界面張力を、各物質に関する体積占有率を表わす関数と所定の定数とを記憶したデータ記憶部を備えた情報処理装置により計算する情報処理方法において、
エネルギー積分を、前記各物質について、前記関数の勾配の絶対値の体積積分に前記所定の定数を乗じた値の和と定義し、
前記データ記憶部に記憶された前記関数及び定数を用いて、前記エネルギー積分を空間方向に対する汎関数微分により得られる微分方程式を差分化した式により、前記３種類の物質の表面張力を計算することを特徴とする情報処理方法。
【請求項２】
考慮している領域全体で定義されたゼロから１の間の値を有し、空間内各点においてＦ１＋Ｆ２＋Ｆ３＝１を満たす関数Ｆ１とＦ２とＦ３（＝１−Ｆ１−Ｆ２）をそれぞれ前記３種類の物質の体積占有率を表わす関数とし、前記勾配の絶対値を｜∇Ｆ１｜、｜∇Ｆ２｜、｜∇Ｆ３｜とし、前記所定の定数をσ１、σ２、σ３として、前記エネルギー積分を微小体積要素ｄｖｏｌを利用して、
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ＋σ３∫｜∇Ｆ３｜ｄｖｏｌ
と定義することを特徴とする請求項１に記載の情報処理方法。
【請求項３】
特に考慮している領域全体で定義されたゼロから１の間の値を有する定関数Ｇを用いて、各点でＦ１＋Ｆ２＝Ｇとなる場合、すなわち変分変数がＦ１またはＦ２のいずれかひとつに限られる場合において、勾配の絶対値｜∇Ｆ１｜、｜∇Ｆ２｜、２つの実定数σ１、σ２に対して、エネルギー積分を微小体積要素ｄｖｏｌを利用して、
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ
と定義し、当該エネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式により３つの物質の表面張力を計算することを特徴とする請求項２に記載の情報処理方法。
【請求項４】
３つの実定数σ１、σ２、σ３が、物質１と物質２との界面の単位面積あたりの界面エネルギーσ１２、物質２と物質３との界面の単位面積あたりの界面エネルギーσ２３、物質３と物質１との界面の単位面積あたりの界面エネルギーσ３１を用いて、
σ１＝（  σ１２−σ２３＋σ３１）／２
σ２＝（  σ１２＋σ２３−σ３１）／２
σ３＝（−σ１２＋σ２３＋σ３１）／２
とあらわされることを特徴とする請求項２に記載の情報処理方法。
【請求項５】
２つの実定数σ１、σ２が、物質１と物質２との界面の単位面積あたりの界面エネルギーσ１２、物質２と物質３との界面の単位面積あたりの界面エネルギーσ２３、物質３と物質１との界面の単位面積あたりの界面エネルギーσ３１を用いて、
σ１＝（σ１２−σ２３＋σ３１）／２
σ２＝（σ１２＋σ２３−σ３１）／２
とあらわされることを特徴とする請求項３に記載の情報処理方法。
【請求項６】
２次元以上の空間を占める３つの物質の界面張力を近似的に計算する計算方法において、Ｆ１＋Ｆ２＋Ｆ３＝１となるように、前記空間の各点における物質１の体積占有率を表す関数Ｆ１のデータ、物質２の体積占有率を表す関数Ｆ２のデータ、物質３の体積占有率を表す関数Ｆ３のデータを保持する保持手段と、
エネルギー積分Ｅを｜∇Ｆ１｜、｜∇Ｆ２｜、｜∇Ｆ３｜と３つの実定数σ１、σ２、σ３と微小体積要素dvolとを用いて、
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ＋σ３∫｜∇Ｆ３｜ｄｖｏｌと定義した場合に、当該エネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式に従って、前記保持手段が保持するデータを用いて計算処理を行うことで、前記３つの物質の表面張力を近似的に計算する計算手段と
を備えることを特徴とする情報処理装置。
【請求項７】
コンピュータに請求項１乃至５の何れか１項に記載の情報処理方法を実行させることを特徴とするプログラム。
【請求項８】
請求項７に記載のプログラムを格納したことを特徴とする、コンピュータ読み取り可能な記憶媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、界面張力の計算技術に関するものである。
【背景技術】
【０００２】
流体力学によって記述される自然現象を計算機上でより正確に模倣し、予測する事を目指す学問分野として計算流体力学と言われる分野がある。この計算流体力学には、流体における表面張力の計算も含まれている（特許文献１を参照」）。
【０００３】
非特許文献１、２は、広い意味での流体体積法（VOF法と以下呼ぶ）という方法を利用したものであると考えてよい。領域全体で定義された状態の占有率関数F(X)を考える。占有率関数F(X)は物質１と物質２のみが存在する領域において、物質１の占有率を表現したものである。この関数に対して、エネルギー汎関数として、
Ｅ１２＝σ１２∫｜∇Ｆ｜ｄｖｏｌ      （１）
を考えることで、表面張力を表現できることを上記非特許文献１、２の著者らは発見した。σ１２は、物質１と物質２の単位面積あたりの界面エネルギーである。空間をデカルト座標（ｘ１，ｘ２，ｘ３）で表すこととする。上記式（１）の汎関数に対して、ネーターの定理を適応して得られるi方向の運動量密度は、
ｆｉ＝∇ｉ（｜∇Ｆ｜）−Σｊ∇ｊ（（∇ｊＦ）（∇ｉＦ）／｜∇Ｆ｜）    （２）
となる。ネーターの定理に関しては例えば、非特許文献３に書かれている。スカラー場の場合に相当することを注意すると、上記の式（２）が得られる。これが所謂表面張力である。Ｅ１２は物質１と物質２の界面の表面エネルギーを近似的に表現していることは重要なことであり、表面エネルギーは界面の面積に比例することは重要なことである。即ち、表面張力は表面積を最小にする最小原理によって作られた力である。
【０００４】
他方、３種類の物質が混在する系において上記の考えを拡張することが上記文献に述べられている。３種類の物質が存在する系においては、図９に示すような接触角により、３つの界面が交わる三重点が特徴付けられることが知られている（非特許文献４を参照）
この接触角を表現する目的に対しても上記非特許文献１は第三の物質（以下物質３）が固体である場合に関して、計算方法を与えている。それによると、物質１と物質２の界面が、物質３に交わる接触角を幾何学的に求めて、それを所望の角度に変更することにより３重点の構成を行っている。
【特許文献１】特開平６−３０９２９４号公報
【非特許文献１】J.U.Backbill, D.B.Kothe, C.Zemach  Journal of Computational Physics 100巻  335-354 1992年  ”A continuum Method for Modeling Surface Tension”
【非特許文献２】D.B.Koth, R.C.Mjolsness, M.D.Torrey  Los  Alamos National Laboratory (LANL)  LANL Report Library ((505) 667-4446)  ”RIPPLE: A Computer Program for Imcompressible Flows with Free Surfaces”
【非特許文献３】「量子場を学ぶための場の解析力学入門」高橋康著  講談社サイエンティフック  １９８２年発行  のp.101(7.10)及び(7.3)
【非特許文献４】G.K.バチェラー著、橋本英典，松信八十男他訳の「入門流体力学」東京電気大学出版局  １９７４年発行
【発明の開示】
【発明が解決しようとする課題】
【０００５】
本発明の課題は、計算流体力学において、接触角を如何に表現するかにある。従来技術  述べた方法は、接触角を幾何学的に捕らえるものである。しかしながら、上記非特許文献１に書れているように、ダイナミックスが絡んできた場合等、方法には任意性を伴うものである。平衡点においては、接触角は表面張力係数のみで決定されるが、ダイナミックスを追いかけている際には、接触角に物理的な意味はない。
【０００６】
本発明は以上の問題に鑑みてなされたものであり、より物理的な接触角を表現する技術を提供することを目的とする。
【課題を解決するための手段】
【０００７】
本発明の目的を達成するために、例えば、本発明の情報処理方法は以下の構成を備える。
【０００８】
即ち、２次元以上の空間内の３種類の物質の界面張力を、各物質に関する体積占有率を表わす関数と所定の定数とを記憶したデータ記憶部を備えた情報処理装置により計算する情報処理方法において、
エネルギー積分を、前記各物質について、前記関数の勾配の絶対値の体積積分に前記所定の定数を乗じた値の和と定義し、
前記データ記憶部に記憶された前記関数及び定数を用いて、前記エネルギー積分を空間方向に対する汎関数微分により得られる微分方程式を差分化した式により、前記３種類の物質の表面張力を計算することを特徴とする。
【０００９】
本発明の目的を達成するために、例えば、本発明の情報処理装置は以下の構成を備える。
【００１０】
即ち、２次元以上の空間を占める３つの物質の界面張力を近似的に計算する処理を行う情報処理装置であって、
Ｆ１＋Ｆ２＋Ｆ３＝１となるように、前記空間の各点における物質１の体積占有率を表す関数Ｆ１のデータ、物質２の体積占有率を表す関数Ｆ２のデータ、物質３の体積占有率を表す関数Ｆ３のデータを保持する保持手段と、
エネルギー積分Ｅを｜∇Ｆ１｜、｜∇Ｆ２｜、｜∇Ｆ３｜と３つの実定数σ１、σ２、σ３と微小体積要素dvolとを用いて、
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ＋σ３∫｜∇Ｆ３｜ｄｖｏｌと定義した場合に、当該エネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式に従って、前記保持手段が保持するデータを用いて計算処理を行うことで、前記３つの物質の表面張力を近似的に計算する計算手段と
を備えることを特徴とする。
【発明の効果】
【００１１】
本発明の構成により、より物理的な接触角を表現することができる。
【発明を実施するための最良の形態】
【００１２】
以下添付図面を参照して、本発明を好適な実施形態に従って詳細に説明する。
【００１３】
＜情報処理装置の構成について＞
図７は、以下の実施形態に係る情報処理装置として機能するコンピュータの外観を示す図である。同図に示す如く、このコンピュータは本体１１１、ディスプレイ１１２，キーボード１１５、マウス１１４により構成されており、以下説明する各処理を本体１１１内に搭載されているＣＰＵに実行させるためのプログラムやデータが格納されている記憶媒体（同図ではフロッピー（登録商標）ディスク）１１０は周知の通り、本体１１１に設けられているフロッピー（登録商標）ディスクドライブ装置１１３に挿入するもので、挿入することで、この記憶媒体１１０に記憶されているプログラムやデータは本体１１１内に搭載されているメモリにロードされることになる。
【００１４】
図８は、上記コンピュータの基本構成を示すブロック図である。同図において図７と同じ部分には同じ番号を付けている。
【００１５】
１１２は上述のディスプレイで、ＣＲＴや液晶画面などにより構成されており、ＣＰＵ１１６による処理結果を画像や文字などでもって表示することができる。
【００１６】
１１４、１１５はそれぞれマウス、キーボードで、各種の指示をＣＰＵ１１６に対して入力することができる。
【００１７】
１１６はＣＰＵで、メモリ１１８に記憶されているプログラムやデータを用いてコンピュータ全体の制御を行うと共に、後述する各計算処理を実行する。
【００１８】
１１８はメモリで、ハードディスクドライブ装置１１７やフロッピー（登録商標）ディスクドライブ装置１１３からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ１１６が各種の処理を実行する際に用いるワークエリアを備える。
【００１９】
１１７はハードディスクドライブ装置で、ここにＯＳ（オペレーティングシステム）や、ＣＰＵ１１６に後述する各計算処理を実行させるために必要なプログラムやデータが保存されており、これらの一部もしくは全部はＣＰＵ１１６の制御に従ってメモリ１１８にロードされ、ＣＰＵ１１６の処理対象となる。
【００２０】
１１３はフロッピー（登録商標）ディスクドライブ装置で、図７に示したフロッピー（登録商標）ディスク１１０に記録されているプログラムやデータを読み取る装置であって、ＣＰＵ１１６の制御により動作する。フロッピー（登録商標）ディスクドライブ装置１１３が読み取ったプログラムやデータはメモリ１１８にロードされ、ＣＰＵ１１６の処理対象となる。
【００２１】
以上の構成を有するコンピュータにより、以下説明する各計算処理を実行するのであるが、その実行結果は単に数値の文字列としてディスプレイ１１２に表示しても良いし、計算処理によって得られた数値列に基づいたコンピュータグラフィックスでもってディスプレイ１１２に表示するようにしても良い。
【００２２】
＜計算アルゴリズム＞
以下説明する表面張力の計算アルゴリズムは以下のエネルギー汎関数を基礎としている。
【００２３】
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ＋σ３∫｜∇Ｆ３｜ｄｖｏｌ    （３）
但し、σ１、σ２、σ３は、表面張力係数（正の実数）であり、dvolはルベーグ測度から決まる微小体積要素である。Ｆ１、Ｆ２及びＦ３は［０，１］の値を持つ関数で、全領域で定義されており、各点でＦ１＋Ｆ２＋Ｆ３＝１を満たす。また、Ｆ１は物質１の占有率を表現しているとし、Ｆ１が１の時、その場所は物質１のみが存在しているとする。Ｆ２は物質２の占有率を表現し、Ｆ２が１の時、物質２のみがその点で占有しているとする。Ｆ３は物質３の占有率を表現しており、Ｆ３が１のとき、物質３のみが占有しているとする。また、それぞれが、０と１の中間の数値を取る場合は、界面を近似的に表現しているとする。よって、｜∇Ｆ１｜、｜∇Ｆ２｜及び｜∇Ｆ３｜の台（非零値を取る定義域）は空間的に局在化している。
【００２４】
ここで、Ｆ３が０の領域における物資１と物質２の界面を考える。この場合、Ｆ１＋Ｆ２＝１であるので、
∫｜∇Ｆ２｜ｄｖｏｌ＝∫｜∇（１−Ｆ１）｜ｄｖｏｌ＝∫｜∇Ｆ１｜ｄｖｏｌ    （４）
である。従って、物質１と物質２との間の界面エネルギーは、物質１と物質２との界面の単位面積あたりの界面エネルギーσ１２を用いて、
Ｅ１２＝σ１２∫｜∇Ｆ１｜ｄｖｏｌ
＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ    （５）
＝（σ１＋σ２）∫｜∇Ｆ１｜ｄｖｏｌ
となるので、
σ１＋σ２＝σ１２    （６）
であればよい。同様にして、Ｆ１＝０における物質２と物質３との界面、Ｆ２＝０における物質３と物質１との界面について同様の考察を加えることにより、結局、σ１、σ２、σ３は、
σ１＋σ２＝σ１２
σ２＋σ３＝σ２３      （７）
σ３＋σ１＝σ３１
を満たす。すなわち、
σ１＝（  σ１２−σ２３＋σ３１）／２
σ２＝（  σ１２＋σ２３−σ３１）／２      （８）
σ３＝（−σ１２＋σ２３＋σ３１）／２
と設定すればよいことが分る。ここで、物質１と物質との界面の単位面積あたりの界面エネルギーはσ１２、物質２と物質３との界面の単位面積あたりの界面エネルギーはσ２３、物質３と物質１との界面の単位面積あたりの界面エネルギーはσ３１である。
【００２５】
さらに、何らかの理由で上記エネルギー汎関数がＦ３によらない場合（たとえば物質３は固体であって表面形状は変形し得ない、あるいは拘束されている場合）は、Ｆ３は変分変数ではないので、上記エネルギー汎積分は
Ｅ＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＋σ２∫｜∇Ｆ２｜ｄｖｏｌ      （９）
となり、変分変数はＦ１あるいはＦ２（＝１−Ｆ３−Ｆ１）のいずれかのひとつになる。
【００２６】
ここで「表面張力の名前」について説明する。表面張力を何らかの引っ張りによる力として捕らえると理解が歪んでしまう。例えば、古典幾何光学において、フェルマーの最小原理により、真空中の２点間を結ぶ光線は直線をとなるわけであるが、現在、これを我々は張力が光線（光の軌道）にあるいは光に働いたとは理解せず、光路長を最小にするものと理解している。しかし、歴史的には１６−１７世紀光路に張力が働いているという認識を持ち、弾性曲線と光の問題についての考察がなされていたようである。これを単純な最小原理に置き換えたのは、クラインによるとハミルトンの功績によるものである。（クライン：１９世紀の数学、Felix Klein著 弥永晶吉 監修 足立恒雄・浪川幸彦・監訳 石井省吾・渡辺 弘 訳 共立出版1995年発行。）物理的な理由として、なぜフェルマーの原理において、最小光路長という概念が張力という概念より自然であるかという理由は量子力学の立場から理解することによって明確になる。（ファインマン経路積分 L.S.シュルマン 著 高塚和夫 訳 講談社 95年発行に記述されている。）
同様に、ニュートン力学は、オイラー・ラグランジの最小原理によって、その軌道が決まるわけであるが、それも粒子の運動が何らかの張力が働くから、ニュートン力学に従うとは考えない。もちろん、歴史的にはオイラーは、時間方向の点粒子の動きを時間軸を空間内の弾性棒に置き換えて、張力の影響として理解していたようであるが、我々は慣性力という場合はあるものの、張力によるものとは理解しない。
【００２７】
そのような対比の下で、表面張力の「張力」という言葉は歴史的なものとして理解する必要がある。表面張力はある表面エネルギーを最小化する際に生じる力であり、弾性体のようなフックの法則に従う意味での張力ではない。
【００２８】
本発明の思想としては、最小原理という現代的な理解に従うことで、計算流体力学の課題を解決しようとするものである。
【００２９】
図１は、以上の説明に従った表面張力の計算アルゴリズムを用いた、流体の状態計算処理のフローチャートである。尚、同図のフローチャートに従った処理をＣＰＵ１１６に実行させるためのプログラムやデータはメモリ１１８にロードされており、これをＣＰＵ１１６が用いて処理を実行することで、上記コンピュータは以下説明する各処理を実行することになる。
【００３０】
先ず、以下の処理で用いる各変数を初期化する（ステップＳ１）。次に、予め定められた占有率関数Ｆ、Ｇに対して、以下説明する各計算処理によって表面張力の計算を行う（ステップＳ２）。そして、次に、予め与えられた適当な流体方程式を解くことで上記関数Ｆ、Ｇの物理状態を更新し（ステップＳ３）、その後、時間を更新させる（ステップＳ４）。
【００３１】
そして終了指示がない場合には処理をステップＳ１に戻し、以降の処理を繰り返す。
【００３２】
［第１の実施形態］
本実施形態では、２次元平面内で上記表面張力の計算を行う場合について説明する。詳しくは、２次元平面内に３つの異なる物質が存在する場合の表面張力の計算アルゴリズムについて説明する。まずは上記式（３）のエネルギー汎関数の第一項のみを考える。第一項は
Ｅ１＝σ１∫｜∇Ｆ１｜ｄｖｏｌ＝：σ１∫ｅ１ｄｖｏｌ      （１０）
と書かれる。σ１は、物質１と物質２との界面の単位面積あたりの界面エネルギーσ１２、物質２と物質３との界面の単位面積あたりの界面エネルギーσ２３、物質３と物質１との界面の単位面積あたりの界面エネルギーσ３１を用いて、σ１＝（σ１２−σ２３＋σ３１）／２で与えられる量である。
【００３３】
以下ではこのＥ１についてのみ説明するが、上記式（３）の第２項目、第３項目のそれぞれについても同様の処理を行えばよい。
【００３４】
この汎関数に対して、ネーターの定理を適応して得られるｉ方向の運動量密度は、
ｆｉ＝∇ｉ（｜∇Ｆ１｜）−Σｊ∇ｊ（（∇ｊＦ１）（∇ｉＦ１）／｜∇Ｆ１｜）      （１１）
より、
ｆｉ＝Σｊ∇ｉ（（∇ｊＦ１）（∇ｊＦ１）／｜∇Ｆ１｜）
−Σｊ∇ｊ（（∇ｊＦ１）（∇ｉＦ１）／｜∇Ｆ１｜）      （１２）
となる。よって、流体の慣用に従って、
τｉｉ＝Σｊ≠ｉ（（∇ｊＦ１）（∇ｊＦ１）／｜∇Ｆ１｜）
τｉｊ＝−（（∇ｉＦ１）（∇ｊＦ１）／｜∇Ｆ１｜）      （１３）
と定義する。
【００３５】
ｆｉ  ＝  ∇ｉτｉｉ−Σｊ∇ｊτｉｊ      （１４）
を計算すればよい。本実施形態では、これをスタッガード格子を利用して表現した。
【００３６】
まず、スカラー関数Ｆは、２次元内の正方格子の胞に定義されている。胞は図３に示すように整数（ｎ１，ｎ２）で指定されている。それらの集合を［Ｆ］と書き、スタッガード格子内のスカラー場あるいは単にスカラー場と呼ぶ。
【００３７】
また、胞の周りの辺ｘ１方向のＬ１とｘ２方向のＬ２によって構成され、それぞれの線分Ｌ１，Ｌ２も、整数（ｎ１，ｎ２）によって指定される。Ｌ１とＬ２は図４及び図５に太線で示している。各Ｌ１にはｘ２方向のベクトル場Ｖ２が、Ｌ２にはｘ１方向のベクトル場Ｖ２が図６のように離散化されて定義されている。それらの集合を［Ｖ１，Ｖ２］と記し、スタッガード格子内のベクトル場あるいは単にベクトル場と呼ぶ。
【００３８】
このとき、離散化されたアトレッド(atled)演算子∇は、スカラー場に対してはベクトル場の空間への写像
∇：  ［Ｆ］  →  ［Ｖ１，Ｖ２］      （１５）
として、また、その像は、
Ｖ１（ｎ１，ｎ２）：＝Ｆ（ｎ１＋１，ｎ２）−Ｆ（ｎ１，ｎ２）
Ｖ２（ｎ１，ｎ２）：＝Ｆ（ｎ１，ｎ２＋１）−Ｆ（ｎ１，ｎ２）      （１６）
となるように定義される。この写像の像［Ｖ１，Ｖ２］を［∇Ｆ］と記す。また［Ｖ１，  Ｖ２］から［Ｖ１］及び［Ｖ２］への自然な射影を利用して、［∇１Ｆ］及び［∇２Ｆ］をそれぞれの成分の集合として定義する。同様に、ベクトル場［Ｖ１，Ｖ２］に対するアトレッド(atled)演算子∇を、スカラー場の空間への写像、
∇：［Ｖ１，Ｖ２］→［Ｆ］      （１７）
として、またその像を
Ｆ（ｎ１，ｎ２）：＝Ｖ１（ｎ１，ｎ２）−Ｖ１（ｎ１−１，ｎ２）＋Ｖ２（ｎ１，ｎ２）−Ｖ２（ｎ１，ｎ２−２）      （１８）
と定義する。
【００３９】
また、このベクトル場［Ｖ１，Ｖ２］に対する絶対値｜｜は
｜｜：    ［Ｖ１，Ｖ２］→［Ｆ］      （１９）
となる写像として定義される。但し、
Ｆ（ｎ１，ｎ２）＝（［Ｖ１（ｎ１＋１，ｎ２）＋Ｖ１（ｎ１，ｎ２）］２／４＋［Ｖ２（ｎ１，ｎ２＋１）＋Ｖ２（ｎ１，ｎ２）］２／４）１／２      （２０）
とする。
【００４０】
二つのスカラー場［Ｆ］［Ｇ］の四則演算は、例えば掛け算が
［Ｆ］［Ｇ］（ｎ１，ｎ２）＝［Ｆ］（ｎ１，ｎ２）［Ｇ］（ｎ１，ｎ２）    （２１）
として各定義点での値の掛け算とするように各点の値の四則演算として定義する。
【００４１】
以上の離散化された場において、エネルギー積分は胞上の和として、一意に
Ｅ１＝Σｎ１，ｎ２（｜［∇Ｆ１］｜）（ｎ１，ｎ２）    （２２）
のように定義される。これをスタッガード格子上のエネルギーと呼ぶ。
【００４２】
同時に、ｘ成分摂動を加えた際のネーターの定理を適応した式（１１）の第一項である連続量としての∇ｉ（｜∇Ｆ１｜）は、離散量∇ｉ（｜∇［Ｆ１］｜）としてスタガード格子上で一意に定義される。つまり、これをネーターの定理を適応した式（１１）の第一項のスタッガード格子上での変形物として考える。
【００４３】
ｘ成分摂動を加えた際のネーターの定理を適応した式（１１）の第二項に関しては、ネーターの定理を適用した際の項
−Σｊ∇ｊ（（∇ｉＦ１）δＥ／δ（∇ｊＦ１））      （２３）
から得られる項である。  ここで、この項は正確には
−Σｊ∇ｊ（（∇ｉＦ１）＊δＥ／δ（∇ｊＦ１））      （２４）
となるテンソル積＊が陰に定義されている。そこで、テンソル積の各部分について説明する。まず
δＥ／δ（∇ｊＦ１）      （２５）
を考える。簡単のためにＬ２上の辺の場合のみ考察する。（ｎ１，ｎ２）で指定されるＬ２上では、
（Ｆ１（ｎ１＋２，ｎ２）−Ｆ１（ｎ１，ｎ２））／｜∇Ｆ１｜（ｎ１＋１，ｎ２）＋（Ｆ１（ｎ１＋１，ｎ２）−Ｆ１（ｎ１−１，ｎ２））／｜∇Ｆ１｜（ｎ１，ｎ２）      （２６）
となる。これと同じ定義点を持つ∇１Ｆ１とのテンソル積は通常の意味で定義する。他方、これと直行する［∇２Ｆ１］に関しては、定義点が異なるので、近傍の平均を取ったものをＬ２上のベクトルとしてテンソルとする。つまり、Ｌ２上で定義されるＶ２として近傍の平均を
Ｖ２：＝（Ｖ２（ｎ１，ｎ２）＋Ｖ２（ｎ１＋１，ｎ２）＋Ｖ２（ｎ１，ｎ２−１）＋Ｖ２（ｎ＋１，ｎ２−１））／４      （２７）
と定義する。これにより、（Ｆ１（ｎ１＋１，ｎ２）−Ｆ１（ｎ１−１，ｎ２））／｜∇Ｆ１｜（ｎ１，ｎ２）は［∇２Ｆ１］とテンソル積＊を取れる。これらを通常の意味で（実ベクトル空間として）のテンソル積＊とする。
【００４４】
これにより、［Ｖ１，Ｖ２］に対するアトレッド(atled)演算子∇を適用でき、上記の演算は定義できる。但し、計算した値は胞上で定義されるので、胞から各辺Ｌ１、Ｌ２に近傍の平均化操作で定義点を変えることができ、計算が可能となる。
【００４５】
このようにして、上記式（１０）に対する表面張力の式（１１）及び式（１４）がスタッガード格子上で離散化されて計算できる。式（３）の他の二項についても同様に計算でき、式（３）に対するネーターの定理を利用した運動量場としての表面張力がスタガード格子上で計算でき、離散空間上の計算流体力学の計算スキームの中に取り入れられる。
【００４６】
図２は上記アルゴリズムを利用して計算した表面エネルギーを考慮した形状を示す図である。同図ではＦ及びＧの等値線が描かれている。
【００４７】
［第２の実施形態］
第１の実施形態では説明を簡単にするために２次元平面上における計算アルゴリズムについて説明したが、この計算アルゴリズムはそのまま３次元空間に拡張しても良い。
【００４８】
［第３の実施形態］
上記汎関数を利用すれば、Ｆ，Ｇを有限の関数基底で表現し、該関数基底で張られるベクトル空間内の最小問題として、表面張力を捕らえることが可能となる。よって、有限要素法等々に適応できる。
【００４９】
以上の各実施形態における説明から、三相界面の問題をエネルギー的な視点から捉えることができ、非定常な場合に関しても、表面張力がエネルギー論的な意味から正確に計算できるようになった。これにより、エネルギー的観点から、問題点を克服できることとなった。
【００５０】
［その他の実施形態］
また、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記録媒体（または記憶媒体）を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記録媒体に格納されたプログラムコードを読み出し実行することによっても、達成されることは言うまでもない。この場合、記録媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記録した記録媒体は本発明を構成することになる。
【００５１】
また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているオペレーティングシステム（ＯＳ）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【００５２】
さらに、記録媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張カードやコンピュータに接続された機能拡張ユニットに備わるメモリに書込まれた後、そのプログラムコードの指示に基づき、その機能拡張カードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【００５３】
本発明を上記記録媒体に適用する場合、その記録媒体には、先に説明したフローチャートに対応するプログラムコードが格納されることになる。
【図面の簡単な説明】
【００５４】
【図１】表面張力の計算アルゴリズムを用いた、流体の状態計算処理のフローチャートである。
【図２】本発明の実施形態に係るアルゴリズムを利用して計算した表面エネルギーを考慮した形状を示す図である。
【図３】計算領域の胞について示す図である。
【図４】計算領域の辺Ｌ２について示す図である。
【図５】計算領域の辺Ｌ１について示す図である。
【図６】計算領域の辺Ｌ１、辺Ｌ２で定義されるベクトルについて示す図である。
【図７】本発明の実施形態に係る情報処理装置として機能するコンピュータの外観を示す図である。
【図８】コンピュータの基本構成を示すブロック図である。
【図９】３種類の物質が存在する系における接触角を説明する図である。
【図１】
2006284352.tif 000003
【図２】
2006284352.tif 000004
【図３】
2006284352.tif 000005
【図４】
2006284352.tif 000006
【図５】
2006284352.tif 000007
【図６】
2006284352.tif 000008
【図７】
2006284352.tif 000009
【図８】
2006284352.tif 000010
【図９】
2006284352.tif 000011