(57)【要約】【課題】より正確に自由界面を再構築できるシミュレーション装置等を実現する。【解決手段】空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する法線ベクトル算出部（３１）と、注目する計算セルの液相の質量変化率を算出する質量変化率算出部（３３）と、上記体積率と上記質量変化率算出部が算出した上記質量変化率とを用いて、注目する計算セルの次の時刻の液相の体積率を算出する体積率更新部（３４）とを備え、質量変化率算出部（３３）は、体積率が１または略１でない場合、注目する計算セルの体積率と法線ベクトル算出部（３１）が算出した法線ベクトルから算出した界面の面積に基づいて質量変化率を算出し、体積率が１または略１である場合、注目する計算セルに隣接する計算セルの質量変化率に基づいて注目する計算セルの上記質量変化率を算出する。【選択図】図１
2018036904.tif 000002
【特許請求の範囲】
【請求項１】
気相および液相を含む混相流における、物質移動を伴う界面を追跡するシミュレーション装置であって、
空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する法線ベクトル算出部と、
注目する計算セルの液相の質量変化率を算出する質量変化率算出部と、
上記体積率と上記質量変化率算出部が算出した上記質量変化率とを用いて、注目する計算セルの次の時刻の液相の体積率を算出する体積率更新部とを備え、
上記質量変化率算出部は、
上記体積率が１または略１でない場合、注目する計算セルの上記体積率と上記法線ベクトル算出部が算出した上記法線ベクトルから算出した界面の面積に基づいて、上記質量変化率を算出し、
上記体積率が１または略１である場合、注目する計算セルに隣接する計算セルの上記質量変化率に基づいて、注目する計算セルの上記質量変化率を算出する
ことを特徴とするシミュレーション装置。
【請求項２】
上記液相は、モノマーおよびポリマーの少なくとも一方を溶質として含み、
上記質量変化率算出部は、上記体積率と上記法線ベクトルから算出した界面の面積を、上記溶質の濃度に応じて修正した界面の面積に基づいて、上記質量変化率を算出する
ことを特徴とする請求項１に記載のシミュレーション装置。
【請求項３】
気相および液相を含む混相流における、物質移動を伴う界面を追跡するシミュレーション装置の制御方法であって、
空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する法線ベクトル算出ステップと、
注目する計算セルの液相の質量変化率を算出する質量変化率算出ステップと、
上記体積率と上記質量変化率算出ステップにて算出した上記質量変化率とを用いて、注目する計算セルの次の時刻の液相の体積率を算出する体積率更新ステップとを含み、
上記質量変化率算出ステップでは、
上記体積率が１または略１でない場合、注目する計算セルの上記体積率と上記法線ベクトル算出ステップにて算出した上記法線ベクトルから算出した界面の面積に基づいて、上記質量変化率を算出し、
上記体積率が１または略１である場合、注目する計算セルに隣接する計算セルの上記質量変化率に基づいて、注目する計算セルの上記質量変化率を算出する
ことを特徴とするシミュレーション装置の制御方法。
【請求項４】
請求項１に記載のシミュレーション装置としてコンピュータを機能させるための制御プログラムであって、上記法線ベクトル算出部、上記質量変化率算出部、および上記体積率更新部としてコンピュータを機能させるための制御プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、混相流における自由界面を体積追跡法に従ってシミュレーションするシミュレーション装置等に関する。
【背景技術】
【０００２】
従来から、混相流における自由界面の解析法について、種々の提案がなされている。自由界面の解析法の１つである体積追跡法は、体積占有率を用いて格子状に界面を表現する方法である。
【０００３】
例えば、特許文献１、２には、体積追跡法を用いて自由界面を再構築するシミュレーション装置であって、ＴＨＩＮＣ／ＷＬＩＣ（Tangent Hyperbola for Interface Capturing/Weighted Line Interface Calculation）法において値を固定されていた変数として、ｘ軸とｙ軸との両方の界面の法線ベクトルを用いた重み付けによって算出された数値を利用するシミュレーション装置が記載されている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１３−１８２４０３号公報（２０１３年９月１２日公開）
【特許文献２】特開２０１４−２６３７９号公報（２０１４年２月６日公開）
【非特許文献】
【０００５】
【非特許文献１】太田ら著、「混相流の数値シミュレーション」、丸善出版、2015年
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、複雑な自由界面の流れのシミュレーションは、まだ現実の流れを正確に再現できているとは言えず、検討の余地が残されている。特に、従来技術では、計算格子（計算セル）が液相で満たされており、かつ、当該計算格子内に界面が存在していない場合に、正確に自由界面の流れを再現できない可能性がある。
【０００７】
本発明は、前記の問題点に鑑みてなされたものであり、その目的は、より正確に自由界面を再構築できるシミュレーション装置等を提供することにある。
【課題を解決するための手段】
【０００８】
上記課題を解決するために、本発明に係るシミュレーション装置は、気相および液相を含む混相流における、物質移動を伴う界面を追跡するシミュレーション装置であって、空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する法線ベクトル算出部と、注目する計算セルの液相の質量変化率を算出する質量変化率算出部と、上記体積率と上記質量変化率算出部が算出した上記質量変化率とを用いて、注目する計算セルの次の時刻の液相の体積率を算出する体積率更新部とを備え、上記質量変化率算出部は、上記体積率が１または略１でない場合、注目する計算セルの上記体積率と上記法線ベクトル算出部が算出した上記法線ベクトルから算出した界面の面積に基づいて、上記質量変化率を算出し、上記体積率が１または略１である場合、注目する計算セルに隣接する計算セルの上記質量変化率に基づいて、注目する計算セルの上記質量変化率を算出することを特徴としている。
【０００９】
また、上記課題を解決するために、本発明に係るシミュレーション装置の制御方法は、気相および液相を含む混相流における、物質移動を伴う界面を追跡するシミュレーション装置の制御方法であって、空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する法線ベクトル算出ステップと、注目する計算セルの液相の質量変化率を算出する質量変化率算出ステップと、上記体積率と上記質量変化率算出ステップにて算出した上記質量変化率とを用いて、注目する計算セルの次の時刻の液相の体積率を算出する体積率更新ステップとを含み、上記質量変化率算出ステップでは、上記体積率が１または略１でない場合、注目する計算セルの上記体積率と上記法線ベクトル算出ステップにて算出した上記法線ベクトルから算出した界面の面積に基づいて、上記質量変化率を算出し、上記体積率が１または略１である場合、注目する計算セルに隣接する計算セルの上記質量変化率に基づいて、注目する計算セルの上記質量変化率を算出することを特徴としている。
【００１０】
上記の構成によれば、計算セルにおける体積率が１または略１でない場合、注目する計算セルの体積率と法線ベクトルから算出した界面の面積に基づいて質量変化率を算出し、体積率が１または略１である場合、注目する計算セルに隣接する計算セルの質量変化率に基づいて、質量変化率を算出する。
【００１１】
これにより、従来であれば、適切に質量変化率を算出できなかった、体積率が１または略１である場合であっても、質量変化率を適切に算出することができ、より正確に自由界面を再構築できるシミュレーション装置や、シミュレーション制御方法を提供することができる。
【００１２】
本発明に係るシミュレーション装置では、上記液相は、モノマーおよびポリマーの少なくとも一方を溶質として含み、上記質量変化率算出部は、上記体積率と上記法線ベクトルから算出した界面の面積を、上記溶質の濃度に応じて修正した界面の面積に基づいて、上記質量変化率を算出するものであってもよい。
【００１３】
上記の構成によれば、液相が二成分の場合であっても、適切に質量変化率を算出することができ、より正確に自由界面を再構築できるシミュレーション装置を提供することができる。
【００１４】
本発明の各態様に係るシミュレーション装置は、コンピュータによって実現してもよく、この場合には、コンピュータを上記シミュレーション装置が備える各部（ソフトウェア要素）として動作させることにより上記シミュレーション装置をコンピュータにて実現させるシミュレーション装置の制御プログラム、およびそれを記録したコンピュータ読み取り可能な記録媒体も、本発明の範疇に入る。
【発明の効果】
【００１５】
本発明の一態様によれば、従来であれば、適切に質量変化率を算出できなかった、体積率が１または略１である場合であっても、質量変化率を適切に算出することができ、より正確に自由界面を再構築できるシミュレーション装置等を提供することができるという効果を奏する。
【図面の簡単な説明】
【００１６】
【図１】シミュレーション装置の要部構成を示すブロック図である。
【図２】（ａ）は、体積率分布を説明するための前提となる状態を示す図であり、（ｂ）は、計算セルにおける体積率分布を示す図である。
【図３】ＳＬＩＣ、ＰＬＩＣ、およびＷＬＩＣにおけるアルゴリズムによって界面を追跡した場合を示す図である。
【図４】シミュレーション装置における液相体積率の更新処理の流れを示すフローチャートである。
【図５】計算セルと液相体積率αとの関係を説明するための図である。
【図６】体積率更新の原理を説明するための図であり、（ａ）は、αが１または略１でない場合を示す図であり、（ｂ）はαが１または略１である場合を示す図である。
【図７】計算セルにおける界面の位置が体積率の計算に及ぼす影響を説明するための図である。
【図８】計算セル中におけるポリマー（溶質）と溶媒との関係を模式的に示した図である。
【図９】シミュレーション装置における液相体積率の更新処理の流れを示すフローチャートである。
【図１０】単成分の場合で、当該成分の形状が四角形および円形のときの形状の時間的変化を示す図である。
【図１１】体積の時間的変化と理論値との関係を示す図である。
【図１２】液滴の初期形状の例を示す図である。
【図１３】図１１に示した３つの形状それぞれの液滴の体積の時間的変化と理論値との関係を示す図である。
【図１４】溶媒の蒸発に伴って気液界面の溶質濃度が上昇する様子を示す図である。
【発明を実施するための形態】
【００１７】
〔概要〕
〔気液界面における動き予測〕
一般に、界面における輸送に関しては、計算に用いる計算セル（計算格子）における、液相の体積率αを使用するＶＯＦ（Volume-of-Fluid）法が用いられる。図２を参照して説明する。図２は、計算セルにおける体積率分布を示す図である。なお、本実施形態では、説明の便宜上、適宜２次元の表記で説明するが、シミュレーションは３次元で行われるため、３次元に拡張して理解されたい。
【００１８】
図２に示すように、例えば、セルが完全に液相（液体）で占められている場合、α＝１であり、セルが完全に気相（気体）で占められている場合、α＝０である。αの値が０と１との間であれば、セル内に界面が存在していることを示している。そして、体積率は、輸送方程式に基づいて変化する。
【００１９】
この方法のポイントは、演算を発散させることなく、気液界面の輸送に関する、気液界面の正確な再現および表現ができることである。界面の長さの単位が小さい場合、表面張力は、再現および表現に大きな影響を与える力となる。
【００２０】
したがって、界面の再現、および表面張力の評価は、予測の正確さに影響を与える、鍵となるパラメータとなる。ＶＯＦ法は、十分な計算セル数を必要とするが、体積を適切に維持したまま、大きな密度の割合を取り扱うことができる。例えば、泡や滴の適切な予測のために必要な計算セルの数は、５よりも大きい。本実施形態に係るシミュレーション装置およびシミュレーション制御方法は、上記ＶＯＦ法を用いるものである。
【００２１】
また、界面の追跡アルゴリズムとしては、ＳＬＩＣ（Simple Line Interface Calculation）、ＰＬＩＣ（Piesewise Linear Interface Calculation）、およびＷＬＩＣ（Weighted Line Interface Calculation）がある。図３は、それぞれのアルゴリズムによって解明を追跡した場合を示す図である。図３では、それぞれ左側が実際の界面を示し、右側が追跡アルゴリズムによって再現された界面を示す。以下の実施形態では、ＷＬＩＣを用いた場合について説明するが、ＰＬＩＣやＳＬＩＣを用いることも可能である。
【００２２】
〔シミュレーション装置の構成〕
本実施形態に係るシミュレーション装置１は、気相および液相を含む混相流における、物質移動を伴う界面を追跡する機能を有する。
【００２３】
まず、図１を参照してシミュレーション装置１の要部構成を説明する。図１は、シミュレーション装置１の要部構成を示すブロック図である。図１に示すように、シミュレーション装置１には、入力装置２と表示装置３とが接続されている。
【００２４】
入力装置２は、シミュレーション装置１を操作するための装置である。入力装置２は、シミュレーション装置１を操作するための入力操作を受け付け、受け付けた入力操作の内容を示す操作データをシミュレーション装置１に送信できるものであればよく、例えばマウス、キーボード等を適用することもできる。
【００２５】
表示装置３は、シミュレーション装置１が出力する表示用データを用いて画像を表示する装置である。表示装置３は、画像を表示できるものであればよく、液晶表示装置、有機ＥＬ（Electro Luminescence）表示装置等を適用することもできる。表示装置３には、シミュレーション装置１を操作するための操作画面、シミュレーション結果等が表示される。
【００２６】
シミュレーション装置１は、制御部１０、および記憶部２０を備えている。
【００２７】
制御部１０は、シミュレーション装置１の動作を統括して制御するものである。制御部１０は、記憶部２０に格納されているプログラムを、例えばＲＡＭ（Random Access Memory）等に読み出して実行することによって上記の制御を行う。
【００２８】
記憶部２０は、シミュレーション装置１が動作するために必要なデータやプログラム等を格納している。また、シミュレーションの過程で算出される途中結果、シミュレーションの最終結果なども記憶部２０に格納される。
【００２９】
制御部１０には、シミュレーション算出部３０、シミュレーション制御部４０、および表示制御部５０が含まれている。シミュレーション算出部３０には、法線ベクトル算出部３１、体積率判定部３２、質量変化率算出部３３、および体積率更新部３４が含まれている。
【００３０】
法線ベクトル算出部３１は、注目する計算セルに含まれる界面の法線ベクトル（法線方向）を算出する。具体的には、空間に設定した計算セルごとに決定された、現在の時刻の液相の体積率のうち、注目する計算セルに隣接する計算セルの体積率を用いて、注目する計算セルに含まれる界面の法線ベクトルを算出する。法線ベクトルの算出方法の詳細については後述する。
【００３１】
体積率判定部３２は、注目する計算セルの液相の体積率が１、または略１（ほとんど１）か否かを判定する。
【００３２】
質量変化率算出部３３は、注目する計算セルの液相の質量変化率を算出する。質量変化率の算出方法の詳細については後述する。
【００３３】
体積率更新部３４は、質量変化率算出部３３が算出した質量変化率を用いて、計算セルにおける液相の体積率を更新し、注目する計算セルの次の時刻の液相の体積率を算出する。液相の体積率の更新方法の詳細については後述する。
【００３４】
シミュレーション制御部４０は、シミュレーション算出部３０を用いて、各時刻、各計算セルについての計算を行うことにより、流体のシミュレーションを実行し、実行結果を表示制御部５０に送信する。なお、シミュレーション制御部４０は、シミュレーションの実行結果を記憶部２０に格納してもよい。
【００３５】
より詳細には、シミュレーション制御部４０は、シミュレーション算出部３０による、ある計算セルにおける液相体積率の算出が終了すると、対象の計算セルを更新して、更新後の計算セルにおける液相体積率の算出をシミュレーション算出部３０に実行させる。これを、全ての計算セル（ｉ、ｊ、ｋの全ての組み合わせ）について計算が終了するまで、繰り返しシミュレーション算出部３０に実行させる。次に、シミュレーション制御部４０は、時刻を更新して、上述した、全ての計算セルにおける液相体積率の計算をシミュレーション算出部３０に実行させる。そして、シミュレーション制御部４０は、以上の処理を、シミュレーション終了時刻まで時刻を更新しながら、繰り返す。
【００３６】
表示制御部５０は、シミュレーション制御部４０が実行したシミュレーション結果を表示装置３に表示させる。また、表示制御部５０は、表示させるデータを記憶部２０から読み出し、表示装置３にこれを表示させる。すなわち、表示制御部５０は、上述したシミュレーション装置１によるシミュレーション結果を表示装置３に表示させたり、シミュレーション装置１を操作するための操作画面データ等を記憶部２０から読み出し、表示装置３に表示させる。
【００３７】
〔シミュレーション装置における液相の体積率の更新処理の流れ〕
次に、図４を参照して、シミュレーション装置１における液相の体積率αの更新処理の流れについて説明する。図４は、シミュレーション装置１における液相の体積率αの更新処理の流れを示すフローチャートである。
【００３８】
処理を開始すると、最初に、シミュレーション制御部４０が、計算の対象である注目する計算セル（対象セル）をシミュレーション空間の先頭の計算セルに設定する（Ｓ１０１）。
【００３９】
次に、図４に示すように、まず、法線ベクトル算出部３１が、現在の時刻の液相の体積率を用いて、注目する計算セルに含まれる界面の法線ベクトル（法線方向）を算出する（Ｓ１０２、法線ベクトル算出ステップ）。具体的に、図５を参照して説明する。図５は、計算セルと液相の体積率αとの関係を説明するための図である。図５に示すように、ｘ方向の計算セル位置をｉ、ｙ方向の計算セル位置をｊとしたときの計算セル（ｉ，ｊ）における液相の体積率をαｉ，ｊ（０≦α≦１）とする。そして、以下の式（１）により、ｘおよびｙについて、αの勾配をとることにより、界面の法線ベクトル（法線方向）ｎを算出する。
【００４０】
【数１】
2018036904.tif 000003
【００４１】
次に、体積率判定部３２は、注目する計算セルにおける液相の体積率αが１または略１であるか（１、または、ほとんど１であるか）を判定する。液相の体積率αが１、または、ほとんど１ではない場合（Ｓ１０３でＮＯ）、ステップＳ１０４へ進み、液相の体積率αが１、または、ほとんど１である場合（Ｓ１０３でＹＥＳ）、ステップＳ１０６へ進む。
【００４２】
ここで、液相の体積率αが１、または、ほとんど１である状態、および、１、または、ほとんど１ではない状態、換言すれば、αが１または略１である場合、およびαが１または略１ではない場合について、図６を参照して説明する。図６は、計算セルと液相の体積率との関係を説明するための図である。例えば、図６（ａ）に示す例のように、計算セルＣ１、およびＣ４では、対角線近傍に界面が存在しており、計算セルＣ１、およびＣ４におけるαは０．５前後となる。したがって、液相の体積率が１、または、ほとんど１ではない。すなわち、Ｓ１０３において「ＮＯ」に該当する。また、計算セルＣ２は、ほぼ気相で占められているのでα≒０となり、Ｓ１０３において「ＮＯ」に該当する。これに対して、計算セルＣ３はほぼ液相で占められているのでα≒１となり、Ｓ１０３において「ＹＥＳ」に該当する。
【００４３】
計算セル中に界面が明確に存在する場合（０＜α＜１）の場合に比べて、αが１、または、ほとんど１である場合（αが１または略１である場合）、すなわち計算セル中に界面が存在せず、計算セルが液相に占められている場合、従来の方法では、当該計算セル中に界面が存在しないので、界面の移動に伴う液相の質量変化率を適切に計算できない。具体的に、図７を用いて説明する。図７は、計算セルにおける界面の位置が質量変化率の計算に及ぼす影響を説明するための図である。図７（ａ）に示すように界面が存在した場合、計算セルＣ１およびＣ４では、計算セル中に界面が存在するので、従来の方法で、適切に質量変化率を計算することができる。
【００４４】
一方、計算セルＣ３では、本来であれば、図７（ｂ）に示すように、界面の位置が移動した場合（理想的には、法線方向にｌ（エル）だけ下がる場合）の液相の質量変化率を計算すべきであるが、計算セルＣ３に界面が存在しないため、図７（ｃ）に示すような形で液相の質量変化率が計算されてしまう。すなわち、計算セル中に界面が存在しないので、界面の移動も無く、したがって、液相の割合も変化しない。そこで、本実施形態では、計算セル中の液相の体積率αが１、または略１である場合、従来の方法と異なる計算方法を行うことによって、適切に液相の質量変化率を計算する。具体的には、図７（ｂ）と（ｃ）との差異である三角領域の面積（図７（ｄ））を計算することによって行う。
【００４５】
注目する計算セルの液相の体積率αが１、または、ほとんど１でない場合（Ｓ１０３でＮＯ）、ステップＳ１０４では、質量変化率算出部３３が従来の物質移動計算を行う（質量変化率算出ステップ）。すなわち、注目する計算セルの体積率と法線ベクトル算出部３１が算出した法線ベクトルから界面の面積を算出し、算出した界面の面積に基づいて、質量変化率を算出する。具体的には、以下の式（２）を用いて、液相の質量変化率Λを求める。
Λ＝ｋＳint（Ｃ＊−Ｃ∞）ＭＬ  （２）
ここで、ｋは物質移動係数、Ｃ＊は飽和濃度、Ｃ∞は気相のバルク濃度、ＭＬは液相の分子量である。またＳintは、界面の面積であり、例えばＰＬＩＣ法を用いて、界面を直線と推定することによって、液相の体積率αと法線ベクトルｎとから求めることができる。また、長さも計算することができる。なお、ＰＬＩＣ法は公知の技術であるので、その説明は割愛する。適宜、非特許文献１等を参照されたい。
【００４６】
一方、注目する計算セルの液相の体積率αが１、または、ほとんど１である場合（Ｓ１０３でＹＥＳ）、注目する計算セルに隣接する計算セル（周辺セル）における液相の体積率αが更新されているか否かを判定する（Ｓ１０６）。具体的には、後述するステップＳ１０７において、注目する計算セルの液相の質量変化率Λを算出するために必要となる注目する計算セルの上および右の計算セルにおける液相の体積率αが更新されているか否かを判定する。そして、更新されていれば（Ｓ１０６でＹＥＳ）、ステップＳ１０７へ進む。一方、更新されていなければ（Ｓ１０６でＮＯ）、ステップＳ１０５へ進む。
【００４７】
ステップＳ１０７では、注目する計算セルに隣接する計算セルの質量変化率に基づいて、注目する計算セルの質量変化率を算出する（質量変化率算出ステップ）。具体的には、以下の方法により、液相の体積率αが１、または、ほとんど１である計算セルにおける、液相の質量変化率Λを算出する。αｉ，ｊがほとんど１である場合（αｉ，ｊ≒１）、当該計算セルの上および右の計算セル（αｉ，ｊ＋１およびαｉ＋１，ｊ）は、０＜αｉ，ｊ＋１＜１、０＜αｉ＋１，ｊ＜１となる。
【００４８】
ここで、周囲の物質移動計算状況から，次式（３）により、界面の法線方向の移動量ｌの平均を算出する。
ｌ＝［average(Λ)／ρ］  （３）
ここで、ρは液相密度である。
【００４９】
本実施形態では、
average(Λ)＝（１／２）×（Λｉ，ｊ＋１＋Λｉ＋１，ｊ）    （４）
である。Λｉ，ｊ＋１及びΛｉ＋１，ｊは、上記式（２）で求められる。
【００５０】
そして、界面の勾配から、次式（５）を用いて界面とＸ軸との角度を求める。
φ＝ｔａｎ−１（ｎｙ／ｎｘ）    （５）
ここで、ｎｙはｙ方向の法線ベクトルであり、ｎｘはｘ方向の法線ベクトルであり、これらは、上記式（１）で求められる。
【００５１】
次に、上述した三角領域の面積θＴを次式（６）で求める。
θＴ＝ｌ２／ｓｉｎ２φ  （６）
最後に次式（７）を用いて、液相の質量変化率Λを求める。
Λ＝ρθＴ／Δｔ  （７）
ここで、ｔは時刻である。
【００５２】
つづいて、体積率更新部３４が、ステップＳ１０４またはステップＳ１０７で求めたΛを用いて、計算セルの液相の体積率αを更新する（Ｓ１０５、体積率更新ステップ）。具体的には、計算セルの格子体積をθ、液相密度をρとすると、
Ｄα／Ｄｔ＝−Λ／ρθ  （８）
であるため、オイラーの陽式スキームを適用して、次式（９）より、新しい時刻（次の時刻）の体積率αｎ＋１を算出する。
αｎ＋１＝αｎ＋Δｔ（−Λ／ρθ）  （９）
ここで、αｎは、既知（現在）の時刻における液相の体積率である。そして、ステップＳ１０８へ進む。
【００５３】
ステップＳ１０８では、全ての計算セルについて液相の体積率αが更新された否かを判定し、更新されていなければ（Ｓ１０８でＮＯ）、対象セルを更新して（Ｓ１０９）、ステップＳ１０２に戻る。一方、全ての計算セルについて液相の体積率αが更新されていれば（Ｓ１０８でＹＥＳ）、液相の体積率αの更新処理を終了する。
【００５４】
以上が、液相の体積率αの更新処理の流れである。
【００５５】
〔二成分液体の数的モデル化〕
上記では、単成分液体の数的モデル化について説明したが、以下では、二成分液体の数的モデル化について説明する。
【００５６】
乾燥過程における乾燥についての異なる処理条件と、膜の最終形状との関係を調査するためには、純粋な液体の代わりに多成分液体を考慮することが必要である。これまで説明してきた計算方法は、液相に１つの成分を含む単成分液体のためのものであるが、液相にポリマー等の溶質および溶媒を含む２成分液体に拡張することができる。具体的には、ポリマー濃度の輸送方程式（１０）を解くことによって実行できる。
【００５７】
【数２】
2018036904.tif 000004
【００５８】
ここで、Ｃは、液体中のポリマー濃度（mol／m3）であり、Ｄｐは、液体中のポリマーの拡散係数であり、Ｖは、液体中のポリマーの速度である。ポリマーは気体側へ移動しない、およびポリマー量は保存されるということを確実にするために、上式（１０）は、液体領域でのみ解かれる。
【００５９】
式（１０）の計算セルにおける液体領域において体積の積分を取り、計算セルで平均すると以下の式（１１）のようになる。
【００６０】
【数３】
2018036904.tif 000005
【００６１】
ここで、θＬは、液相の計算セルの格子体積である。
【００６２】
Ｓint，０／θ０（Ｓint，０は、ｔが０のときのＳintであり、θ０は、ｔが０のときの計算セルの格子体積である。）の右辺に発散定理を適用すると、上式（１１）は以下の式（１２）となる。
【００６３】
【数４】
2018036904.tif 000006
【００６４】
ここで、ηＬは、流体中における液相の体積割合であり、ＳＬは、界面の面積である。左辺におけるＣの対流項は、気液界面の輸送に関して上記式（２）または上記式（７）で計算されたΛを用いて容易に計算できる。右辺の界面での発散項は次式（１３）で計算できる。
【００６５】
【数５】
2018036904.tif 000007
【００６６】
ここで、Ｘは、界面での移動量であり、ＡＦは、面積割合である。
【００６７】
図８に示すように、界面におけるポリマーの存在は、溶媒の蒸発率に影響をあたえるものとして考慮される。図８は、計算セル中におけるポリマー（溶質）と溶媒との関係を模式的に示した図である。例えば、界面がポリマーによって完全に覆われている場合、溶媒の蒸発は起こらない。ポリマーの存在の効果は、考慮に入れなければならず、これは、界面でのポリマー濃度は、液体におけるポリマー濃度と同じと仮定することによって、実行される。
【００６８】
液相の質量変化率の計算において、上記式（２）のΛにおける界面の面積Ｓintは、次式で計算できるＳＬに置き換えられる。
ＳＬ＝Ｓint（１−Ｃ）  （１４）     
なお、上記式（１４）ではＳintに（１−Ｃ）を掛けたが、（１−Ｃ）の代わりに、液体中のポリマー濃度Ｃに応じて、式を適宜選択して用いることができる。
【００６９】
このように、液相に溶質と溶媒を含む２成分液体の場合は、溶質の濃度Ｃを算出し、質量変化率算出部３３にて、体積率と法線ベクトルから算出した界面の面積Ｓintを溶質の濃度に応じて修正した界面の面積ＳＬに基づいて、質量変化率を算出する。液相の体積率の更新処理の流れとしては、図９に示すとおり、図４のステップに、溶質の濃度を算出するステップ（Ｓ１１０）が加わる。
【００７０】
上記では、溶質としてポリマーを含む液相について述べたが、本実施形態は、ポリマーに限定されず、モノマーおよびポリマーの少なくとも一方を溶質として含むものであればよい。
【００７１】
〔シミュレーション結果〕
上述した方法を用いてシミュレーションした結果を、図１０〜図１４に示す。図１０は、単成分の場合で、当該成分の形状が四角形および円形のときの形状の時間的変化を示す図である。また、図１１は、体積の時間的変化と理論値との関係を示す図である。図１１に示すように四角形でも円形でも理論値とよく一致していることがわかる。
【００７２】
図１２は、液滴の初期形状の例を示す図である。ここでは、３つのφの値、すなわち６０°、９０°、および１２０°が示されている。図１３は、図１２に示した３つの形状それぞれの液滴の体積の時間的変化と理論値との関係を示す図である。図１３に示すように、何れの形状においても、理論値とよく一致していることがわかる。
【００７３】
最後に、図１４を参照して、蒸発により気液界面の溶質濃度が上昇する様子を説明する。図１４は、溶媒（水）の蒸発に伴って気液界面の溶質濃度が上昇する様子を示す図である。図１４では、実線で気液界面を表し、上から下にかけて経時的に変化する様子を示している。図１４に示すように、気液界面は時間の経過とともに溶質濃度が高くなっていくことが分かる。
【００７４】
〔まとめ〕
以上のように、本実施形態では、注目する計算セルの液相の体積率が１または略１である場合、隣接する計算セルの質量変化率に基づいて、注目する計算セルの質量変化率を算出する。これにより、液相の体積率が１または略１である場合でも、界面を高精度に追跡することができる。
【００７５】
なお、本実施形態は本発明の範囲を限定するものではなく、本発明の範囲内で種々の変更が可能であり、例えば、以下のように構成することができる。
【００７６】
本発明に係る計算方法は、溶質が溶媒に混ざっている溶液の計算方法であって、溶液の空間領域は物質移動を伴う境界面を含み、その境界面は溶媒の物質移動を含む情報を与えるものであり、この境界面を持った成分の生成、消散および空間移動現象を輸送方程式によって計算する方法であり、かつ計算する要素の形状に依存しない非物理的な誤差を回避する方法であってもよい。
【００７７】
さらに、本発明に係る計算方法は、境界面が空間的に移動を伴いながら可逆的あるいは非可逆的に変化する計算を含み、溶液の質量、運動量、エネルギー状態と相互干渉させながら計算する方法であってもよい。
【００７８】
さらに、本発明に係る計算方法は、溶液中の溶質がモノマーやポリマーを含む物質から構成されており、自由エネルギーの状態に応じて数ｎｍ〜数ｍｍにおよぶ微小な固体の結晶として析出する現象を計算する方法であってもよい。
【００７９】
〔ソフトウェアによる実現例〕
シミュレーション装置１としてコンピュータを機能させるための制御プログラム（特に制御部１０（シミュレーション算出部３０（法線ベクトル算出部３１、体積率判定部３２、質量変化率算出部３３、体積率更新部３４）、シミュレーション制御部４０、および表示制御部５０））は、集積回路（ＩＣチップ）等に形成された論理回路（ハードウェア）によって実現してもよいし、ＣＰＵ（Central Processing Unit）を用いてソフトウェアによって実現してもよい。
【００８０】
後者の場合、シミュレーション装置１は、各機能を実現するソフトウェアであるプログラムの命令を実行するＣＰＵ、上記プログラムおよび各種データがコンピュータ（またはＣＰＵ）で読み取り可能に記録されたＲＯＭ（Read Only Memory）または記憶装置（これらを「記録媒体」と称する）、上記プログラムを展開するＲＡＭ（Random Access Memory）などを備えている。そして、コンピュータ（またはＣＰＵ）が上記プログラムを上記記録媒体から読み取って実行することにより、本発明の目的が達成される。上記記録媒体としては、「一時的でない有形の媒体」、例えば、テープ、ディスク、カード、半導体メモリ、プログラマブルな論理回路などを用いることができる。また、上記プログラムは、該プログラムを伝送可能な任意の伝送媒体（通信ネットワークや放送波等）を介して上記コンピュータに供給されてもよい。なお、本発明は、上記プログラムが電子的な伝送によって具現化された、搬送波に埋め込まれたデータ信号の形態でも実現され得る。
【００８１】
本発明は上述した各実施形態に限定されるものではなく、請求項に示した範囲で種々の変更が可能であり、異なる実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても本発明の技術的範囲に含まれる。さらに、各実施形態にそれぞれ開示された技術的手段を組み合わせることにより、新しい技術的特徴を形成することができる。
【符号の説明】
【００８２】
１  シミュレーション装置
２  入力装置
３  表示装置
１０  制御部
２０  記憶部
３０  シミュレーション算出部
３１  法線ベクトル算出部
３２  体積率判定部
３３  質量変化率算出部
３４  体積率更新部
４０  シミュレーション制御部
５０  表示制御部
【図１】
2018036904.tif 000008
【図２】
2018036904.tif 000009
【図３】
2018036904.tif 000010
【図４】
2018036904.tif 000011
【図５】
2018036904.tif 000012
【図６】
2018036904.tif 000013
【図７】
2018036904.tif 000014
【図８】
2018036904.tif 000015
【図９】
2018036904.tif 000016
【図１０】
2018036904.tif 000017
【図１１】
2018036904.tif 000018
【図１２】
2018036904.tif 000019
【図１３】
2018036904.tif 000020
【図１４】
2018036904.tif 000021