(57)【要約】【課題】、固体壁に付着した液滴が気体から外力を受け、液滴の固体壁に対する接触角が変化する場合の気液二相流の挙動のシミュレーションにおいて、シミュレーション精度を向上する。【解決手段】気液二相流のシミュレーション方法であって、固相に付着する液相の、固相に対する初期接触角を設定するステップと、初期接触角を用いて、界面速度を算出するステップと、固相に対する液相の随時接触角を、所定の条件を満たす、界面速度の関数として算出するステップと、初期接触角と随時接触角との差が閾値以下の場合、計算が収束したと判定するステップと、を有する。所定の条件は、界面の前進時、随時接触角は、液相が移動を開始するときの第１臨界角以下の範囲で増加し、界面の後退時、随時接触角は、液相が移動を開始するときの第２臨界角以上の範囲で減少し、界面速度が０では、随時接触角は初期接触角から変化しない、との条件である。【選択図】図３
2017138123.tif 000002
【特許請求の範囲】
【請求項１】
気液二相流のシミュレーション方法であって、
固相に付着する液相の、前記固相に対する初期接触角を設定するステップと、
前記初期接触角を用いて、前記固相に接触する前記液相と気相との界面における前記液相の速度である界面速度を算出する、第一の算出ステップと、
前記固相に対する前記液相の随時接触角を、所定の条件を満たす、前記界面速度の関数として算出する、第二の算出ステップと、
前記初期接触角と前記随時接触角との差が閾値以下の場合、計算が収束したと判定するステップと、
を有し、
前記所定の条件は、
前記界面が前進しているとき、前記随時接触角は、前記液相が移動を開始するときの前記固相に対する前記液相の移動方向前側の接触角である第１臨界角以下の範囲で増加し、
前記界面が後退しているとき、前記随時接触角は、前記液相が移動を開始するときの前記固相に対する前記液相の前記移動方向後側の接触角である第２臨界角以上の範囲で減少し、
前記界面速度が０の場合、前記随時接触角は前記初期接触角から変化しない、
との条件であり、
前記判定するステップにおいて、前記初期接触角と前記随時接触角との差が閾値よりも大きい場合、前記設定するステップにおいて前記随時接触角を前記初期接触角に再設定し、前記第一の算出ステップにおいて前記界面速度を再計算し、前記第二の算出ステップにおいて前記随時接触角を再計算する
気液二相流のシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、気液二相流のシミュレーション方法に関する。
【背景技術】
【０００２】
従来から、混相流のシミュレーションを行う技術が知られている。
【０００３】
この種のシミュレーション技術としては、基板の状態毎に基板の領域を区分し、領域ごとに接触角θの関数である流体の移動速度式を設定することにより、液滴の挙動を時系列で把握し、計算時間の短縮化を図るものがある（特許文献１参照）。
【０００４】
また、インク放出のような流体の流れを解くために臨界角の大小で三重点の移動を判断し、表面張力を修正する技術も知られている（特許文献２参照）。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２０１０−０６７０２９号公報
【特許文献２】特開２００３−２８５４２８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
固体壁に付着した液滴が気体から外力を受ける状況では、液滴が気体から受ける力が臨界値を超えるまで、液滴は動き始めないことが実験的に確認されている。これは、外力によって液滴の形状が変化、すなわち、液滴の固体壁に対する接触角が変化することで、液滴に外力の向きと逆向きの抗力が働くためである。しかしながら、従来の手法では、液滴の固体壁に対する接触角に固定値を用いており、液滴に働く抗力が設定されないため、液滴が気体から受ける力が臨界値を超えるまで液滴が静止するという現象をシミュレーションにおいて再現できず、シミュレーション結果と実験結果とに乖離が生じる可能性があった。
【０００７】
そこで、本明細書開示の気液二相流のシミュレーション方法は、固体壁に付着した液滴が気体から外力を受け、液滴の固体壁に対する接触角が変化する場合の気液二相流の挙動のシミュレーションにおいて、シミュレーション精度を向上することを課題とする。
【課題を解決するための手段】
【０００８】
かかる課題を解決するために、本明細書に開示された気液二相流のシミュレーション方法は、固相に付着する液相の、前記固相に対する初期接触角を設定するステップと、前記初期接触角を用いて、前記固相に接触する前記液相と気相との界面における前記液相の速度である界面速度を算出する、第一の算出ステップと、前記固相に対する前記液相の随時接触角を、所定の条件を満たす、前記界面速度の関数として算出する、第二の算出ステップと、前記初期接触角と前記随時接触角との差が閾値以下の場合、計算が収束したと判定するステップと、を有し、前記所定の条件は、前記界面が前進しているとき、前記随時接触角は、前記液相が移動を開始するときの前記固相に対する前記液相の移動方向前側の接触角である第１臨界角以下の範囲で増加し、前記界面が後退しているとき、前記随時接触角は、前記液相が移動を開始するときの前記固相に対する前記液相の前記移動方向後側の接触角である第２臨界角以上の範囲で減少し、前記界面速度が０の場合、前記随時接触角は前記初期接触角から変化しない、との条件であり、前記判定するステップにおいて、前記初期接触角と前記随時接触角との差が閾値よりも大きい場合、前記設定するステップにおいて前記随時接触角を前記初期接触角に再設定し、前記第一の算出ステップにおいて前記界面速度を再計算し、前記第二の算出ステップにおいて前記随時接触角を再計算する気液二相流のシミュレーション方法である。
【発明の効果】
【０００９】
本明細書開示の気液二相流のシミュレーション方法によれば、固体壁に付着した液滴が気体から外力を受け、液滴の固体壁に対する接触角が変化する場合の気液二相流の挙動のシミュレーションにおいて、シミュレーション精度を向上することができる。
【図面の簡単な説明】
【００１０】
【図１】図１は、実施形態に係るシミュレーション方法を実行するシミュレーション装置のブロック図である。
【図２】図２（ａ）は、固体壁に付着した液滴が静止している状態を示す模式図であり、図２（ｂ）は、固体壁に付着した液滴が外力を受けている状態を示す模式図である。
【図３】図３は、実施形態に係るシミュレーション方法のフローチャートである。
【図４】図４は、固体壁に付着した液滴における界面の前進・後退について説明する図である。
【図５】図５は、前進接触角及び後退接触角について説明するための図である。
【図６】図６は、流路上に存在する液滴に圧力をかけた場合に液滴が動き出す圧力について、シミュレーション結果と実験結果とを比較した図である。
【図７】図７は、臨界圧力と溝流路の代表長さとの関係を示す図である。
【発明を実施するための形態】
【００１１】
以下、本発明の実施形態に係るシミュレーション方法について、添付図面を参照しつつ説明する。
【００１２】
まず、図１を参照し、シミュレーション装置の一例について説明する。図１は、本実施形態に係るシミュレーション方法を実行するシミュレーション装置１００のブロック図である。
【００１３】
図１に示すように、シミュレーション装置１００は、全体の動作を制御するＣＰＵ（Central Processing Unit）１１０、ＲＡＭ（Random Access Memory）１２０、ＲＯＭ（Read Only Memory）１３０、ＨＤＤ（Hard Disk Drive）１４０、出力インタフェース（Ｉ／Ｆ）１５０、及び入力インタフェース（Ｉ／Ｆ）１６０を備える。また、シミュレーション装置１００は、出力Ｉ／Ｆ１５０に接続された表示部１７０、入力Ｉ／Ｆ１６０に接続された入力デバイスとしてのキーボード１６２、及びポインティングデバイスとしてのマウス１６４を備える。ＨＤＤ１４０には、シミュレーションを実行するためのシミュレーションプログラムが格納されている。
【００１４】
シミュレーション装置１００は、例えば、燃料電池の単セルが有する流路における気液二相流の挙動をシミュレーションするためのシミュレーション装置であるが、これに限定されない。シミュレーション装置１００において、キーボード１６２及びマウス１６４を用いてシミュレーションプログラムを起動すると、ＣＰＵ１１０がシミュレーションプログラムをＲＡＭ１２０に読み出して実行することにより、流路における気液二相流の挙動がシミュレーションされる。本実施形態において、シミュレーション装置１００は、境界面の解析方法として公知の技術であるＶＯＦ（Volume of Fluid）法を採用する。
【００１５】
次に、図２（ａ）及び図２（ｂ）を参照して、液滴に働く、外力に反する抗力について説明する。図２（ａ）は、固体壁Ｗに付着した液滴Ｌが静止している状態を示す模式図であり、図２（ｂ）は、固体壁Ｗに付着した液滴Ｌが気体Ｇから外力Ｆを受けている状態を示す模式図である。なお、固体壁Ｗは固相の一例であり、液滴Ｌは液相の一例であり、気体Ｇは、気相の一例である。
【００１６】
図２（ａ）に示すように、固体壁Ｗに付着した液滴Ｌが静止している状態では、点Ａ及び点Ｂにおいて、接触角はいずれもθ０であり、表面張力の壁面平行成分は、いずれも、式(１)で表される。
【数１】
2017138123.tif 000003
なお、
【数２】
2017138123.tif 000004
は、気体Ｇと液滴Ｌとによって定まる定数であり、その向きは、固体壁Ｗ、液滴Ｌの物性値や、液滴Ｌが受ける外力等によって決定される。
【００１７】
液滴Ｌに働く、外力に反する抗力は、表面張力の壁面平行成分を、三相界面全周で積分することにより求められる。図２（ａ）に示すように、外力がなく、液滴Ｌが静止している場合には、液滴Ｌの全周で接触角は均一なため、抗力は０となる。
【００１８】
一方、図２（ｂ）に示すように、液滴Ｌが気体Ｇから外力Ｆを受けている状態では、液滴Ｌの形状が変形し、点Ａ´における接触角はθ１＞θ０となり、点Ｂ´における接触角はθ２＜θ０となる。このとき、点Ａ´における表面張力の壁面平行成分は式（２）で表され、点Ｂ´における表面張力の壁面平行成分は式（３）で表され、式（４）の関係にある。
【数３】
2017138123.tif 000005
外力Ｆの向きを図２（ｂ）の＋側と定義すると、外力Ｆに反する抗力Ｒは、式（５）で表されるように、負の値となり、図２（ｂ）に示すように、液滴Ｌの形状の変化に伴う接触角の変化によって、外力Ｆに反する抗力Ｒが表現される。
【数４】
2017138123.tif 000006
【００１９】
従来のシミュレーション手法では、接触角を図２（ａ）に示した静止状態の接触角から変化させないため、抗力Ｒを設定することができず、シミュレーションにおいて、小さい外力Ｆでも液滴Ｌが動き出してしまい、シミュレーション結果と実験結果とに乖離が生じる可能性があった。例えば、燃料電池セルの溝流路のように、断面積の小さい流路内の気液二相流においては、抗力が気液二相流の挙動に与える影響が大きいことが知られており、抗力が設定されていないことが、燃料電池セルの生成水の排水挙動のシミュレーション精度が低下する要因のひとつであると考えられる。
【００２０】
そこで、本実施形態に係るシミュレーション方法では、接触角を界面速度の関数として与えることで、外力Ｆによる接触角の変化を計算モデルに織り込む。これにより、接触角の変化による外力Ｆに反する抗力Ｒの変化が計算モデルに織り込まれることになり、シミュレーション精度が向上する。
【００２１】
以下、本実施形態に係るシミュレーション方法の具体的な計算処理について、図３及びその他図面を適宜参照して説明する。図３は、シミュレーション装置１００により実行される本実施形態に係るシミュレーション方法のフローチャートである。
【００２２】
図３の処理では、まず、要素毎に初期条件を設定する（ステップＳ１１）。初期条件として設定されるデータには、各要素に占める第一流体の体積率を表すＶＯＦ値ＮＶＯＦ、初期接触角、液滴、気体、及び固体壁の物性値等がある。
【００２３】
ＶＯＦ値ＮＶＯＦは、その要素における第一流体の体積率を示し、０〜１の値をとる。本実施形態では、ＮＶＯＦ＝０のとき、第二流体のみであることを示し、ＮＶＯＦ＝１のとき、第一流体のみであることを示すとする。なお、本実施形態において、第一流体を気体、第二流体を液滴とする。
【００２４】
また、本実施形態では、ステップＳ１１において、初期接触角に静止接触角（液滴が静止している状態での接触角）を設定する。なお、ステップＳ１１は、固相に付着する液相の、固相に対する初期接触角を設定するステップの一例である。
【００２５】
次に、ｔ＝０を設定する（ステップＳ１３）。
【００２６】
次に、三次元領域において、気液状態を読み込む（ステップＳ１５）。
【００２７】
次に、ステップＳ１１で設定された初期条件に基づいて、ＶＯＦ法による二相流計算を実行する（ステップＳ１７）。これにより、要素毎に、界面速度、圧力分布、及び流速分布などの流体解析結果が出力される。なお、ステップＳ１７は、初期接触角を用いて、固相に接触する液相と気相との界面における液相の速度である界面速度を算出する、第一の算出ステップの一例である。
【００２８】
次に、固体壁面に液滴及び気体の少なくとも一方が接する領域（流体領域）を読み込む（ステップＳ１９）。そして、ＶＯＦ値ＮＶＯＦに基づき、ステップＳ１９で読み込んだ流体領域から、気液固三相界面を抽出する（ステップＳ２１）。
【００２９】
次に、ステップＳ２１で抽出した気液固三相界面において、ＶＯＦ値ＮＶＯＦの勾配と界面速度ｖとに基づいて、界面が前進しているのか、後退しているのかを判定する（ステップＳ２３）。
【００３０】
ここで、図４に示すように、∇ＮＶＯＦは、気液界面に垂直で、液滴Ｌから気体Ｇへと向かう方向のベクトルとなる。本実施形態に係るシミュレーション方法では、気液固三相界面において、∇ＮＶＯＦと界面速度ｖとがなす角が９０度未満の場合、当該界面が前進していると定義し、９０度以上の場合、当該界面が後退していると定義し、界面の前進及び後退を、式（６）により判定する。
【数５】
2017138123.tif 000007
なお、界面速度ｖは、ステップＳ１７で計算された流速分布のうち、ステップＳ１９で読み込んだ流体領域の気体と液滴との界面における液滴の流速である。
【００３１】
次に、ステップＳ２３における界面の前進後退の判定結果と界面速度ｖとから、接触角の変化分を設定する（ステップＳ２５）。本実施形態に係るシミュレーション方法では、接触角の変化分Δθを界面速度ｖの関数として与える。Δθは、以下の条件１及び２を満たす界面速度ｖの関数であればよく、条件に応じて収束しやすい関数を探索するのが望ましい。
条件１：Δθは、界面が前進する場合に正の値となり、後退する場合に負の値となる。
条件２：∇ＮＶＯＦ・ｖ＝０のとき、Δθ＝０。
【００３２】
なお、収束性を向上させるために、関数は、以下の条件３を満たすことが望ましい。
条件３：｜ｖ｜が大きいほど、｜Δθ｜が大きい。
【００３３】
本実施形態に係るシミュレーション方法では、接触角の変化分Δθを、式（７）を用いて設定する。なお、式（７）において、ａは傾きを示す定数である。
【数６】
2017138123.tif 000008
【００３４】
続いて、ステップＳ２５で設定した接触角の変化分Δθと、接触角の最大値θａ及び最小値θｒとを用いて、随時接触角を算出する（ステップＳ２７）。随時接触角は、ステップＳ１７で計算された液滴の状態において実現される接触角である。なお、本ステップでは、接触角の最大値θａと最小値θｒとを、固体壁の物性として予め設定しておく。接触角の最大値θａ及び最小値θｒは、液滴が移動を開始する第１臨界角及び第２臨界角に相当し、例えば、図５に示すように、傾いた平面に置かれた液滴Ｌが移動を開始する際の前進接触角及び後退接触角としてそれぞれ測定できる。本実施形態に係るシミュレーション方法では、式（８）により、随時接触角を算出する。ただし、式（８）において、θ（ｔ）は現時刻における接触角（初期接触角）、θ（ｔ＋Δｔ）は次のタイムステップにおける接触角（随時接触角）を表す。
【数７】
2017138123.tif 000009
【００３５】
なお、ステップＳ２５及びＳ２７は、固相に対する液相の随時接触角を、界面が前進しているとき、随時接触角は、液相が移動を開始するときの固相に対する液相の移動方向前側の接触角である第１臨界角以下の範囲で増加し、界面が後退しているとき、随時接触角は、液相が移動を開始するときの固相に対する液相の移動方向後側の接触角である第２臨界角以上の範囲で減少し、界面速度が０の場合、随時接触角は初期接触角から変化しない、との条件を満たす界面速度の関数として算出する、第二の算出ステップの一例である。
【００３６】
次に、計算が収束したか否かを判定する（ステップＳ３１）。本実施形態では、初期接触角と随時接触角との差が、予め定めた閾値以下であった場合に、計算が収束したと判定する。ステップＳ３１は、初期接触角と随時接触角との差が閾値以下の場合、計算が収束したと判定するステップの一例である。
【００３７】
収束判定がＮＯの場合、値更新を行う（ステップＳ３３）。例えば、各要素の初期条件のデータを、ステップＳ１７の二相流計算により求められた値で更新するとともに、初期接触角をステップＳ２７で算出した随時接触角で更新する。ステップＳ３３の終了後は、ステップＳ３３において更新した値に基づき、ステップＳ１５からの処理が再度実行され、随時接触角が再計算される。なお、ステップＳ３３は、固相に付着する液相の、固相に対する初期接触角を設定するステップの一例である。
【００３８】
収束判定がＹＥＳの場合、シミュレーション終了条件（例えば、計算時間の上限）を満たしているか否かに基づいて、シミュレーションを終了するか否かを判定する（ステップＳ３５）。シミュレーション終了条件を満たしていない場合（ステップＳ３５／ＮＯ）、ｔ＝ｔ＋Δｔとし（ステップＳ３７）、次のタイムステップについて、ステップＳ１５からの処理を実行する。一方、シミュレーション終了条件を満たしている場合（ステップＳ３５／ＹＥＳ）には、図３の処理を終了する。
【００３９】
図３の処理によれば、随時接触角が界面速度ｖの関数として与えられるため、図２（ｂ）に示すように、固体壁Ｗに付着した液滴Ｌが気体Ｇから受ける外力Ｆによる接触角の変化が計算モデルに織り込まれる。これにより、接触角の変化による外力Ｆに反する抗力Ｒの変化が計算モデルに織り込まれることになり、固体壁に付着した液滴が気体から外力を受け、液滴の固体壁に対する接触角が変化する場合の気液二相流の挙動のシミュレーションにおいて、シミュレーション精度が向上する。
【００４０】
上述したシミュレーション方法及び従来手法により、流路上に存在する液滴に圧力をかけた場合に液滴が動き出す圧力（臨界圧力）をシミュレーションし、シミュレーション結果を実験結果と比較した。図６は、流路上に存在する液滴に圧力をかけた場合に液滴が動き出す圧力について、シミュレーション結果と実験結果とを比較した図である。
【００４１】
図６に示すように、従来手法では、低い臨界圧力で液滴が動き出してしまい、実験結果との乖離が大きい。一方、本実施形態に係るシミュレーション方法では、実験結果との乖離が小さい。すなわち、本実施形態に係るシミュレーション方法では、一定の圧力まで液滴が動き出さないという実験と同様の現象がシミュレーションにおいて再現されていると考えられ、シミュレーション精度が向上していることが確認できる。
【００４２】
次に、本実施形態に係るシミュレーション方法が有効な場合について検討する。臨界圧力の大きさは、溝流路の幅と深さとに依存する。例えば、正方形断面を仮定し、代表長さ（断面の一辺の長さ）に対して、臨界圧力をプロットすると、図７に実線で示す曲線となる。溝流路の排水性能を計算する場合、条件にも依存するが、臨界圧力がおよそ１００Ｐａ以上のとき、臨界圧力が排水挙動に及ぼす影響が大きくなる。当該事実と図６とにより、本実施形態にかかるシミュレーション方法は、代表長さ４ｍｍ以下の溝流路の排水性能を計算する際に有効であると考えられる。
【００４３】
なお、上記実施形態において、式（７）を、（∇ＮＶＯＦ・ｖ）の一次関数から二次関数、三次関数等にすることで、低流速領域における急激な接触角変化を抑えて、計算安定性を増すことができる。
【００４４】
上記実施形態は本発明を実施するための例にすぎず、本発明はこれらに限定されるものではなく、これらの実施例を種々変形することは本発明の範囲内であり、更に本発明の範囲内において、他の様々な実施例が可能であることは上記記載から自明である。
【符号の説明】
【００４５】
Ｌ：液滴（液相）
Ｇ：気体（気相）
Ｗ：固体壁（固相）
【図１】
2017138123.tif 000010
【図２】
2017138123.tif 000011
【図３】
2017138123.tif 000012
【図４】
2017138123.tif 000013
【図５】
2017138123.tif 000014
【図６】
2017138123.tif 000015
【図７】
2017138123.tif 000016