(57)【要約】【課題】  計算流体力学において、接触角を表現する。【解決手段】  インクジェット等の微小流体素子を数値的に模擬する際の２次元以上の空間内の３種類の物質の界面張力を近似的に計算する計算方法において、エネルギー積分を定義し、その定義した積分を離散化した後に、変分を取ることでエネルギー論的に無矛盾な離散空間上の表面張力を導出することを特徴とする表面張力計算方法を提供する。【選択図】  図２
2006030060.tif 000002
【特許請求の範囲】
【請求項１】
２次元以上の空間内の３種類の物質の界面張力を近似的に計算する計算方法において、考慮している領域全体で定義されたゼロから１の間の値を有する関数FとGと、該関数ＦとＧとの勾配の絶対値｜∇F｜及び|∇G|と、３つの実定数σ１、σ２、σ３に対して、エネルギー積分を微小体積要素dvolを利用して、
E=σ１∫ |∇F|G dvol+σ２∫|∇G|F dvol+σ３∫|∇G|(1-F)dvol    （数式１）
と定義し、該エネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式により３つの物質の表面張力を近似的に表現することを特徴とする表面張力計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
流体力学によって記述される現象を計算機上で模倣する計算技術に関する。また、自然現象を可視化するアニメーション技術等に関するものである。
【背景技術】
【０００２】
流体力学によって記述される自然現象を計算機上でより正確に模倣し、予測する事を目指す学問分野として計算流体力学と言われる分野がある。
【０００３】
本発明は、流体における表面張力の計算流体力学での取り扱いに関するものである（例えば、特許文献１参照）。
【０００４】
本発明では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"を元に従来の技術について述べる。
【０００５】
これらは、広い意味での流体体積法（VOF法と以下呼ぶ）という方法を利用したものであると考えてよい。領域全体で定義された状態の占有率関数F(X)を考える。占有率関数F(X)は物質１と物質２のみが存在する領域において、物質１の占有率を表現したものである。この関数に対して、エネルギー汎関数として、
E12=σ12∫|∇F|dvol      （数式２）
を考えることで、表面張力を表現できることを上記著者たちは発見した。σ12は、物質１と物質２の界面に働く表面エネルギーの係数、表面張力係数である。空間をデカルト座標（ｘ１，ｘ２，ｘ３）で表すこととする。（数式２）の汎関数に対して、ネーターの定理を適応して得られるi方向の運動量密度は、
fi =  ∇i (|∇F|)-Σｊ  ∇j ((∇j F)(∇i F)/|∇F|)      （数式３）
となる。ネーターの定理に関しては例えば、「量子場を学ぶための場の解析力学入門」高橋康著  講談社サイエンティフック  １９８２年発行のp.101(7.10)及び(7.3)に書かれている。スカラー場の場合に相当することを注意すると、上記の式（数式３）が得られる。これが所謂表面張力である。E12は物質１と物質２の界面の表面エネルギーを近似的に表現していることは重要なことであり、表面エネルギーは界面の面積に比例することは重要なことである。
【０００６】
即ち、表面張力は表面積を最小にする最小原理によって作られた力である。
【０００７】
他方、三種類の物質が混在する系において上記の考えを拡張することが上記文献に述べられている。三種類の物質が存在する系においては、図９に示すような接触角により、３つの界面が交わる三重点が特徴付けられることが知られている。（例えば、G.K.バチェラー著、橋本英典，松信八十男他訳の「入門流体力学」東京電気大学出版局  １９７４年発行に書かれている。）
この接触角を表現するのが本発明の目的である。本目的に対しても上記文献（J.U.Backbill, D.B.Kothe, C.ZemachらのJournal of Computational Physics 100巻  335-354 1992, p.341 D節）は第三の物質（以下物質３）が固体である場合に関して、計算方法を与えている。それによると、物質１と物質２の界面が、物質３に交わる接触角を幾何学的に求めて、それを所望の角度に変更することにより３重点の構成を行っている。
【特許文献１】特開平６-３０９２９４
【発明の開示】
【発明が解決しようとする課題】
【０００８】
本発明の課題は、計算流体力学において、接触角を如何に表現するかにある。従来技術で述べた方法は、接触角を幾何学的に捕らえるものである。しかしながら、上記文献（J.U.Backbill, D.B.Kothe, C.ZemachらのJournal of Computational Physics 100巻  335-354 1992, p.341 D節）に書れているように、ダイナミックスが絡んできた場合等、方法には任意性を伴うものである。平衡点においては、接触角は表面張力係数のみで決定されるが、ダイナミックスを追いかけている際には、接触角に物理的な意味はない。
【０００９】
また、図１０に示すように、重力中の３物質の状態において、下方向に重力が掛かっている際の接触角２１０１と２１０２は本来同じ接触角を持つことはないと考えるのが物理的に自然である。しかし、従来の方法ではこの差を取り込むことができない。
従って、より物理的な接触角を表現する方法が必要である。
【課題を解決するための手段】
【００１０】
本発明によると、２次元以上の空間内の３種類の物質の界面張力を近似的に計算する計算方法において、考慮している領域全体で定義されたゼロから１の間の値を有する関数FとGと、該関数ＦとＧとの勾配の絶対値｜∇F｜及び|∇G|と、３つの実定数σ１、σ２、σ３に対して、エネルギー積分を微小体積要素dvolを利用して、
E=σ１∫|∇F|G dvol+σ２∫|∇G|F dvol+σ３∫|∇G|(1-F)dvol      （数式４）
と定義し、該エネルギー積分Ｅを空間方向に対する汎関数微分により得られる微分方程式を差分化した式により３つの物質の表面張力を近似的に表現することを特徴とする表面張力計算アルゴリズムが提供される。
【発明の効果】
【００１１】
本発明によって、三相界面の問題をエネルギー的な視点から捉えることができ、非定常な場合に関しても、表面張力がエネルギー論的な意味から正確に計算できるようになった。例えば、図１０における重力下での接触角の同一性の矛盾に関しても、表面張力係数と重力係数との比として、２１０１部と２１０２部に相当するする部分の接触角のずれが生じるようにできる。これにより、エネルギー的観点から、問題点を克服できることとなった。
【発明を実施するための最良の形態】
【００１２】
本実施例は図１に示したフローチャートに従って本発明を計算機上で行う場合である。ここでいう計算機とは図７及び図８に示すようにCPU（１１６）やメモリ（１１８）等からなる演算処理部、ハードディスク（１１７）やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部そしてディスプレイ１１２等からなる出力表示部などで構成されているシステムである。各演算部等は、本体１１１に収まっている。
【００１３】
計算機はメモリ上に図１に示したフローチャートに沿って設計、コード化されたプログラムをロードし、必要な計算領域を確保し、適当な方法で入力された画像やパラメータに基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク等に書き込んで記憶させたり、ディスプレイ上に表示させる。
【００１４】
本発明の表面張力計算アルゴリズムは、次のエネルギー汎関数を基礎としている。
【００１５】
E=σ12∫|∇F|Gｄvol+σ13∫|∇G|Fｄvol+σ23∫|∇G|(1-F)ｄvol      （数式５）
但し、σ12、σ13、σ23は、それぞれ、表面張力係数であり、正の実数であり、dvolはルベーグ測度から決まる微小体積要素である。Ｆ及びGは[0,1]の値を持つ関数で、全領域で定義されている。また、ＦＧは物質１の占有率を表現しているとし、ＦＧが１の時、その場所は物質1が存在しているとする。（１−Ｆ）Ｇは物質２の占有率を表現し、（１−Ｆ）Ｇが１の時、、物質２がその点で占有しているとする。Ｇは物質３の占有率を表現しており、Ｇが０のとき、物質３が占有しているとする。また、それぞれが、０と１の中間の数値を取る場合は、界面を近似的に表現しているとする。よって、|∇F|及び|∇G|の台（非零値を取る定義域）は空間的に局在化している。
【００１６】
上記のエネルギー汎関数（数式５）の第一項の積分はＧが非零の領域（物質３の存在しない領域）で、物質１と物質２の界面の表面積を近似的に表現しており、それが小さい方がエネルギーが低いことを意味している。同様に、第二項の積分は、Ｆが非零の領域（物質１が存在する領域）で、物質１と物質３の界面の表面積を近似的に表現しており、それが小さい方がエネルギーが低いことを意味している。同様に第三項の積分も理解される。これらの項の影響はその係数σ12、σ13、σ23の大きさで決定され、それらが表面張力係数となっている。
【００１７】
ここで注意して置きたいこととして、表面張力の名前についてである。表面張力を何らかの引っ張りによる力として捕らえると理解が歪んでしまう。例えば、古典幾何光学において、フェルマーの最小原理により、真空中の２点間を結ぶ光線は直線をとなるわけであるが、現在、これを我々は張力が光線（光の軌道）にあるいは光に働いたとは理解せず、光路長を最小にするものと理解している。
【００１８】
しかし、歴史的には１６−１７世紀光路に張力が働いているという認識を持ち、弾性曲線と光の問題についての考察がなされていたようである。これを単純な最小原理に置き換えたのは、クラインによるとハミルトンの功績によるものである。（クライン：１９世紀の数学、Felix Klein著  弥永晶吉監修  足立恒雄・浪川幸彦・監訳  石井省吾・渡辺弘訳  共立出版1995年発行。）
物理的な理由として、なぜフェルマーの原理において、最小光路長という概念が張力という概念より自然であるかという理由は量子力学の立場から理解することによって明確になる。（ファインマン経路積分L.S.シュルマン著  高塚和夫訳  講談社95年発行に記述されている。）
同様に、ニュートン力学は、オイラー・ラグランジの最小原理によって、その軌道が決まるわけであるが、それも粒子の運動が何らかの張力が働くから、ニュートン力学に従うとは考えない。
【００１９】
もちろん、歴史的にはオイラーは、時間方向の点粒子の動きを時間軸を空間内の弾性棒に置き換えて、張力の影響として理解していたようであるが、我々は慣性力という場合はあるものの、張力によるものとは理解しない。
【００２０】
そのような対比の下で、表面張力の「張力」という言葉は歴史的なものとして理解する必要がある。表面張力はある表面エネルギーを最小化する際に生じる力であり、弾性体のようなフックの法則に従う意味での張力ではない。
【００２１】
本発明の思想としては、最小原理という現代的な理解に従うことで、計算流体力学の課題を解決しようとするものである。
【００２２】
本実施例の流れ図、図1について述べる。本実施例においては、流体のアルゴリズムに則り、本発明の表面張力計算アルゴリズムを取り入れた計算を行う。
【００２３】
まず、プログラムを開始すると、ステップS１で初期化を行う。
【００２４】
次にステップS2で、与えられた占有率関数FとGに対して、以下に述べるように表面張力を計算する。
【００２５】
また、ステップS3で適当な流体方程式を解き、FとGの物理状態を変更させる。
【００２６】
ステップS4で終了のチェックを行い、終了する場合はそのまま終える。終了しない場合は、ステップS1に戻る。これらを有限回行うことで、本発明の表面張力計算アルゴリズムを利用した流体の状態を模擬する。
【００２７】
＜実施形態１＞
本実施例では２次元の場合の本発明の表面計算アルゴリズムについて述べる。まず簡単のために、（数式５）のエネルギー汎関数の第一項のみを考える。第一項は
E12=σ12∫|∇F|G dvol=:σ12∫e12 dvol      （数式６）
と書かれる。σ12は、物質１と物質２の界面に働く表面エネルギーの係数、表面張力係数である。この項について本発明の表面張力に関する発明をのべれば、後の項に関しては同様に扱える。この汎関数に対して、ネーターの定理を適応して得られるi方向の運動量密度は、
fi =  ∇i (|∇F|G)-Σｊ  ∇j ((∇j F)(∇i F)G/|∇F|)      （数式７）
より、
fi =  Σｊ  ∇i ((∇j F)(∇j F)G/|∇F|)-Σｊ  ∇j ((∇j F)(∇i F)G/|∇F|)    （数式８）
となる。よって、流体の慣用に従って、
τii=  Σｊ≠i  ((∇j F)(∇j F)G/|∇F|)
τij=  -((∇i F)(∇j F)G/|∇F|)              （数式９）
と定義する。
【００２８】
fi =  ∇iτii- Σｊ  ∇jτij                  （数式１０）
を計算すればよい。本実施例では、これをスタッガード格子を利用して表現した。
【００２９】
まず、スカラー関数Fは２次元内の正方格子の胞に定義されている。胞は図３に示すように整数(n1,n2)で指定されている。それらの集合を[F]と書き、スタッガード格子内のスカラー場あるいは単にスカラー場と呼ぶ。また、胞の周りの辺x1方向のL1とx2方向のL2によって構成され、それぞれの線分L1,L2も、整数(n1,n2)によって指定される。Ｌ１とＬ２は図４及び図５に太線で示している。各L1にはx2方向のベクトル場V2が、L2にはｘ１方向のベクトル場V2が図６のように離散化されて定義されている。それらの集合を[V1,V2]と記し、スタッガード格子内のベクトル場あるいは単にベクトル場と呼ぶ。
【００３０】
このとき、離散化されたアトレッド(atled)演算子∇は、スカラー場に対してはベクトル場の空間への写像
∇：  [F]→  [V1,V2]                    （数式１１）
として、また、その像は、
V1(n1,n2):=F(n1+1,n2)-F(n1,n2)
V2(n1,n2):=F(n1,n2+1)-F(n1,n2)      （数式１２）
となるように定義される。この写像の像[V1，V2]を[∇F]と記す。また[V1，V2]から[V1]及び[V2]への自然な射影を利用して、[∇1F]及び[∇2F]をそれぞれの成分の集合として定義する。同様に、ベクトル場[V1，V2]に対するアトレッド(atled)演算子∇を、スカラー場の空間への写像、
∇：  [V1，V2]→  [F]      （数式１３）
として、またその像を
F(n1,n2):=V1(n1,n2)-V1(n1-1,n2)+V2(n1,n2)-V2(n1,n2-2)      （数式１４）
と定義する。
【００３１】
また、該ベクトル場[V1,V2]に対する絶対値｜｜は
｜｜：    [V1，V2]→  [F]      （数式１５）
となる写像として定義される。但し、
F(n1,n2)=([V1(n1+1,n2)+V1(n1,n2)]2/4+[V2(n1,n2+1)+V2(n1,n2)]2/4)1/2    （数式１６）
とする。
二つのスカラー場[F][G]の四則演算は、例えば掛け算が
[F][G]（n1,n2）=[F](n1,n2)[G]（n1,n2）      （数式１７）
として各定義点での値の掛け算とするように各点の値の四則演算として定義する。
これによりまず、エネルギー積分が胞上の和として、一意に
E=  Σn1,n2  (|[∇F]|[G])(n1,n2)      （数式１８）
のように定義される。これをスタッガード格子上のエネルギーと呼ぶ。
【００３２】
同時に、ｘ成分摂動を加えた際のネーターの定理を適応した式（数式１７）の第一項とである連続量としての∇i(|∇F|G)は、離散量∇i(|∇[F]|[G])としてスタガード格子上で一意に定義される。つまり、これをネーターの定理を適応した式（数式１７）の第一項のスタッガード格子上での変形物として考える。
【００３３】
ｘ成分摂動を加えた際のネーターの定理を適応した式（数式１７）の第二項に関しては、ネーターの定理を適用した際の項
-Σｊ  ∇j ((∇i F)δE/δ(∇j F))      （数式１９）
から得られる項である。ここで、この項は正確には
【００３４】
【数１】
2006030060.tif 000003
そこで、テンソル積の各部分について説明する。まず
δE/δ(∇j F)      （数式２０）
を考える。簡単のためにL2上の辺の場合のみ考察する。（n1,n2）で指定されるL2上では、
（F(n1+2,n2)-F(n1,n2)）G(n1+1,n2)/|∇F|(n1+1,n2)
+（F(n1+1,n2)-F(n1-1,n2)）G(n1,n2)/|∇F|(n1,n2)      （数式２１）
となる。これと同じ定義点を持つ∇１Fとのテンソル積は通常の意味で定義する。他方、これと直行する[∇２F]に関しては、定義点が異なるので、近傍の平均を取ったものをL2上のベクトルとしてテンソルする。つまり、L2上で定義されるV2として近傍の平均を
（数式２２）
V2：＝（V2（n1,n2）+V2(n1+1,n2)+V2(n1,n2-1)+V2(n+1,n2-1)）/4
と定義する。これにより、（F(n1+1,n2)-F(n1-1,n2)）G(n1,n2)/|∇F|(n1,n2)は[∇２F]とテンソル積
【００３５】
【数２】
2006030060.tif 000004
を取れる。これらを通常の意味で（実ベクトル空間として）のテンソル積
【００３６】
【数３】
2006030060.tif 000005
とする。
【００３７】
これにより、[V1,V2]に対するアトレッド(atled)演算子∇を適用でき、上記の演算は定義できる。但し、計算した値は胞上で定義されるので、胞から各辺L1、L2に近傍の平均化操作で定義点を変えることができ、計算が可能となる。
【００３８】
このようにして、（数式６）に対する表面張力（数式７）及び（数式１０）がスタッガード格子上で離散化されて計算できる。（数式５）の他の二項についても同様に計算でき、（数式５）に対するネーターの定理を利用した運動量場としての表面張力がスタガード格子上で計算でき、離散空間上の計算流体力学の計算スキームの中に取り入れられる。
【００３９】
図２は本発明のアルゴリズムを利用して、計算した表面エネルギーを考慮した形状の図である。Ｆ及びＧの等値線が描かれている。
【００４０】
＜実施形態２＞
同様にして３次元の計算に応用できる。
【００４１】
＜実施形態３＞
本発明の汎関数を利用すれば、F,Gを有限の関数基底で表現し、該関数基底で張られるベクトル空間内の最小問題として、表面張力を捕らえることが可能となる。よって、有限要素法等々に適応できる。
【図面の簡単な説明】
【００４２】
【図１】本実施例で使用した流れを示した流れ図。
【図２】本発明のアルゴリズムを利用して計算した本実施例１の計算結果。
【図３】本実施例１の計算領域の胞についての説明図。
【図４】本実施例１の計算領域の辺Ｌ２についての説明図。
【図５】本実施例１の計算領域の辺Ｌ１についての説明図。
【図６】本実施例１の計算領域の辺Ｌ１、辺Ｌ２で定義されるベクトルについての説明図。
【図７】本発明の亀裂幅算出方法を搭載した装置の説明図。
【図８】本発明の亀裂幅算出方法を搭載した装置の説明図。
【図９】従来の技術を説明する説明図。
【図１０】従来の技術を説明する説明図。
【符号の説明】
【００４３】
１１１  計算機の本体
１１２  ディスプレイ
１１３  フロッピー（登録商標）ディスクドライブ
１１４  マウス
１１５  キーボード
１１０  フロッピー（登録商標）ディスク
１１６  CPU
１１７  ハードディスク
１１８  メモリ
２１０１  接触角
２１０２  接触角
【図１】
2006030060.tif 000006
【図２】
2006030060.tif 000007
【図３】
2006030060.tif 000008
【図４】
2006030060.tif 000009
【図５】
2006030060.tif 000010
【図６】
2006030060.tif 000011
【図７】
2006030060.tif 000012
【図８】
2006030060.tif 000013
【図９】
2006030060.tif 000014
【図１０】
2006030060.tif 000015