固液界面における壁粘性応力計算方法
【課題】 格子に対して斜行した固液界面であっても、また固液界面が移動する場合であっても、格子に依存しない、高精度な壁粘性応力の計算方法を提供する。
【解決手段】 固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとをもって固液界面における壁粘性応力を計算する。
【解決手段】 固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとをもって固液界面における壁粘性応力を計算する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、離散格子点で構成された固体・液体の二相が共存する流れを数値的に模擬する方法に関し、特に固液界面における壁粘性応力を近似的に計算する方法に関する。
【背景技術】
【０００２】
異なる相が混在する流れを数値的に解く手法のひとつとして、占有率関数を用いる手法が知られている。占有率関数とは、各格子点における気体、液体、固体のそれぞれの占有率をあらわしたものである。この占有率関数を用いた手法は、ＶＯＦ法（C. W. Hirt and B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comp. Phys., 39, 201(1981)）の流れをくむものであり、液滴の分裂・合体といった問題や、界面が大変形するような問題に対しても適用可能であることから、現在に至っても種々の改良を施されながら用いられている。
【０００３】
占有率関数を用いた手法のもとでの壁粘性応力の計算方法の中で、特にレベルセット関数を併用して計算精度の向上をはかる手法が知られている。代表的なものはＣＬＳＶＯＦ法（M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comp. Phys. 187(2003), 110-136）である。図９はＣＬＳＶＯＦ法における壁粘性応力の計算方法を説明する図である。９０１は計算格子の胞、９０２は壁、９０３は固液界面、９０４は離散速度および壁粘性応力のｘ方向成分が定義される辺、９０５は前記定義点（９０４）から固液界面（９０３）までの距離ｈ、９０６は前記胞（９０１）内部の固液界面の面積Ａｓｕｒｆをそれぞれ表す。ＣＬＳＶＯＦ法における壁粘性横領の計算方法では、図９にて定義された量を用いて、ｘ、ｙ方向の壁粘性応力ｗｓｘ、ｗｓｙをそれぞれ
【０００４】
【数１】
【０００５】
と定義する。ここでηは固液界面における粘性応力係数である。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開2006-170655号公報
【非特許文献】
【０００７】
【非特許文献１】C. W. Hirt and B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comp. Phys., 39, 201(1981)
【非特許文献２】M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comp. Phys. 187(2003), 110-136
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上記従来例における手法においては、主要な二つの問題点が知られている。第一に、格子に対して斜行した固液界面に対する壁粘性応力の計算の問題がある。ＣＬＳＶＯＦ法における壁粘性応力の計算方法は、固液界面近傍の離散流速が固液界面の法線ベクトルと直交していることを仮定しているが、一般に離散的な系ではこの仮定が妥当でない場合が多く、計算精度が悪化する場合がある。第二に、壁が移動する場合の壁粘性応力の計算方法の問題がある。上述したように、ＣＬＳＶＯＦ法における壁粘性応力の計算においては固液界面近傍の離散流速が固液界面の法線ベクトルと直交していることを仮定しているが、この仮定により壁の移動に起因する流速成分をうまく取り込むことができない。
【０００９】
本発明は以上のような課題に鑑みてなされたもので、その目的は、格子に対して斜行した固液界面であっても、また固液界面が移動する場合であっても、格子に依存しない、高精度な壁粘性応力の計算方法を提供することにある。
【課題を解決するための手段】
【００１０】
上記課題を解決するため、本発明では、離散格子点で構成された領域上で固体・液体の二相が共存する流れを数値的に模擬する方法において、前記固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとをもって固液界面における壁粘性応力を計算するとした。
【発明の効果】
【００１１】
以上に説明したように、本発明により、格子に対して斜行した固液界面であっても、また固液界面が移動する場合であっても、格子に依存せず、安定に精度良く壁粘性応力を計算できるようになる。
【図面の簡単な説明】
【００１２】
【図１】本実施例にて使用した流体計算手法の流れ図
【図２】本実施例にて使用した計算機のブロック図
【図３】本実施例にて使用したスタガード格子の模式図
【図４】格子上において壁粘性応力が定義される点を示した図
【図５】胞中心に補間された流速u*を示す図
【図６】壁面の法線ベクトルと直交するように射影された流速を示す図
【図７】厚さτを有する固液界面領域Ｖと射影された流速のｘ成分から仮想的な胞内の壁粘性応力の平均値のｘ成分を求めるステップの模式図
【図８】辺の開口率Axとｘ方向の壁粘性応力の平均値からｘ方向の壁粘性応力ＷＳi+1/2,jを求めるステップの模式図
【図９】ＣＬＳＶＯＦ法における壁粘性応力の計算方法を説明する図
【発明を実施するための形態】
【００１３】
以下、本発明の原理について説明する。一般に壁粘性応力は、固液界面の十分近傍において、固液界面の法線ベクトルと直交する流速成分（すなわち固液界面と平行な流速成分）に比例するような力として計算される。このとき、前記成分以外の他の流速成分からの壁粘性応力への寄与は０であるとされる。一方離散格子点で構成された領域上では、固液界面から最も近い定義点における離散流速であっても、一般には固液界面に平行な流速成分と他の方向の流速成分とが混じり合っている。ゆえに、離散流速を固液界面の法線ベクトルに直交するように射影し、所望の成分のみを抽出するステップを経ることが必要である。なおここで、固液界面が任意の方向に任意の速度で移動する場合には、移動する固液界面から見た相対的な離散流速を射影すればよい。また壁粘性応力は、先に述べたように、固液界面の十分近傍の領域のみでデルタ関数的に働く力である。しかし本発明は離散格子点で構成された領域上でのものであるため、デルタ関数を厳密には扱えない。よって、固液界面に任意の厚さを持たせた固液界面領域を定義し、この領域内に壁粘性応力が働くとするのが自然である。上述のような着目点に留意して考察を進めた結果、本発明の発明者らは、以下に示す数式２・数式３を用いて、任意の壁面形状近傍での壁粘性応力を格子によらず精度良く計算できることを見いだした。
【００１４】
【数２】
【００１５】
【数３】
【００１６】
ここでηは粘性係数、Ｖは固液界面領域の体積、ｈは壁面から流速定義点までの距離、
【００１７】
【数４】
【００１８】
は射影された流速、
【００１９】
【数５】
【００２０】
は仮想的な単位胞における壁粘性応力の平均値、Ａは開口率、Δｘ、Δｙは格子間隔、そしてＷＳが壁粘性応力を表す。
【実施例１】
【００２１】
図１に示したようなフローチャートに従い、本発明を計算機上にて行う第一の実施例について説明する。本実施例では、おおまかにわけて、開始（１０１）、初期化（１０２）、壁粘性応力を計算するブロック（１０３）、壁粘性応力を含む流体計算を行うブロック（１０４）、時間更新を行うブロック（１０５）、終了判定を行うブロック（１０６）、からなる。終了判定が真ならば終了（１０７）であり、偽ならば初期化（１０２）へ戻る。なおここでいう計算機は、図２に示すような構成からなり、ＣＰＵ（２０１）、メモリ（２０２）、ハードディスクドライブ（２０３）、出力装置（２０４）、入力装置（２０５）らがデータバス（２０６）を介して結合されたものであって、たとえば周知のパーソナルコンピュータやワークステーション等に相当する。以後は図１のフローチャート中の、壁粘性応力の計算のステップについて説明する。
【００２２】
いま、系として２次元の系を考え、ここにスタガード格子が配されているとする。図３に格子の様子を示す。胞（３０１）には圧力等のスカラー量が定義される。この胞の位置をｘ方向にｍ番目、ｙ方向にｎ番目というラベル付けを行う。よってこの胞の位置は（ｍ，ｎ）のように表される。同様に辺（３０２）には速度や運動量などのベクトル量がそれぞれ定義される。この辺の位置は（ｍ＋１／２，ｎ），（ｍ，ｎ＋１／２）のように表される。また節（３０３）には、必要に応じて種々の物理量が補間されて配置される。この節の位置は（ｍ＋１／２，ｎ＋１／２）と表される。格子間隔はｘ方向、ｙ方向共に一様であるとし、それぞれΔｘ（３０４）、Δｙ（３０５）と表す。また、本発明にて計算される壁粘性応力の定義点を図４に示す。ｘ方向の壁粘性応力をＷＳi+1/2,j（４０２）、ｙ方向の壁粘性応力をＷＳi,j+1/2（４０４）とし、それぞれの定義点を含む仮想的な胞（４０１）、（４０３）を定義する。
【００２３】
以後、ｘ方向の壁粘性応力を計算する方法について詳細に説明する。ｙ方向については、以下のｘ方向についての説明と全く同様にして計算することができるため、説明は省略する。上述したように、本実施例で用いるスタガード格子では、流速が辺に定義されている。これをもとに、胞中心における流速を補間により求める。この様子を図５に示す。ここで５０１は壁、５０２は胞中心における流速u*i,j、５０３は胞中心から壁面までの距離ｈである。胞中心における流速u*i,jは、数式４に示すように線形補間によって求められる。
【００２４】
【数６】
【００２５】
続いて、求めた胞中心における流速u*i,jを、壁面の法線ベクトルと直交するように射影する。この様子を図６に示す。６０１は壁面の法線ベクトルと直交するように射影された胞中心における流速
【００２６】
【数７】
【００２７】
である。なおここで、壁面が任意の速度v=(vx, vy)で移動している場合には、u*i,jを射影するのではなく、壁面に対する相対的な流速u*i,j−vを射影すればよい。
【００２８】
このようにして壁面の法線ベクトルと直交するように射影された流速
【００２９】
【数８】
【００３０】
を用い、発明の実施の形態の項で説明した数式２によって壁粘性応力の平均値が計算される。この様子を図７に示す。７０１は射影された胞中心における流速
【００３１】
【数９】
【００３２】
のｘ方向成分、７０５は射影された胞中心における流速
【００３３】
【数１０】
【００３４】
のｘ方向成分、７０２は着目する固液界面領域であってその体積をＶとし、７０３は仮想的な胞内における平均の壁粘性応力
【００３５】
【数１１】
【００３６】
、７０４は固液界面領域の厚さτである。本実施例では、固液界面領域の厚さτは格子間隔よりも小さいとした。上記にて定義された量を用い、数式２によって仮想的な胞内における平均の壁粘性応力
【００３７】
【数１２】
【００３８】
を計算する。すなわち、
【００３９】
【数１３】
【００４０】
である。なおこの数式２による平均の壁粘性応力の計算は、固液界面領域Ｖが存在するすべての仮想的な胞において行われる。
【００４１】
こうして求めた平均の壁粘性応力
【００４２】
【数１４】
【００４３】
を、辺上での壁粘性応力ＷＳi+1/2,jに変換する。この様子を図８に示す。８０１は開口率Ａｘである。この開口率と壁粘性応力の平均値から、数式３に従って、壁粘性応力のｘ成分ＷＳi+1/2,jが計算される。すなわち、
【００４４】
【数１５】
【００４５】
である。
【００４６】
以上のようにして、壁粘性応力を二次元スタガード格子上で計算でき、図１で示したような計算スキームに組み込むことが可能になった。
【００４７】
また、本実施例は二次元スタガード格子上のものであったので、添え字の組（ｍ，ｎ）＝｛１，２｝であったが、これを三次元化するには添え字の組を（ｍ，ｎ）＝｛１，２，３｝とするだけでよい。差分化についても、上記実施例から容易に類推することが可能であるため、三次元の計算への適用も容易に行うことができる。
【実施例２】
【００４８】
次に、本発明の第二の実施例について述べる。なお、以後においては第一の実施例との相違点にのみ注目して説明し、その他第一の実施例と同様の部分については説明を省略する。第二の実施例では、固液界面領域の厚さτを、格子間隔の数倍〜十数倍程度と大きく取る。ＶＯＦ法に代表される、物体を占有率関数を用いて表現する計算手法においては、計算の安定性を確保するために占有率関数の遷移領域幅を格子間隔より広く取り、数格子点分の遷移領域を取ることがある。このような場合には、この占有率関数の遷移領域幅と固液界面領域の厚さとを同じ程度の長さに取ることで、計算の安定性及び精度の向上を図ることができる。なおこの際は、開口率Ａｘ，Ａｙは経験的に決めてよい。
【００４９】
こうして、固液界面領域の厚さを大きくする他は、第一の実施例とまったく同様にして二次元スタガード格子上での差分化を行い、図１で示したような計算スキームに組み込むことが可能になる。以上のようにして、実際に第一の実施例と同様の計算結果を得ることができる。
【００５０】
なお、本発明の思想に添ったものであれば、上記実施例に限定されるものではなく、たとえば物理量の定義点によらない。更には、差分法のみならず有限体積法や有限要素法に対しても本発明を適用することが可能である。
【符号の説明】
【００５１】
１０１ 計算開始時の処理を行う部分
１０２ 系の初期化を行う部分
１０３ 壁粘性応力を計算するブロック
１０４ 壁粘性応力を含む流体計算を行うブロック
１０５ 時間を更新するブロック
１０６ 終了判定を行う部分
１０７ 計算終了時の処理を行う部分
２０１ ＣＰＵ
２０２ メモリ
２０３ ハードディスクドライブ（ＨＤＤ）
２０４ 出力装置
２０５ 入力装置
２０６ データバス
３０１ 胞
３０２ 辺
３０３ 節
３０４ ｘ方向格子間隔Δｘ
３０５ ｙ方向格子間隔Δｙ
４０１ ｘ方向に働く壁粘性応力ＷＳi+1/2,jが定義される点を中心とする仮想的な胞
４０２ ｘ方向に働く壁粘性応力ＷＳi+1/2,jが定義される点
４０３ ｙ方向に働く壁粘性応力ＷＳi,j+1/2が定義される点を中心とする仮想的な胞
４０４ ｙ方向に働く壁粘性応力ＷＳi,j+1/2が定義される点
５０１ 壁
５０２ 胞中心(i,j)に補間された流速u*i,j
５０３ 壁面から胞中心までの距離h
６０１ u*i,jを壁面の法線ベクトルと直交するように射影した流速
７０１ ｘ成分
７０２ 固液界面領域の体積Ｖ
７０３ 仮想的な単位胞４０１内における壁粘性応力のｘ方向の平均値
７０４ 固液界面領域の厚みτ
７０５ ｘ成分
８０１ 開口率Ax
９０１ 胞
９０２ 壁
９０３ 固液界面
９０４ 離散流速のｘ方向成分および壁粘性応力のｘ方向成分の定義点
９０５ 離散流速のｘ方向成分および壁粘性応力のｘ方向成分の定義点から固液界面までの距離
９０６ 胞内の固液界面の面積Ａｓｕｒｆ
【技術分野】
【０００１】
本発明は、離散格子点で構成された固体・液体の二相が共存する流れを数値的に模擬する方法に関し、特に固液界面における壁粘性応力を近似的に計算する方法に関する。
【背景技術】
【０００２】
異なる相が混在する流れを数値的に解く手法のひとつとして、占有率関数を用いる手法が知られている。占有率関数とは、各格子点における気体、液体、固体のそれぞれの占有率をあらわしたものである。この占有率関数を用いた手法は、ＶＯＦ法（C. W. Hirt and B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comp. Phys., 39, 201(1981)）の流れをくむものであり、液滴の分裂・合体といった問題や、界面が大変形するような問題に対しても適用可能であることから、現在に至っても種々の改良を施されながら用いられている。
【０００３】
占有率関数を用いた手法のもとでの壁粘性応力の計算方法の中で、特にレベルセット関数を併用して計算精度の向上をはかる手法が知られている。代表的なものはＣＬＳＶＯＦ法（M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comp. Phys. 187(2003), 110-136）である。図９はＣＬＳＶＯＦ法における壁粘性応力の計算方法を説明する図である。９０１は計算格子の胞、９０２は壁、９０３は固液界面、９０４は離散速度および壁粘性応力のｘ方向成分が定義される辺、９０５は前記定義点（９０４）から固液界面（９０３）までの距離ｈ、９０６は前記胞（９０１）内部の固液界面の面積Ａｓｕｒｆをそれぞれ表す。ＣＬＳＶＯＦ法における壁粘性横領の計算方法では、図９にて定義された量を用いて、ｘ、ｙ方向の壁粘性応力ｗｓｘ、ｗｓｙをそれぞれ
【０００４】
【数１】
【０００５】
と定義する。ここでηは固液界面における粘性応力係数である。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開2006-170655号公報
【非特許文献】
【０００７】
【非特許文献１】C. W. Hirt and B. D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comp. Phys., 39, 201(1981)
【非特許文献２】M. Sussman, A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles, J. Comp. Phys. 187(2003), 110-136
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上記従来例における手法においては、主要な二つの問題点が知られている。第一に、格子に対して斜行した固液界面に対する壁粘性応力の計算の問題がある。ＣＬＳＶＯＦ法における壁粘性応力の計算方法は、固液界面近傍の離散流速が固液界面の法線ベクトルと直交していることを仮定しているが、一般に離散的な系ではこの仮定が妥当でない場合が多く、計算精度が悪化する場合がある。第二に、壁が移動する場合の壁粘性応力の計算方法の問題がある。上述したように、ＣＬＳＶＯＦ法における壁粘性応力の計算においては固液界面近傍の離散流速が固液界面の法線ベクトルと直交していることを仮定しているが、この仮定により壁の移動に起因する流速成分をうまく取り込むことができない。
【０００９】
本発明は以上のような課題に鑑みてなされたもので、その目的は、格子に対して斜行した固液界面であっても、また固液界面が移動する場合であっても、格子に依存しない、高精度な壁粘性応力の計算方法を提供することにある。
【課題を解決するための手段】
【００１０】
上記課題を解決するため、本発明では、離散格子点で構成された領域上で固体・液体の二相が共存する流れを数値的に模擬する方法において、前記固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとをもって固液界面における壁粘性応力を計算するとした。
【発明の効果】
【００１１】
以上に説明したように、本発明により、格子に対して斜行した固液界面であっても、また固液界面が移動する場合であっても、格子に依存せず、安定に精度良く壁粘性応力を計算できるようになる。
【図面の簡単な説明】
【００１２】
【図１】本実施例にて使用した流体計算手法の流れ図
【図２】本実施例にて使用した計算機のブロック図
【図３】本実施例にて使用したスタガード格子の模式図
【図４】格子上において壁粘性応力が定義される点を示した図
【図５】胞中心に補間された流速u*を示す図
【図６】壁面の法線ベクトルと直交するように射影された流速を示す図
【図７】厚さτを有する固液界面領域Ｖと射影された流速のｘ成分から仮想的な胞内の壁粘性応力の平均値のｘ成分を求めるステップの模式図
【図８】辺の開口率Axとｘ方向の壁粘性応力の平均値からｘ方向の壁粘性応力ＷＳi+1/2,jを求めるステップの模式図
【図９】ＣＬＳＶＯＦ法における壁粘性応力の計算方法を説明する図
【発明を実施するための形態】
【００１３】
以下、本発明の原理について説明する。一般に壁粘性応力は、固液界面の十分近傍において、固液界面の法線ベクトルと直交する流速成分（すなわち固液界面と平行な流速成分）に比例するような力として計算される。このとき、前記成分以外の他の流速成分からの壁粘性応力への寄与は０であるとされる。一方離散格子点で構成された領域上では、固液界面から最も近い定義点における離散流速であっても、一般には固液界面に平行な流速成分と他の方向の流速成分とが混じり合っている。ゆえに、離散流速を固液界面の法線ベクトルに直交するように射影し、所望の成分のみを抽出するステップを経ることが必要である。なおここで、固液界面が任意の方向に任意の速度で移動する場合には、移動する固液界面から見た相対的な離散流速を射影すればよい。また壁粘性応力は、先に述べたように、固液界面の十分近傍の領域のみでデルタ関数的に働く力である。しかし本発明は離散格子点で構成された領域上でのものであるため、デルタ関数を厳密には扱えない。よって、固液界面に任意の厚さを持たせた固液界面領域を定義し、この領域内に壁粘性応力が働くとするのが自然である。上述のような着目点に留意して考察を進めた結果、本発明の発明者らは、以下に示す数式２・数式３を用いて、任意の壁面形状近傍での壁粘性応力を格子によらず精度良く計算できることを見いだした。
【００１４】
【数２】
【００１５】
【数３】
【００１６】
ここでηは粘性係数、Ｖは固液界面領域の体積、ｈは壁面から流速定義点までの距離、
【００１７】
【数４】
【００１８】
は射影された流速、
【００１９】
【数５】
【００２０】
は仮想的な単位胞における壁粘性応力の平均値、Ａは開口率、Δｘ、Δｙは格子間隔、そしてＷＳが壁粘性応力を表す。
【実施例１】
【００２１】
図１に示したようなフローチャートに従い、本発明を計算機上にて行う第一の実施例について説明する。本実施例では、おおまかにわけて、開始（１０１）、初期化（１０２）、壁粘性応力を計算するブロック（１０３）、壁粘性応力を含む流体計算を行うブロック（１０４）、時間更新を行うブロック（１０５）、終了判定を行うブロック（１０６）、からなる。終了判定が真ならば終了（１０７）であり、偽ならば初期化（１０２）へ戻る。なおここでいう計算機は、図２に示すような構成からなり、ＣＰＵ（２０１）、メモリ（２０２）、ハードディスクドライブ（２０３）、出力装置（２０４）、入力装置（２０５）らがデータバス（２０６）を介して結合されたものであって、たとえば周知のパーソナルコンピュータやワークステーション等に相当する。以後は図１のフローチャート中の、壁粘性応力の計算のステップについて説明する。
【００２２】
いま、系として２次元の系を考え、ここにスタガード格子が配されているとする。図３に格子の様子を示す。胞（３０１）には圧力等のスカラー量が定義される。この胞の位置をｘ方向にｍ番目、ｙ方向にｎ番目というラベル付けを行う。よってこの胞の位置は（ｍ，ｎ）のように表される。同様に辺（３０２）には速度や運動量などのベクトル量がそれぞれ定義される。この辺の位置は（ｍ＋１／２，ｎ），（ｍ，ｎ＋１／２）のように表される。また節（３０３）には、必要に応じて種々の物理量が補間されて配置される。この節の位置は（ｍ＋１／２，ｎ＋１／２）と表される。格子間隔はｘ方向、ｙ方向共に一様であるとし、それぞれΔｘ（３０４）、Δｙ（３０５）と表す。また、本発明にて計算される壁粘性応力の定義点を図４に示す。ｘ方向の壁粘性応力をＷＳi+1/2,j（４０２）、ｙ方向の壁粘性応力をＷＳi,j+1/2（４０４）とし、それぞれの定義点を含む仮想的な胞（４０１）、（４０３）を定義する。
【００２３】
以後、ｘ方向の壁粘性応力を計算する方法について詳細に説明する。ｙ方向については、以下のｘ方向についての説明と全く同様にして計算することができるため、説明は省略する。上述したように、本実施例で用いるスタガード格子では、流速が辺に定義されている。これをもとに、胞中心における流速を補間により求める。この様子を図５に示す。ここで５０１は壁、５０２は胞中心における流速u*i,j、５０３は胞中心から壁面までの距離ｈである。胞中心における流速u*i,jは、数式４に示すように線形補間によって求められる。
【００２４】
【数６】
【００２５】
続いて、求めた胞中心における流速u*i,jを、壁面の法線ベクトルと直交するように射影する。この様子を図６に示す。６０１は壁面の法線ベクトルと直交するように射影された胞中心における流速
【００２６】
【数７】
【００２７】
である。なおここで、壁面が任意の速度v=(vx, vy)で移動している場合には、u*i,jを射影するのではなく、壁面に対する相対的な流速u*i,j−vを射影すればよい。
【００２８】
このようにして壁面の法線ベクトルと直交するように射影された流速
【００２９】
【数８】
【００３０】
を用い、発明の実施の形態の項で説明した数式２によって壁粘性応力の平均値が計算される。この様子を図７に示す。７０１は射影された胞中心における流速
【００３１】
【数９】
【００３２】
のｘ方向成分、７０５は射影された胞中心における流速
【００３３】
【数１０】
【００３４】
のｘ方向成分、７０２は着目する固液界面領域であってその体積をＶとし、７０３は仮想的な胞内における平均の壁粘性応力
【００３５】
【数１１】
【００３６】
、７０４は固液界面領域の厚さτである。本実施例では、固液界面領域の厚さτは格子間隔よりも小さいとした。上記にて定義された量を用い、数式２によって仮想的な胞内における平均の壁粘性応力
【００３７】
【数１２】
【００３８】
を計算する。すなわち、
【００３９】
【数１３】
【００４０】
である。なおこの数式２による平均の壁粘性応力の計算は、固液界面領域Ｖが存在するすべての仮想的な胞において行われる。
【００４１】
こうして求めた平均の壁粘性応力
【００４２】
【数１４】
【００４３】
を、辺上での壁粘性応力ＷＳi+1/2,jに変換する。この様子を図８に示す。８０１は開口率Ａｘである。この開口率と壁粘性応力の平均値から、数式３に従って、壁粘性応力のｘ成分ＷＳi+1/2,jが計算される。すなわち、
【００４４】
【数１５】
【００４５】
である。
【００４６】
以上のようにして、壁粘性応力を二次元スタガード格子上で計算でき、図１で示したような計算スキームに組み込むことが可能になった。
【００４７】
また、本実施例は二次元スタガード格子上のものであったので、添え字の組（ｍ，ｎ）＝｛１，２｝であったが、これを三次元化するには添え字の組を（ｍ，ｎ）＝｛１，２，３｝とするだけでよい。差分化についても、上記実施例から容易に類推することが可能であるため、三次元の計算への適用も容易に行うことができる。
【実施例２】
【００４８】
次に、本発明の第二の実施例について述べる。なお、以後においては第一の実施例との相違点にのみ注目して説明し、その他第一の実施例と同様の部分については説明を省略する。第二の実施例では、固液界面領域の厚さτを、格子間隔の数倍〜十数倍程度と大きく取る。ＶＯＦ法に代表される、物体を占有率関数を用いて表現する計算手法においては、計算の安定性を確保するために占有率関数の遷移領域幅を格子間隔より広く取り、数格子点分の遷移領域を取ることがある。このような場合には、この占有率関数の遷移領域幅と固液界面領域の厚さとを同じ程度の長さに取ることで、計算の安定性及び精度の向上を図ることができる。なおこの際は、開口率Ａｘ，Ａｙは経験的に決めてよい。
【００４９】
こうして、固液界面領域の厚さを大きくする他は、第一の実施例とまったく同様にして二次元スタガード格子上での差分化を行い、図１で示したような計算スキームに組み込むことが可能になる。以上のようにして、実際に第一の実施例と同様の計算結果を得ることができる。
【００５０】
なお、本発明の思想に添ったものであれば、上記実施例に限定されるものではなく、たとえば物理量の定義点によらない。更には、差分法のみならず有限体積法や有限要素法に対しても本発明を適用することが可能である。
【符号の説明】
【００５１】
１０１ 計算開始時の処理を行う部分
１０２ 系の初期化を行う部分
１０３ 壁粘性応力を計算するブロック
１０４ 壁粘性応力を含む流体計算を行うブロック
１０５ 時間を更新するブロック
１０６ 終了判定を行う部分
１０７ 計算終了時の処理を行う部分
２０１ ＣＰＵ
２０２ メモリ
２０３ ハードディスクドライブ（ＨＤＤ）
２０４ 出力装置
２０５ 入力装置
２０６ データバス
３０１ 胞
３０２ 辺
３０３ 節
３０４ ｘ方向格子間隔Δｘ
３０５ ｙ方向格子間隔Δｙ
４０１ ｘ方向に働く壁粘性応力ＷＳi+1/2,jが定義される点を中心とする仮想的な胞
４０２ ｘ方向に働く壁粘性応力ＷＳi+1/2,jが定義される点
４０３ ｙ方向に働く壁粘性応力ＷＳi,j+1/2が定義される点を中心とする仮想的な胞
４０４ ｙ方向に働く壁粘性応力ＷＳi,j+1/2が定義される点
５０１ 壁
５０２ 胞中心(i,j)に補間された流速u*i,j
５０３ 壁面から胞中心までの距離h
６０１ u*i,jを壁面の法線ベクトルと直交するように射影した流速
７０１ ｘ成分
７０２ 固液界面領域の体積Ｖ
７０３ 仮想的な単位胞４０１内における壁粘性応力のｘ方向の平均値
７０４ 固液界面領域の厚みτ
７０５ ｘ成分
８０１ 開口率Ax
９０１ 胞
９０２ 壁
９０３ 固液界面
９０４ 離散流速のｘ方向成分および壁粘性応力のｘ方向成分の定義点
９０５ 離散流速のｘ方向成分および壁粘性応力のｘ方向成分の定義点から固液界面までの距離
９０６ 胞内の固液界面の面積Ａｓｕｒｆ
【特許請求の範囲】
【請求項１】
離散格子点で構成された領域上で固体・液体の二相が共存する流れを数値的に模擬する固液界面における壁粘性応力計算方法において、前記固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとを有することを特徴とする固液界面における壁粘性応力計算方法。
【請求項１】
離散格子点で構成された領域上で固体・液体の二相が共存する流れを数値的に模擬する固液界面における壁粘性応力計算方法において、前記固体・液体が接する面として定義される固液界面に厚みを持たせた固液界面領域を定義するステップと、離散流速の定義点を中心とする胞を定義するステップと、前記胞内における前記固液界面領域の体積を計算するステップと、前記胞における離散流速を前記固液界面の法線ベクトルと直交するように射影するステップと、前記射影された離散流速および前記離散流速の定義点と前記固液界面との距離および前記固液界面領域の体積をもって前記胞における壁粘性応力の平均値を計算するステップと、前記壁粘性応力の平均値を前記離散流速の定義点における圧力に変換するステップとを有することを特徴とする固液界面における壁粘性応力計算方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１０−１９７２８７（Ｐ２０１０−１９７２８７Ａ）
【公開日】平成２２年９月９日（２０１０．９．９）
【国際特許分類】
【出願番号】特願２００９−４３９３０（Ｐ２００９−４３９３０）
【出願日】平成２１年２月２６日（２００９．２．２６）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【公開日】平成２２年９月９日（２０１０．９．９）
【国際特許分類】
【出願日】平成２１年２月２６日（２００９．２．２６）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
[ Back to top ]