特許 7386481 気液二相流のシミュレーション方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2023-11-16

(45)【発行日】2023-11-27

(54)【発明の名称】気液二相流のシミュレーション方法

(51)【国際特許分類】

   G01N 13/02 20060101AFI20231117BHJP

【ＦＩ】

G01N13/02

【請求項の数】 5

(21)【出願番号】P 2020133864

(22)【出願日】2020-08-06

(65)【公開番号】P2022030097

(43)【公開日】2022-02-18

【審査請求日】2023-01-27

(73)【特許権者】

【識別番号】504139662

【氏名又は名称】国立大学法人東海国立大学機構

(73)【特許権者】

【識別番号】501089759

【氏名又は名称】デンソーテクノ株式会社

(74)【代理人】

【識別番号】110000578

【氏名又は名称】名古屋国際弁理士法人

(72)【発明者】

【氏名】辻 義之

(72)【発明者】

【氏名】伊藤 高啓

(72)【発明者】

【氏名】佐藤 隆哉

(72)【発明者】

【氏名】大山 武士

(72)【発明者】

【氏名】甲村 圭司

【審査官】北条 弥作子

(56)【参考文献】

【文献】特開２０１７－１３８１２３（ＪＰ，Ａ）

【文献】特開２０１２－２１２８３３（ＪＰ，Ａ）

【文献】米国特許出願公開第２０１２／０２５３７６７（ＵＳ，Ａ１）

【文献】佐藤 隆哉，せん断流中の壁面上水滴挙動の動的接触角のモデル化，自動車技術会論文集，2019年05月，Vol. 50, No. 3，pp. 926-931

【文献】１０．境界層と物体まわりの流れ，流体の力学 基礎編，2019年09月13日，pp. 1-7

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｎ １１／００～１３／０４

ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

固体壁面に付着している液滴に対して気体の流れによる外力が働いていない状態における、前記固体壁面に対する前記液滴の接触角を用いて、前記液滴及び前記気体それぞれの運動方程式により、前記液滴と前記気体との間での運動量の交換則に基づいて気液二相流の質量流量を算出する、質量流量算出ステップ（３１，Ｓ１０３）と、
前記質量流量算出ステップにて算出された前記質量流量、及び前記それぞれの運動方程式により算出された前記液滴と前記気体との速度差を用いて、前記外力により生じる、前記液滴と前記気体との界面に働く抗力を算出する、抗力算出ステップ（３１，Ｓ１０４）と、
あらかじめ導出された、前記接触角と前記抗力との関係式により、前記抗力算出ステップにて算出された前記抗力を用いて、前記液滴に対して前記外力が働いている状態における前記接触角を算出する、接触角算出ステップ（３１，Ｓ１０５）と、
を備える、気液二相流のシミュレーション方法。

【請求項2】

請求項１に記載の気液二相流のシミュレーション方法であって、
前記質量流量算出ステップでは、前記液滴及び前記気体を含む計算領域が複数に分割された計算要素ごとに、前記計算要素における気液二相流の運動量交換係数と前記計算要素の体積との積を前記質量流量として算出し、
前記運動量交換係数は、前記液滴と前記気体との体積比率に基づいて算出された前記計算要素における気液二相流の密度及び粘度を用いて算出される、気液二相流のシミュレーション方法。

【請求項3】

請求項１又は請求項２に記載の気液二相流のシミュレーション方法であって、
前記関係式は、
前記気体の速度ベクトルが示す方向に対して前側で成す前記接触角についての前記関係式である前進関係式と、
前記方向に対して後ろ側で成す前記接触角についての前記関係式である後退関係式と、
を含み、
前記接触角算出ステップは、前記前側に位置している前記界面について、前記前進関係式により、前記液滴に対して前記外力が働いている状態における前記接触角を算出し、前記後ろ側に位置している前記界面について、前記後退関係式により、前記液滴に対して前記外力が働いている状態における前記接触角を算出する、気液二相流のシミュレーション方法。

【請求項4】

請求項３に記載の気液二相流のシミュレーション方法であって、
前記接触角算出ステップでは、前記気体の速度ベクトルと前記界面の法線ベクトルとの成す角に基づいて、前記界面が前記前側に位置しているか前記後ろ側に位置しているかを判定する、気液二相流のシミュレーション方法。

【請求項5】

請求項３又は請求項４に記載の気液二相流のシミュレーション方法であって、
前記前進関係式は、前記液滴に対して前記外力が働いている状態における前記接触角のうちの最も大きい前記接触角と前記抗力との前記関係式であり、
前記後退関係式は、前記液滴に対して前記外力が働いている状態における前記接触角のうちの最も小さい前記接触角と前記抗力との前記関係式である、気液二相流のシミュレーション方法。

【発明の詳細な説明】

【技術分野】

【0001】

本開示は、気液二相流のシミュレーション方法に関する。

【背景技術】

【0002】

固体壁面に付着した液滴に対して気体の流れによる外力が働いている状態においては、外力が臨界値を超えるまで、液滴は気液界面の形状を変化させながら、すなわち固体壁面に対する接触角を変化させながら固体壁面上で移動を開始しないことが、実験的に確認されている。しかし、従来の気液二相流のシミュレーション方法では、液滴の接触角として固定値を用いていたため、上記外力が臨界値を超えるまで液滴は移動を開始しないという現象をシミュレーションで再現することが難しかった。

【0003】

また、気液界面の形状が変化している状態における液滴の接触角に関しては、例えば非特許文献１に接触角の予測式が開示されているが、この予測式は、固体壁面上を移動している状態における液滴の接触角に注目して算出されたものである。つまり、非特許文献１に記載の接触角の予測式は、上記のような固体壁面上で移動を開始する前後における液滴の接触角について注目して算出されたものではなかった。

【0004】

これらに対し、特許文献１では、固体壁面上で移動を開始する前後における液滴の接触角に注目し、気液界面の速度に基づいて液滴の接触角条件を変更する、気液二相流のシミュレーション方法が開示されている。

【先行技術文献】

【特許文献】

【0005】

【文献】特開２０１７－１３８１２３号公報

【非特許文献】

【0006】

【文献】Kensuke Yokoi, Damien Vadillo, “Numerical studies of the influence of the dynamic contact angle on a droplet impacting on a dry surface”, Physics of Fluids, (米), 2009, Vol. 21, Issue 7, p. 072102

【発明の概要】

【発明が解決しようとする課題】

【0007】

しかしながら、発明者の詳細な検討の結果、特許文献１に記載された気液二相流のシミュレーション方法では、初期接触角と、気液界面の速度の関数として与えられる随時接触角と、の差が所定の閾値以下となるまで収束計算が必要であるため、多くの計算時間を要するという課題が見出された。

【0008】

本開示の一局面は、固体壁面に付着している液滴に対して気体の流れによる外力が働き、固体壁面に対する液滴の接触角が変化する場合における、気液二相流のシミュレーションにおいて、計算時間を短縮するための技術を提供する。

【課題を解決するための手段】

【0009】

本開示の一態様は、気液二相流のシミュレーション方法であって、質量流量算出ステップ（３１，Ｓ１０３）と、抗力算出ステップ（３１，Ｓ１０４）と、接触角算出ステップ（３１，Ｓ１０５）と、を備える。質量流量算出ステップでは、固体壁面に付着している液滴に対して気体の流れによる外力が働いていない状態における、固体壁面に対する液滴の接触角を用いて、液滴及び気体それぞれの運動方程式により、液滴と気体との間での運動量の交換則に基づいて気液二相流の質量流量を算出する。抗力算出ステップでは、質量流量算出ステップにて算出された質量流量、及びそれぞれの運動方程式により算出された液滴と気体との速度差を用いて、外力により生じる、液滴と気体との界面に働く抗力を算出する。接触角算出ステップでは、あらかじめ導出された、接触角と抗力との関係式により、抗力算出ステップにて算出された抗力を用いて、液滴に対して外力が働いている状態における接触角を算出する。

【0010】

このような構成によれば、固体壁面に付着している液滴に対して気体の流れによる外力が働き、固体壁面に対する液滴の接触角が変化する場合における、気液二相流のシミュレーションにおいて、計算時間を短縮することができる。

【図面の簡単な説明】

【0011】

【図1】シミュレーション装置のブロック図である。

【図2】固体壁面に付着している液滴が気体による外力を受けている状態を示す模式図である。

【図3】液滴の高さと排除厚さとを比較したグラフである。

【図4】気体の流れによる外力によって生じる、気液界面に働く抗力を示す模式図である。

【図5】数値計算により得られた、排除厚さに対する液滴の初期高さの比率と抗力係数との関係を示すグラフである。

【図6】固体壁面に対する液滴の接触角のモデル化により得られた、液滴の動的接触角に関するグラフである。

【図7】シミュレーション装置がシミュレーションプログラムに従って実行する処理を示すフローチャートである。

【図8】固体壁面に付着している液滴が固体壁面上で移動を開始するときの気体の速度について、シミュレーション結果と実験結果とを比較したグラフである。

【発明を実施するための形態】

【0012】

以下、本開示の例示的な実施形態について図面を参照しながら説明する。
［１．構成］
図１に示すシミュレーション装置１は、図２に示すような、固体壁面Ｗに付着している液滴Ｌに対して気体Ｇの流れによる外力が働き、液滴と気体との界面（以下、気液界面ともいう。）の形状が変化している場合における、気液二相流のシミュレーションを実行する装置である。シミュレーション装置１は、例えば、車載カメラのレンズに付着した雨水等の水滴を気流により除去する場合や、熱交換器内で凝縮した水滴を気流により排出する場合等における、気液二相流のシミュレーションを実行するために用いられる。

【0013】

シミュレーション装置１は、ＡＮＳＹＳ社の汎用流体解析ソフトウエアＦＬＵＥＮＴに対して、本開示の機能をユーザ定義関数ＵＤＦ（User Defined Function）として組み込むことにより実現した。ただし、本開示が適用されるソフトウエアは、ＦＬＵＥＮＴに限定されるものではない。

【0014】

図１に示すように、シミュレーション装置１は、入力部１０、表示部２０、及び制御部３０を備える。
入力部１０は、シミュレーションの条件等が入力されるデバイスである。入力部１０には、例えば、キーボードやマウス等が該当する。

【0015】

表示部２０は、シミュレーションの結果等が表示されるデバイスである。表示部２０には、例えば、ディスプレイ等が該当する。
制御部３０は、ＣＰＵ３１、ＲＡＭ３２、ＲＯＭ３３、及びＨＤＤ３４を備える。ＨＤＤ３４には、シミュレーションを実行するためのシミュレーションプログラムが格納されている。シミュレーション装置１において、入力部１０への入力に応じてシミュレーションプログラムが起動されると、ＣＰＵ３１がシミュレーションプログラムをＲＡＭ３２に読み出して実行する。これにより、当該シミュレーションプログラムに対応するシミュレーション方法が実行され、気液二相流のシミュレーションが実行される。そして、その結果が表示部２０に表示される。

【0016】

シミュレーション装置１は、境界面の解析方法として公知の技術であるＶＯＦ（Volume of Fluid）法を採用する。この際、シミュレーション装置１は、後述するとおり、液滴と気体との間での運動量の交換則に基づいて気液二相流のシミュレーションを実行するため、気相及び液相それぞれの運動方程式を計算するモデルを採用する。ＦＬＵＥＮＴでは、multi-Fluid ＶＯＦ法が当該モデルを採用したＶＯＦ法に該当する。

【0017】

［２．動的接触角についての関係式の導出］
シミュレーション装置１は、後述するように、ＦＬＵＥＮＴのＵＤＦとして組み込まれた、液滴の動的接触角についての関係式を用いて、気液二相流のシミュレーションを実行する。本開示における液滴の動的接触角とは、固体壁面に付着している液滴に対して気体の流れによる外力が働き、気液界面の形状が変化している状態における、固体壁面に対する液滴の接触角を意味する。

【0018】

そこで、まず、上記関係式の導出を目的とした、液滴の動的接触角のモデル化について説明する。なお、当該モデル化について、本発明者らは、研究論文（“せん断流中の壁面上水滴挙動の動的接触角のモデル化”，自動車技術会論文集，公益社団法人自動車技術会，2019年3月，第50巻，第3号，p. 926－931）にて詳細を報告している。

【0019】

本開示におけるシミュレーションで注目している、固体壁面に付着している液滴は、図３に示すように、十分に発達した気体の境界層δ内に位置する。したがって、液滴が気体から受ける、気体の流れによる外力、すなわち気体の運動エネルギーに基づく外力は、気体の中心速度Ｕ_０ではなく境界層δ内における気体の速度に基づいて考える必要がある。

【0020】

気体の境界層δを特徴づけるパラメータとして、式（１）で表される排除厚さδ^＊が挙げられる。式（１）におけるＵ（ｙ）は、固体壁面からｙ［ｍｍ］離れた位置における気体の速度である。

【0021】

式（２）で表される、排除厚さδ^＊に対する液滴の初期高さｈ_０の比率は、液滴が気体の流れの影響をどの程度受けるかを示すパラメータであると考えられる。液滴の初期高さｈ_０とは、液滴に対して気体の流れによる外力が働いていない状態における、固体壁面からの液滴の高さである。

【0022】

【数1】

【0023】

気体の流れによる外力によって生じる、気液界面に働く抗力Ｆ_Ｄは、図４に示すように、気液界面にかかる圧力Ｐ及びせん断応力τのそれぞれにおける気体の速度ベクトルが示す方向の成分を表面積分する数値計算により、算出することができる。数値計算に用いられるソフトとしては、例えば、ＳＩＥＭＥＮＳ社の汎用熱流体解析ソフトSimcenter ＳＴＡＲ－ＣＣＭ＋が挙げられる。しかし、このように算出される抗力は、気体の速度や液滴の大きさに依存する。そこで、本実施形態では、汎用性の観点から、上記抗力を無次元化した抗力係数を定義する。具体的には、気体の運動エネルギー（１／２ρＵ_０ ^２）が気液界面の投影面積Ｓに作用した場合の力（１／２ρＵ_０ ^２Ｓ）に対する、気液界面に働く抗力Ｆ_Ｄの比率として、式（３）で表される抗力係数Ｃ_Ｄを定義する。式（３）では、投影面積Ｓが液滴の初期高さｈ_０と単位幅長さ１との積で表されている。なお、ここでいう気液界面の投影面積Ｓとは、気体の速度ベクトルが示す方向に沿って見たときの気液界面の面積を指す。また、式（３）におけるρは、気体の密度である。

【0024】

【数2】

【0025】

本発明者らが、Simcenter ＳＴＡＲ－ＣＣＭ＋を用いて、液滴の初期高さｈ_０及び気体の中心速度Ｕ_０の条件を複数変更して二次元定常解析を行ったところ、図５に示すように、式（２）と式（３）との間には式（４）で表される関係式が成立することが判明した。式（４）におけるｆは、関数を意味する。

【0026】

【数3】

【0027】

ここで、S. R. Annapragadaらの報告（“Droplet Retention on an Incline”, International Journal of Heat and Mass Transfer, (蘭), 2012, Vol. 55, p. 1457－1465）によれば、傾斜した固体壁面上の液滴が重力により滑落している状態、すなわち傾斜した固体壁面上を液滴が移動している状態における液滴の接触角θは、式（５）で表される。式（５）におけるβ、Ｂｏ、ρ、ｇ、Ｄ、及びσは、それぞれ、水平に対する固体壁面の傾斜角度、ボンド数、液滴の密度、重力加速度、液滴の直径、及び気液界面の界面張力である。また、式（５）におけるｆは、関数を意味する。

【0028】

【数4】

【0029】

本発明者らは、式（５）におけるボンド数Ｂｏの分子（ρｇＤ^２sinβ）を重力による推進力、上述した抗力Ｆ_Ｄを気体の流れによる推進力と捉え、式（５）におけるボンド数Ｂｏの分子を抗力Ｆ_Ｄに置換可能であると考えた。その結果、式（４）及び式（５）より式（６）が得られ、液滴の動的接触角のモデル化が達成された。式（６）で表されるように、固体壁面に対する液滴の接触角θは、気液界面の界面張力σに対する抗力Ｆ_Ｄの比率の関数として与えられる。

【0030】

【数5】

【0031】

図６に、式（６）及び実験にて得られた計測値を用いて作成された、液滴の動的接触角θ_ｄに関するグラフを示す。計測値とは、具体的には、液滴の初期接触角θ_０、動的接触角θ_ｄ、及び初期高さｈ_０、並びに気体の中心速度Ｕ_０の数値である。実験条件を限定するものではないが、図６に示すグラフでは、アクリル板上の水滴に対して所定の速度の空気を作用させる実験にて得られた計測値が用いられている。

【0032】

液滴の動的接触角には、液滴の前進接触角及び後退接触角が含まれる。本開示における液滴の前進接触角とは、気体の速度ベクトルが示す方向に対して前側で成す液滴の動的接触角を指す。本開示における液滴の後退接触角とは、気体の速度ベクトルが示す方向に対して後ろ側で成す液滴の動的接触角を指す。図６において、四角形のプロットは、動的接触角として前進接触角を用いたものであり、三角形のプロットは、動的接触角として後退接触角を用いたものである。また、四角形又は三角形のプロットのうち、白抜きのプロットは、実験にて液滴が固体壁面上で移動を開始したことが観察されるまでの測定値を用いたものであり、黒塗りのプロットは、液滴が固体壁面上で移動を開始したことが観察された後の測定値を用いたものである。

【0033】

式（６）を用いた、液滴の動的接触角に関するグラフの作成にあたっては、液滴におけるいずれの位置で計測された動的接触角が用いられてもよい。図６に示すグラフの作成にあたっては、気体の速度ベクトルが示す方向に対して、最も前側で成す液滴の動的接触角、及び、最も後ろ側で成す液滴の動的接触角が用いられている。換言すれば、動的接触角のうちの最も大きい動的接触角及び最も小さい動的接触角が用いられている。

【0034】

図６に示すグラフを基に、前進接触角θ_１及び後退接触角θ_２のそれぞれについて、式（７）で表される近似式が導出される。式（７）は、液滴の動的接触角についての関係式に該当する。具体的には、式（７）は、液滴の動的接触角θ_ｄと抗力Ｆ_Ｄとの関係式である。当該関係式には、関数ｆ_１で表される、前進接触角θ_１についての関係式である前進関係式と、関数ｆ_２で表される、後退接触角θ_２についての関係式である後退関係式と、が含まれる。

【0035】

【数6】

【0036】

導出された前進関係式及び後退関係式のそれぞれは、ＦＬＵＥＮＴのＵＤＦとしてシミュレーション装置１に組み込まれる。すなわち、本実施形態では、式（７）がＦＬＵＥＮＴのＵＤＦとしてシミュレーション装置１に組み込まれる。

【0037】

［３．処理］
次に、ＨＤＤ３４に格納されたシミュレーションプログラムに従ってＣＰＵ３１が実行する具体的な処理について、図７を用いて説明する。

【0038】

まず、Ｓ１０１で、ＣＰＵ３１は、ＦＬＵＥＮＴであらかじめ定められた計算要素ごとに初期条件を設定する。計算要素とは、気体及び固体壁面に付着している液滴を含む計算領域が複数個に分割された、複数の領域のそれぞれを指す。初期条件として設定されるデータには、各計算要素に占める第１流体の体積率を表すＶＯＦ値Ｎ_VOF、初期接触角、液滴及び気体の物性値等がある。

【0039】

ＶＯＦ値Ｎ_VOFは、０～１の値を取る。ＶＯＦ値Ｎ_VOFは、Ｎ_VOF＝１のときに、その要素に第１流体のみが存在することを示し、Ｎ_VOF＝０のときに、その要素に第２流体のみが存在することを示す。本実施形態では、第１流体を気体、第２流体を液滴とする。

【0040】

また、本実施形態では、初期接触角として静止接触角を設定する。静止接触角とは、固体壁面に付着している液滴に対して気体の流れによる外力が働いていない状態における、固体壁面に対する液滴の接触角である。

【0041】

続いて、Ｓ１０２で、ＣＰＵ３１は、計算時間ｔを０に設定する。
続いて、Ｓ１０３で、ＣＰＵ３１は、Ｓ１０１で設定された初期条件に基づいて、液滴及び気体それぞれの運動方程式により、気液界面を包含する計算要素（以下、界面要素ともいう。）ごとに、液滴と気体との間での運動量の交換則に基づいて質量流量を算出する。界面要素は、液滴及び気体の双方が存在する計算要素とも言い換えることができる。界面要素におけるＶＯＦ値Ｎ_VOFは、０よりも大きく１未満である。

【0042】

具体的には、ＣＰＵ３１は、まず、ＦＬＵＥＮＴのmulti-Fluid ＶＯＦ法による液滴及び気体の二相流計算を実行し、式（８）で表される液滴の運動方程式、及び式（９）で表される気体の運動方程式をそれぞれ計算する。式（８）及び式（９）におけるα、ρ、ｕ、ｐ、τ及びｇは、それぞれ、体積割合、密度、速度、圧力、せん断応力及び重力加速度である。式（８）及び式（９）における添え字ｄ及びａは、それぞれ液滴及び気体を意味する。また、式（８）におけるＦは、液滴に働く外力の和を表す外力項である。ここでいう外力には、上述した抗力に加え、重力や揚力なども含まれる。

【0043】

【数7】

【0044】

次に、ＣＰＵ３１は、界面要素ごとに、界面要素における気液二相流の運動量交換係数を算出する。当該運動量交換係数の算出にあたっては、Schiller-Naumannの抗力モデルを基に、液滴と気体との体積比率に基づいて算出された界面要素における気液二相流の密度及び粘度が用いられる。そして、ＣＰＵ３１は、界面要素ごとに、界面要素における上記運動量交換係数と、界面要素の体積と、の積を質量流量として算出する。

【0045】

また、ＣＰＵ３１は、各界面要素について、液滴の運動方程式から液滴の速度を、気体の運動方程式から気体の速度を、それぞれ算出する。なお、ＣＰＵ３１は、液滴及び気体それぞれの運動方程式により、各計算要素における速度分布や圧力分布等も算出可能である。

【0046】

続いて、Ｓ１０４で、ＣＰＵ３１は、ＦＬＵＥＮＴのＵＤＦにより、気体の流れによる外力によって生じる、気液界面に働く抗力を算出する。具体的には、ＣＰＵ３１は、まず、界面要素ごとに、Ｓ１０３で算出された質量流量と、液滴と気体との速度差と、の積を界面要素における上記抗力として算出する。そして、ＣＰＵ３１は、複数の界面要素それぞれにおける抗力の総和を、気液界面全体に働く抗力として算出する。

【0047】

続いて、Ｓ１０５で、ＣＰＵ３１は、ＦＬＵＥＮＴのＵＤＦにより、Ｓ１０４で算出された抗力を用いて液滴の動的接触角を算出する。
具体的には、まず、ＣＰＵ３１は、固体壁面から第１層目の界面要素、換言すれば固体壁面に接する界面要素、更に換言すれば三重線を包含する界面要素を抽出する。三重線とは、固体壁面、液滴及び気体の界面が形成する線である。以下、三重線を包含する界面要素を三重線要素ともいう。

【0048】

そして、ＣＰＵ３１は、三重線要素ごとに、ＶＯＦ値Ｎ_VOFの勾配▽Ｎ_VOFと気体の速度ベクトルとの成す角に基づいて、気液界面が前進しているか後退しているかを判定する。ここでいう前進している気液界面とは、固体壁面に付着している液滴に対して気体の流れによる外力が働いている状態において、気体の速度ベクトルが示す方向に対して前側に位置している気液界面を指す。後退している気液界面とは、上記状態において、気体の速度ベクトルが示す方向に対して後ろ側に位置している気液界面を指す。つまり、ＣＰＵ３１は、三重線要素ごとに、気体の速度ベクトルが示す方向に対して、気液界面が前側に位置しているか後ろ側に位置しているかを判定する。

【0049】

ここで、▽Ｎ_VOFは、気液界面に垂直、且つ、液滴から気体へと向かう方向のベクトルである。そこで、本実施形態では、ＣＰＵ３１は、▽Ｎ_VOFと気体の速度ベクトルとが成す角が９０度未満である場合に、気液界面が前進していると判定し、上記成す角が９０度よりも大きい場合に、気液界面が後退していると判定する。上記成す角が９０度である場合には、ＣＰＵ３１は、気液界面が前進も後退もしていない、すなわち平衡状態であると判定する。判定にあたっては、式（１０）に示すように、▽Ｎ_VOFと気体の速度ベクトルとの内積が用いられる。式（１０）におけるｖは、気体の速度ベクトルである。

【0050】

【数8】

【0051】

次に、ＣＰＵ３１は、Ｓ１０４で算出された気液界面全体に働く抗力を用いて、ＦＬＵＥＮＴのＵＤＦに組み込まれた、動的接触角と抗力との関係式により、動的接触角を算出する。上述したとおり、本実施形態では、動的接触角と抗力との関係式として、式（７）で表される前進関係式及び後退関係式がＵＤＦに組み込まれている。したがって、式（７）における前進関係式及び後退関係式のそれぞれに対して、抗力Ｆ_Ｄの値として、Ｓ１０４で算出された気液界面全体に働く抗力を代入することにより、ＣＰＵ３１は前進接触角及び後退接触角を算出する。

【0052】

続いて、Ｓ１０６で、ＣＰＵ３１は、三重線要素ごとに接触角条件を変更する。具体的には、Ｓ１０５で気液界面が前進していると判定した三重線要素について、ＣＰＵ３１は、接触角条件をＳ１０５で算出した前進接触角に変更する。また、Ｓ１０５で気液界面が後退していると判定した三重線要素について、ＣＰＵ３１は、接触角条件をＳ１０５で算出した後退接触角に変更する。Ｓ１０５で気液界面が平衡状態であると判定した三重線要素については、ＣＰＵ３１は、接触角条件を変更せず、Ｓ１０１で設定した初期接触角を維持する。

【0053】

続いて、Ｓ１０７で、ＣＰＵ３１は、計算時間ｔがあらかじめ定められた終了時間に達したかを判定する。後述するように、計算時間ｔが終了時間に達していないと判定した場合、ＣＰＵ３１は、次のタイムステップについて、Ｓ１０３以降の処理を実行する。換言すれば、ＣＰＵ３１は、各タイムステップについて、Ｓ１０３以降の処理を１回ずつ実行する。そこで、例えば、シミュレーション装置１のユーザが液滴の挙動を継続して確認したい時間幅を基に、終了時間は定められる。例えば、シミュレーション装置１のユーザが１０秒間にわたって液滴の挙動を確認したい場合、終了時間は１０秒に定められる。

【0054】

Ｓ１０７で、計算時間ｔが終了時間に達していないと判定した場合、ＣＰＵ３１は、Ｓ１０８に移行し、計算時間ｔをｔ＋Δｔに設定する。Δｔは、タイムステップである。そして、ＣＰＵ３１は、本処理をＳ１０３に戻す。

【0055】

一方、Ｓ１０７で、計算時間ｔが終了時間に達したと判定した場合、ＣＰＵ３１は、本処理を終了する。
［４．効果］
以上詳述した実施形態によれば、以下の効果が得られる。

【0056】

（４ａ）上記実施形態に係る気液二相流のシミュレーション方法では、あらかじめ導出された、固体壁面に対する液滴の接触角と抗力との関係式を用いる。ここでいう抗力は、気体の流れによる外力により生じる、気液界面に働く抗力である。

【0057】

したがって、固体壁面に付着している液滴に対して気体の流れによる外力が働き、固体壁面に対する液滴の接触角が変化する場合における、気液二相流のシミュレーションにおいて、計算時間を短縮することができる。

【0058】

また、上記関係式を用いることにより、上記抗力によって液滴の接触角が変化する状態をシミュレーション上で再現することができる。これにより、液滴の接触角が変化することによって生じる気液界面の界面張力の抵抗成分、すなわち気液界面の界面張力における気体の速度ベクトルが示す方向の成分が、上記抗力に抵抗する状態を再現することができる。したがって、固体壁面に付着している液滴に対して気体の流れによる外力が働き、固体壁面に対する液滴の接触角が変化している状態において、固体壁面上で液滴が移動を開始しない状態を再現することができる。

【0059】

さらに、上記実施形態に係る気液二相流のシミュレーション方法では、液滴及び気体それぞれの運動方程式を用いるため、液滴及び気体を１つの混合流体として扱う１つの運動方程式を用いる場合と比較して、液滴と気体との速度差をより正確に算出することができ、上記抗力をより正確に算出することができる。したがって、比較的高いシミュレーション精度を得ることができる。

【0060】

図８は、上記実施形態に係る気液二相流のシミュレーション方法によるシミュレーション結果と実験結果とを比較したグラフである。図８では、固体壁面に付着している液滴に対して気体の流れによる外力が働いている状態において、固体壁面上で液滴が移動を開始するときの気体の速度について、シミュレーション結果と実験結果とを比較している。また、図８には、２つの従来手法によるシミュレーション結果も併せて示す。ここでいう２つの従来手法、すなわち図８における従来手法Ａ及び従来手法Ｂは、それぞれ、固体壁面に対する液滴の接触角として固定値を用いる手法、及び、固体壁面上を移動している状態における液滴の接触角に注目して算出された接触角の予測式（非特許文献１参照）を用いる手法である。

【0061】

図８に示すように、従来手法Ａ及び従来手法Ｂでは、気体の速度が遅い条件において液滴が移動を開始してしまうシミュレーション結果となり、実験結果との乖離が大きかった。これに対して、本実施形態に係る気液二相流のシミュレーション方法では、従来手法Ａ及び従来手法Ｂと比較して、実験結果との乖離が小さいシミュレーション結果を得た。すなわち、本実施形態に係る気液二相流のシミュレーション方法では、固体壁面上で液滴が移動を開始しない状態を再現することができ、比較的高いシミュレーション精度を得ることができる。

【0062】

（４ｂ）上記実施形態に係る気液二相流のシミュレーション方法では、計算要素ごとに質量流量を算出する。したがって、質量流量の算出に要する計算時間を短縮することができる。

【0063】

（４ｃ）上記実施形態に係る気液二相流のシミュレーション方法では、液滴の接触角と上記抗力との関係式は、前進関係式と後退関係式とを含む。そして、気体の速度ベクトルが示す方向に対して前側に位置している気液界面について、前進関係式により液滴の接触角を算出し、後ろ側に位置している気液界面について、後退関係式により液滴の接触角を算出する。

【0064】

例えば、固体壁面に付着している液滴に対して気体の流れによる外力が働いている状態では、固体壁面上で液滴が移動を開始するまで、気体の速度ベクトルに対して前側に位置している気液界面では、固体壁面に対する液滴の接触角が徐々に大きくなり、後ろ側に位置している気液界面では、固体壁面に対する液滴の接触角が徐々に小さくなる。つまり、気液界面の位置によって、液滴の接触角の変化が異なる。このため、気体の速度ベクトルに対して前側及び後ろ側それぞれについて異なる関係式、具体的には前進関係式又は後退関係式を用いることとすれば、これら関係式は、比較的簡単な計算式で記述することができる。したがって、計算時間をより短縮することができる。

【0065】

また、例えば、前進関係式又は後退関係式のいずれか一方のみを用いて、液滴の接触角を算出することも考えられるところ、上記実施形態に係る気液二相流のシミュレーション方法では、前進関係式及び後退関係式のそれぞれを用いて液滴の接触角を算出する。したがって、より高いシミュレーション精度を得ることができる。

【0066】

（４ｄ）上記実施形態に係る気液二相流のシミュレーション方法では、気体の速度ベクトルと気液界面の法線ベクトルとの成す角に基づいて、気体の速度ベクトルが示す方向に対して気液界面が前側に位置しているか後ろ側に位置しているかを判定する。

【0067】

ここで、気体の速度ベクトルが示す方向に対して気液界面が前側に位置しているか後ろ側に位置しているかの判定方法としては、例えば、気液界面の速度ベクトルと気液界面の法線ベクトルとの成す角に基づく判定方法も考えられる。しかし、この判定方法では、気液界面の速度ベクトルが０である場合、すなわち気液界面の速さが０である場合、気液界面の速度ベクトルと気液界面の法線ベクトルとの成す角が０となってしまう。このため、上記判定方法では、気液界面の速さが０である場合に、気体の速度ベクトルが示す方向に対して気液界面が前側に位置しているか後ろ側に位置しているかを適切に判定することが難しい。

【0068】

これに対して、上記実施形態に係る気液二相流のシミュレーション方法では、気液界面の速さによらず、気体の速度ベクトルが示す方向に対して気液界面が前側に位置しているか後ろ側に位置しているかを適切に判定することができる。

【0069】

（４ｅ）上記実施形態に係る気液二相流のシミュレーション方法では、前進関係式は、最も大きい動的接触角と上記抗力との関係式であり、後退関係式は、最も小さい動的接触角と上記抗力との関係式である。すなわち、算出された上記抗力を用いて、最も大きい動的接触角及び最も小さい動的接触角を算出する。

【0070】

ここで、固体壁面に付着している液滴に対して気体の流れによる外力が働いている状態においては、最も大きい動的接触角と最も小さい動的接触角との差が所定の閾値を上回ったときに固体壁面上で液滴は移動を開始することが、実験的に確認されている。このため、最も大きい動的接触角及び最も小さい動的接触角を算出すれば、これら以外の動的接触角については算出しなくても、固体壁面上で移動を開始する前後における液滴の挙動のシミュレーションを行うことは可能である。上記実施形態に係る気液二相流のシミュレーション方法では、最も大きい動的接触角及び最も小さい動的接触角以外の動的接触角についての計算を削減することができ、計算時間をより短縮することができる。

【0071】

なお、上記実施形態では、Ｓ１０３が質量流量算出ステップに相当し、Ｓ１０４が抗力算出ステップに相当し、Ｓ１０５が接触角算出ステップに相当する。
［５．他の実施形態］
以上、本開示の実施形態について説明したが、本開示は、上記実施形態に限定されることなく、種々の形態を採り得ることは言うまでもない。

【0072】

（５ａ）上記実施形態では、図６に示す動的接触角に関するグラフを基に、液滴の動的接触角についての関係式として、式（７）で表される２つの近似式を算出するが、例えば、当該グラフに含まれるプロットに対して１つの近似式を算出してもよい。この場合、上述したＳ１０５で、Ｓ１０４で算出された抗力を当該１つの近似式に代入することにより２つの動的接触角が算出され得るところ、例えば、上述したＳ１０６で、Ｓ１０５における気液界面が前進しているか後退しているかの判定結果に基づいて、いずれの動的接触角を接触角条件として採用するかを決定してもよい。

【0073】

（５ｂ）上記実施形態では、上述したＳ１０６で、Ｓ１０５で前進していると判定した三重線要素について、接触角条件をＳ１０５で算出した前進接触角に変更し、Ｓ１０５で後退していると判定した三重線要素について、接触角条件をＳ１０５で算出した後退接触角に変更する。Ｓ１０５で算出される前進接触角及び後退接触角は、上記（４ｅ）で述べたとおり、最も大きい動的接触角及び最も小さい動的接触角であるが、接触角条件として採用される接触角は、最も大きい動的接触角及び最も小さい動的接触角に限定されるものではない。例えば、Ｓ１０５で算出された最も大きい動的接触角及び最も小さい動的接触角のそれぞれを用いて、補間計算により、これら以外の動的接触角を１つ以上算出し、▽Ｎ_VOFと気体の速度ベクトルとが成す角に基づいて、算出された動的接触角のいずれかを接触角条件として採用してもよい。

【0074】

（５ｃ）上記実施形態における１つの構成要素が有する機能を複数の構成要素として分散させたり、複数の構成要素が有する機能を１つの構成要素に統合したりしてもよい。また、上記実施形態の構成の一部を省略してもよい。また、上記実施形態の構成の少なくとも一部を、他の上記実施形態の構成に対して付加、置換等してもよい。

【0075】

（５ｄ）本開示は、上述した気液二相流のシミュレーション方法の他、上述した気液二相流のシミュレーションを実行するためのプログラム、このプログラムを記録した媒体及びシミュレーション装置１など、種々の形態で実現することができる。

【符号の説明】

【0076】

１…シミュレーション装置、１０…入力部、２０…表示部、３０…制御部、３１…ＣＰＵ、３２…ＲＡＭ、３３…ＲＯＭ、３４…ＨＤＤ、Ｗ…固体壁面、Ｌ…液滴、Ｇ…気体。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】