(57)【要約】【課題】  離散化された格子上の計算流体問題において、固体近傍で粘性応力が巧く計算されないことが知られており、それを改善する計算する方法を提供する。【解決手段】  速度差の関数として粘性流を導入し、定義した粘性流の流入量が粘性応力であるとする定式化を行う。【選択図】  図１
2006170655.tif 000002
【特許請求の範囲】
【請求項１】
離散化された格子上での計算流体問題における、離散化された該格子上の速度の隣接する格子間の差と、該速度の格子間の差に付随する補正項により、粘性流を定義し、該粘性流の着目する格子胞への流入量を粘性力とする粘性応力計算方法において、該粘性流を二次元の場合は該補正項がゼロとし、三次元の場合は該粘性流の該補正項が粘性流の方向とは異なる速度成分を持たない事を特徴とする粘性応力計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
流体力学によって記述される現象を計算機上で模倣する計算技術に関する。また、自然現象を可視化するアニメーション技術等に関するものである。
【背景技術】
【０００２】
流体力学によって記述される自然現象を計算機上でより正確に模倣し、予測する事を目指す学問分野として計算流体力学と言われる分野がある。
【０００３】
本発明は、流体における粘性応力の計算流体力学での取り扱いに関するものである。特に、固体近傍での取り扱いに関するものである。以下では、固体と液体との界面を流体にとって壁として捉え、流体領域側の固液界面近傍のことを壁近傍と呼ぶ。壁近傍での流体方程式の取り扱いについては、計算流体力学の分野における標準としてFlow Scienece社作成のシミュレータFlow-3D内で採用されたFAVOR(Fractional-Area-Volume Obstracle Representation)法が基本である。よって、従来例としてFAVOR法をHirt, C.W., "Volume-fractiontechniques: powerful tools for wind engineering" Journal of Wind Engineering and Industrial Aerodynamics, 46 &47 (1993) 327-338の中身に沿ってその概略を述べることとする。
【０００４】
まず、粘性応力の標準的な非圧縮性流体力学での取り扱いについて述べる。（例えば、「流体力学」  巽友正著  培風館  1982年発行にある。）  流体力学では連続体力学としての取り扱いの観点から、通常テンソル場を対称テンソル場及び反対称テンソル場とで表現することが慣例となっており、粘性応力も
【０００５】
【数１】
2006170655.tif 000003
として表現されている。ここで、uiは流体の速度場であり、ｄは現在対象としている空間の次元であり、２か３であるとする。ηは粘性係数とする。各添え字、i, j, kは次元に応じて、１からｄまでを走る。またδijは、クロネッカーのδ指標で、i=jの時δij  =1で、それ以外のッ場合はゼロを返すものである。右辺の第一項と第二項とで対称テンソル場を表現している。
【０００６】
また、第三項は粘性応力による速度場の変化が非圧縮性流体の非圧縮性という性質と矛盾しないようにするために、体積保存の観点から導入される項である。よって、この式はＦi=η Σｊ∇τijと通常、表現されたりする。τij は、粘性応力テンソルでτij＝ｊ[∇ｊui +∇iuｊ -δijΣk[2/d] ∇kuk ]として与えられる。
【０００７】
上記文献においても、（数式１）の連続的な表式を基礎として粘性応力の取り扱いを行っている。それについて述べるために以下でまず、上記文献の概要について簡単に述べることとする。
【０００８】
上記文献に用いられている方法は流体体積法（VOF法と以下呼ぶ）を基礎としている。領域全体で定義された状態の占有率関数F(X)を考える。占有率関数F(X)は物質１と物質２のみが存在する領域において、物質１の占有率を表現したものであるが、壁近傍においてもどれくらい流体が占有するかをもこの考え方に準じて表される。
【０００９】
この上記文献での壁の取り扱いについて、図１１−図１３に沿って、その概念を述べる。まず、図１１は対象としている系を表している。２０００１は上記で述べた固体部分であり、以下壁と呼ぶことする。対して２０００２は対象とする流体部分である。２０００１は動かないものとする。これを離散力学系として、構造格子で表現し、流体部分の動力学を模擬することを考える。これが上記文献で扱った計算流体力学である。
【００１０】
簡単のために２次元空間の場合について説明を行う。図１２に示すように、構造格子に切って、図１１の一部を表現しているとする。この時、壁（固体部分）がどれくらい、その格子を横切るかの指標として開口率(Ax, Ay)を導入する。図１２に示しているのはその開口率の数値の分布図である。Axはｙ方向の節に張り付いたｘ方向の流れに対する開口率を表現するもので、Ayはｘ方向の節に張り付いてｙ方向の流れに対する開口率を表現するものである。開口率とは、もしも対象する１つの胞を構成する１つの節に対して、その節が壁を横切らなければ１の数字を当て、もしもその節が壁を横切った場合には壁でない領域の長さを節の全長で割ることで得られる数のことである（３次元では節を面と置き換えれば全く同様に議論できる。）。
【００１１】
例えば、速度は示強的な数であるが、運動量は示量的な数である。そのために図１２に示すように、密度１の流体の運動量を各節で評価する際には、速度場(ux, uy)に対して、（Ax・ux, Ay・uy）という量を扱う必要がある。
【００１２】
従って、流れ量の差が、示量的なものの変化となり、示量的な物理量と示強的な物理量との間にある種の関係を仮定する時には、上記の示量的な流れ量の差を評価する必要がある。連続の力学系ではこれは発散項div として表現される。しかしながら、上記に述べたように離散力学系においては、それぞれの節に対する流入量を正確に評価する必要がある。従って、例えば、流れ(jx, jy)に対して（但し、jxはｙ方向の節で、jyはｘ方向の節で定義されている数値列であるが）、 連続的な（div j）は、
【００１３】
【数２】
2006170655.tif 000004
となる。但し、Diは差分作用素で、上記のベクトル場に対しては
【００１４】
【数３】
2006170655.tif 000005
とする差分作用素である。
【００１５】
このような背景において、前記文献では粘性応力を
【００１６】
【数４】
2006170655.tif 000006
として表現している。簡便のためにここでは格子の添え字(m,n)は省略した。また、τij は、離散化された粘性応力テンソルでτij＝[Ｄｊui +Ｄiuｊ -δijΣk[2/d] Ｄkuk ]として差分系で与えられる（例えば、特許文献１参照）。
【特許文献１】特開平６−１９９５４号公報
【発明の開示】
【発明が解決しようとする課題】
【００１７】
近年、数ｍｍから数μｍ以下の微小素子内においてある種の液体等を利用して、新たな機能の発現を目指した試みが成されてきている。体積は代表長さの３乗で、また表面積は２乗でその大きさが変化することより、一般的に系のサイズが小さくなると界面の効果が重要となり、それらの素子の開発あるいは研究を行うに当たっては、壁の周りの粘性等の効果は重要性を増してくると考えられる。
【００１８】
他方、流体力学は連続力学として弾性体力学との対応の観点から歴史的な表記がなされてきているが、近年これらを現代数学的な観点から見直す試みが成されている。この視点から考えると従来の連続形式で定式化された表式は必ずしも離散力学系として適切なものとなっていない場合がある。
【００１９】
本発明の課題は、離散力学系として無矛盾な粘性応力項を導出し、計算流体力学ソルバとして利用し、より適切な力学系として、微小領域の流体力学系を模擬することになる。
【課題を解決するための手段】
【００２０】
本発明によると離散化された格子上での計算流体問題における、離散化された該格子上の速度の隣接する格子間の差と、該速度の格子間の差に付随する補正項により、粘性流を定義し、該粘性流の着目する格子胞への流入量を粘性力とする粘性応力計算方法において、該粘性流を二次元の場合は該補正項がゼロとし、三次元の場合は該粘性流の該補正項が粘性流の方向とは異なる速度成分を持たない事を特徴とする粘性応力計算方法が提供される。
【発明の効果】
【００２１】
本発明によって、粘性応力がフーリエの法則に従い離散化できたため、計算が基本法則と無矛盾となり、計算の安定性を生むこととなる。
【発明を実施するための最良の形態】
【００２２】
＜実施形態１＞
本実施例は図１に示したフローチャートに従って本発明を計算機上で行う場合である。ここでいう計算機とは図９及び図１０に示すようにCPU（１１６）やメモリ（１１８）等からなる演算処理部、ハードディスク（１１７）やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部そしてディスプレイ１１２等からなる出力表示部などで構成されているシステムである。各演算部等は、本体１１１に収まっている。
【００２３】
計算機はメモリ上に図１に示したフローチャートに沿って設計、コード化されたプログラムをロードし、必要な計算領域を確保し、適当な方法で入力された画像やパラメータに基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク等に書き込んで記憶させたり、ディスプレイ上に表示させる。
【００２４】
本発明の粘性応力計算アルゴリズムは、ハミルトン力学系ではない一般的な力学系において適応されているフーリエの法則をベースとするものである。
【００２５】
従来の差分方程式は、一般には図６（a）に示すように、連続の力学系に対して考察を行い、微分方程式や積分方程式のような連続の方程式を得た後に、それを差分近似することで差分方程式を得るものであるが、本発明の基本的な原理は、図６（ｂ）に示すように、連続の力学系に対する議論を離散力学系として解釈し直し、得られた離散力学系から離散的な方程式を得ようとするものである。これら二つの考え方には大きな隔たりがある。
【００２６】
本発明が基礎とするフーリエの法則についてここで述べる。フーリエの法則は一般には熱伝導方程式を得るためにフーリエが発見したものとして知られているが、この考えは一般の非ハミルトン力学系に対する方程式の導出の基本原理として一般に認識されている。特に非平衡統計力学系の関わる力学系において原理として十分に機能している。
【００２７】
ここでフーリエの法則について、熱の伝導現象を元に説明する。図７のように温度が離散的に空間上に分布している状況を考える。図７の着目する胞４００１において、フーリエの法則は、次のように定期式化される。隣との示強変数である温度の差比例して、単位時間に流れ込む熱の流れ量が矢印４００２や４００３のように与えられる。この比例定数が熱伝導率に関係する。矢印４００２と矢印４００３の差が着目する胞４００１に単位時間内に流れ込む熱の量ｄQ/dtである。つまり、各胞に左より番号をつけ、節にi+1/2という半整数を与えると
【００２８】
【数５】
2006170655.tif 000007
となる。熱量差ΔQと温度差ΔTとは比熱cを当して、ΔQ=cΔTという関係があるので、これを通して、温度の変化が起き、温度の変化あるいは熱の変化として熱伝導現象が説明される。
【００２９】
本発明においては、粘性応力について同様の考察を行うことから始める。粘性応力とは示強変数である速度差が空間的に生じた際に、それを緩和するために起こる現象であり、上記のフーリエの法則が適応できるものである。
【００３０】
本実施例では、これをスタッガード格子を利用して表現した。スタッガード格子とは圧力等のスカラー関数は２次元内の正方格子の胞に定義されている。胞は図２に示すように整数(n1,n2)で指定されている。また、胞の周りの辺x1方向のL1とx2方向のL2によって構成され、それぞれの線分L1,L2も、整数(n1,n2)によって指定される。Ｌ１とＬ２は図３及び図４に太線で示している。各L1にはx2方向のベクトル場V2が、L2にはｘ１方向のベクトル場V2が図５のように離散化されて定義されている。それらの集合を[V1,V2]と記し、スタッガード格子内のベクトル場あるいは単にベクトル場と呼ぶ。以下で利用する速度場は、このようにして定義されている。
【００３１】
まず、図８に示すように黒矢印で速度場（９００２〜６）が定義されているとする。着目する節（９００１）において速度差が生じているとする。点線で書いた胞は双対格子による双対胞である。
【００３２】
この時、示強変数である速度差に応じて、流れを緩和するための示量的な流れ（９００７〜１０）が生じると考えられる。これを粘性流Kj[ui]とここで呼ぶ。図では白抜きの矢印で与えた。これにより、運動量の緩和が起き、運動量を経て速度の変化が起きるというのが粘性応力をフーリエの法則において述べたものである。
【００３３】
これを式で表すとuiに対して、
【００３４】
【数６】
2006170655.tif 000008
となることが判る。Kj[ui]の主要項は差分系においてもＤｊuiであるが、非圧縮条件よる補正を考慮して、余分な項が付くこととなる。この際、連続では（数式１）と一致することを考慮すると、まず（数式１）を
【００３５】
【数７】
2006170655.tif 000009
と変形しておく。この時、[∇ｊui -Σk[(d-2)/d] ∇iuj ]は粘性流Kj[ui]と一致することが判る。よって、これを差分化して
【００３６】
【数８】
2006170655.tif 000010
となる。これは従来の例の（数式３）と異なっていることが判る。Ａは従来の技術の説明時に述べた開口部であり、従来の技術で述べたように図１１、１２、１３のように与えられる。この時、Aを除いて式の変形を行うと（数式７）と連続極限で一致することは自明である。
【００３７】
（数式８）は次元ｄが２の時、
【００３８】
【数９】
2006170655.tif 000011
となり、単純なフーリエの法則を適応したものと一致する。これはｄ＝２の場合の（数式３）と明らかに異なることが判る。
【００３９】
次に、次元d=3の場合を眺める。従来の技術との差を見るために（数式４）を具体的にi=1の場合を書き下すと
【００４０】
【数１０】
2006170655.tif 000012
となるのに対して、本発明の（数式８）は
【００４１】
【数１１】
2006170655.tif 000013
となる。本発明の粘性応力は、各方向に対して対称的に表現されていることがわかる。
【００４２】
より具体的に速度uiの速度差から起因する粘性流Kj[ui]を書き下すと、
【００４３】
【数１２】
2006170655.tif 000014
となり、粘性流の方向ｊに対して、Kj[ui]はuiとuj以外の方向の速度場を含んでいない。これは従来の例の（数式１０）の第一項とは著しく異なる性質である。
【００４４】
このように、本発明の粘性流は方向に対する依存性が少なく、対称性が高いものである。実際、もしも、速度場uiがあるスカラー場φの差分として表現されている場合、即ちui=Ｄiφと書けているという状況を想定する。この時、粘性流はKj[ui]＝（2／３）ＤiＤjφと書けることが判る。よって、
【００４５】
【数１３】
2006170655.tif 000015
で与えられる差分方程式はφの開口Ａを考慮した拡散方程式の差分型と自然な一致を見る。つまり、上記の仮定の元では（数式１３）は、
【００４６】
【数１４】
2006170655.tif 000016
となるので、適当な境界条件の後で、
【００４７】
【数１５】
2006170655.tif 000017
となり、開口Ａを考慮した拡散方程式の差分型となる。これは粘性応力の特徴を強く反映している。このような対称性は従来の技術の（数式４）や（数式１０）にはない対称性である。
【００４８】
ここで本発明の重要な点は、示量変数的な意味の流れとは何かを本発明ではきっちりと定義していることにある。つまり、示強的な流れの垂直方向の面積の積分が差分化された際の開口率Aとの積と近似しているという幾何学的な意味付けを考える際にこの視点は重要である。このような視点から考えると（数式４）におけるクロネッカーδで表現された項は、スカラー量に座標系に依存するクロネッカーδを添加した項であり、歴史的な意味こそあれ、現代微分幾何的視点から眺めると不自然である。これに比較して本発明の（数式８）が自然であるのは、座標によらない現代微分幾何的に記法を利用しても、表現できることである。（座標に依らない現代微分幾何的な表現については、例えば「理論物理学のための幾何学とトポロジー〈１〉」（中原幹夫著訳・佐久間  一浩訳  ピアソン・エデュケーション社、2000年出版）に記述されている。）  両者の差異は連続極限の元では重要ではないが、離散系においてはどこに開口率Ａを挿入するかという意味で重要となるのである。
【００４９】
つまり、従来の例では、図６（a）に示したように連続の力学系より連続の方程式を導出した後に、導出を行ったが、本発明では図６(b)のように離散力学系の従い離散系の方程式を導出した。それにより、従来の例のものとは異なる方程式を得た。
【００５０】
ここで本実施例の流れ図、図１について述べる。本実施例においては、流体のアルゴリズムに則り、本発明の粘性応力計算アルゴリズムを取り入れた計算を行う。
【００５１】
まず、プログラムを開始すると、ステップS１で初期化を行う。
【００５２】
次にステップS2で、与えられた速度場に対して、上記に述べたように粘性応力を計算する。
【００５３】
また、ステップS3で適当な流体方程式を解き、速度の物理状態を変更させる。
【００５４】
ステップS4で終了のチェックを行い、終了する場合はそのまま終える。終了しない場合は、ステップS1に戻る。これらを有限回行うことで、本発明の粘性応力計算アルゴリズムを利用した流体の状態を模擬する。
【００５５】
＜実施形態２＞
粘性定数ηが場所によって変化する場合、（数式８）の代りに
【００５６】
【数１６】
2006170655.tif 000018
を利用すると、速度場に対して、粘性応力がフーリエの法則に従って計算可能となる。
【００５７】
＜実施形態３＞
実施例１ではスタッガード格子に依存したものを実施例として示したが、本発明の本質は、格子系の特徴に強くは依存しない。そのため、他の格子系に対しても同様に本発明の粘性応力を拡張することは容易である。
【図面の簡単な説明】
【００５８】
【図１】本実施例で使用した流れを示した流れ図。
【図２】本実施例１の計算領域の胞についての説明図。
【図３】本実施例１の計算領域の辺Ｌ２についての説明図。
【図４】本実施例１の計算領域の辺Ｌ１についての説明図。
【図５】本実施例１の計算領域の辺Ｌ１、辺Ｌ２で定義されるベクトルについての説明図。
【図６】従来の考え方と本発明の基本的な考え方を示した図。
【図７】本発明の基本となるフーリエの法則についての説明図。
【図８】本発明の説明図。
【図９】本発明の亀裂幅算出方法を搭載した装置の説明図。
【図１０】本発明の亀裂幅算出方法を搭載した装置の説明図。
【図１１】従来の技術を説明する説明図。
【図１２】従来の技術を説明する説明図。
【図１３】従来の技術を説明する説明図。
【符号の説明】
【００５９】
１１１  計算機の本体
１１２  ディスプレイ
１１３  フロッピー（登録商標）ディスクドライブ
１１４  マウス
１１５  キーボード
１１０  フロッピー（登録商標）ディスク
１１６  ＣＰＵ
１１７  ハードディスク
１１８  メモリ
２００１  流体部
２００２  固体部分
【図１】
2006170655.tif 000019
【図２】
2006170655.tif 000020
【図３】
2006170655.tif 000021
【図４】
2006170655.tif 000022
【図５】
2006170655.tif 000023
【図６】
2006170655.tif 000024
【図７】
2006170655.tif 000025
【図８】
2006170655.tif 000026
【図９】
2006170655.tif 000027
【図１０】
2006170655.tif 000028
【図１１】
2006170655.tif 000029
【図１２】
2006170655.tif 000030
【図１３】
2006170655.tif 000031