JP2011022660A - 流体の数値計算装置 - Google Patents

流体の数値計算装置 Download PDF

Info

Publication number
JP2011022660A
JP2011022660A JP2009164887A JP2009164887A JP2011022660A JP 2011022660 A JP2011022660 A JP 2011022660A JP 2009164887 A JP2009164887 A JP 2009164887A JP 2009164887 A JP2009164887 A JP 2009164887A JP 2011022660 A JP2011022660 A JP 2011022660A
Authority
JP
Japan
Prior art keywords
fluid
equation
calculation
numerical
variation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
JP2009164887A
Other languages
English (en)
Inventor
Kota Nakano
浩太 中野
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Canon Inc
Original Assignee
Canon Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Canon Inc filed Critical Canon Inc
Priority to JP2009164887A priority Critical patent/JP2011022660A/ja
Publication of JP2011022660A publication Critical patent/JP2011022660A/ja
Withdrawn legal-status Critical Current

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

【課題】界面を含む圧縮性二相流の低速流れについて、各流体とその界面の運動を、精度良く、高速に、且つ、界面での数値振動を発生することなく、計算することのできる技術を提供する。
【解決手段】数値計算装置は、二流体の初期条件を含むデータを取得する取得手段と、前記データを用い、支配方程式を離散化された空間及び時間に関して解くことにより、前記二流体の運動を数値計算する計算手段と、前記計算手段の計算結果を出力する出力手段と、を備える。数値計算装置は、二流体の界面での速度の連続性を１流速仮定によって担保するとともに、質量保存則・エネルギー保存則・状態方程式によって数値解の探索範囲を限定することで、界面での数値振動の発生を抑制する。
【選択図】図２

Description

本発明は、二流体の運動をコンピュータにより解析するための技術に関する。
気液二相流は、産業上の様々な場所（例えば、原子力発電所、ボイラ、化学プラントなど）に現れ、工学的にも工業的にもきわめて重要な現象であると認識されている。このため、その挙動を予測することは産業上非常に大きな価値があり、その挙動予測の手段としての気液二相流の数値計算には様々な方法・装置が現在までに提案されている。気液界面の挙動予測が製品の機能改善に直結する産業分野も多い。例えばインクジェット技術における液滴挙動や、スクリュー回りのキャビテーション崩壊が代表的である。このような気液二相流に関わる現象を、その界面の運動とともに精度良く計算することは、その装置の性能改善に対して重要な示唆を与えることができる。
流体現象の数値的予測手法は、流体を非圧縮性流体とみなして計算する方法と圧縮性を考慮して計算する手法との二つに大別される。
非圧縮性流体の分野では、気液二相流に対する様々な計算方法が提案されている。しかしながら、当然であるが、非圧縮の仮定が妥当なケースにしか適用することはできないという制限がある。
圧縮性流体の数値解法には、時間積分の方法に大別して二つの方法がある。一つは密度ベース解法と呼ばれるものであり、もう一つは圧力ベース解法と呼ばれるものである。
密度ベース解法とは、圧縮性流体の支配方程式であるEuler方程式を陽的に時間発展さ
せ、求まった密度値から状態方程式を使って圧力値を求める手法である。密度ベース手法には様々な手法（非特許文献１にまとまった報告が記載されている）があるが、一般に、高速流に対して有効な手法であり、音速に比べて極端に遅い流れに適用すると、得られる数値解が非物理的な振動を起こしてしまうという難点がある。また、時間刻み幅が音速で規定されているために、解析対象の流体の速度が音速に比べて極めて遅い場合には、流動現象を追うために非現実的な計算時間を要求されてしまうという難点がある。
圧力ベース解法は、音波の伝播に相当する密度変化と状態方程式を結合させて、密度変化を圧力変化に還元して計算を行なう方法である。通常計算対象となる物質の状態方程式では、密度変化と圧力変化を比べた場合圧力変化の方がダイナミックレンジが広いため、圧力ベース解法は密度に関する数値誤差の影響が出にくいという利点がある。一方、圧力ベース解法では圧力を求めるために連立方程式を解く工程が現れるため、この連立方程式を解くために計算時間の多くを費やす必要がある。よって、１つの時間ステップの計算負荷としては、圧力ベース解法の方が密度ベース解法よりも計算負荷が大きい。しかしながら、時間刻み幅は密度ベース解法のように音速による条件がかからず時間刻み幅を大きく取ることができるため、音波の伝播現象に興味がなく、もっぱら流動現象が興味の対象である場合には、圧力ベース解法の方が総計算時間を小さくすることができる。
圧縮性の気液二相流に対して界面運動を直接計算する手法で今までに提案されているものには、Five-Equation model (非特許文献２)、Ghost fluid method（非特許文献３）、C-CUP法（非特許文献４）、二流体モデル（非特許文献５、６）等がある。
Five-Equation model, Ghost fluid methodは、多くの圧縮性流体解法と同様に密度ベ
ースの手法である。このため、低マッハ流れの計算では数値振動の発生というリスクを伴う。また、流動現象を対象とした時間スケールの計算を実用的な時間内に終えることは難
しい。
C-CUP法は圧縮性流体を圧力ベース解法にて解く解法であり、低速流に関しても数値振
動を起こすことなく計算をすることが可能である。しかしながら、C-CUP法は気液二相流
を一つの流体とみなしており、また、支配方程式として非保存形式を選択しているため、界面を通してやり取りされるエンタルピー収支をうまく計算することができない。界面における圧力仕事は界面運動の計算結果を左右するため、界面運動の様子を正確に表現できない。また、状態方程式を時間微分し、従属変数の１次微小変化に対する効果しか取り入れていないため、従属変数の値が時間的に大きく変化する場合は、得られた数値解が状態方程式から外れている可能性があるというリスクを伴う。C-CUP法における状態方程式か
らの逸脱という問題は、GCUP法（非特許文献８）にて解決がなされているが、１流体（ひとつの状態方程式）の場合である。
二流体モデルを支配方程式として半陰解法によって数値解を求める計算手法が非特許文献５、６に記載されている。二流体モデルは、それぞれの流体の質量保存則、運動方程式、エネルギー方程式を時間または空間積分した積分方程式を支配方程式としたものである。非特許文献６では、空間平均による支配方程式を用いており、また、圧力値は二つの流体で等しいという１圧力モデルを採用している。しかしながら、これらの文献の手法は、空間的に大きなスケール（例えば大きなタンクや管の中）を計算対象とし、内部に発生する細かな気液界面の挙動を直接の予測対象とせず、大きな空間スケールでみえる平均的な挙動をその予測対象としている。そのため、気液界面の効果を構成方程式として取り入れており、気液界面を数値計算格子で直接解像して計算するような計算手法ではないため、界面の運動を問題にするような対象を扱うことはできない。
二流体モデルを界面があるような問題に適用した計算を紹介している従来技術（非特許文献７）も存在するが、その場合は非圧縮性が仮定されており圧縮性が顕著に現れるような解析対象に対して界面運動を計算することはできない。
以上を総合すると、気液二相流についての従来技術を用いる場合、以下のような状況に対して精度の良い流体の挙動予測ができない。つまり、液体中にある蒸気泡の挙動をその気液界面の運動を含めて計算したい、といった場合である。このような状況の場合、蒸気泡に対して非圧縮性の条件を取ることはできない。もし非圧縮性条件を課した場合、蒸気の体積は一定であるから、その成長や収縮といった蒸気泡のダイナミクスは消えてしまうことになる。そのため、非圧縮性流体を仮定した様々な数値計算手法を使用することはできない。一方、従来の圧縮性流体の密度ベース数値計算手法を用いた場合、気液の質量密度比は数千倍にも及ぶため、界面付近で非物理的な数値振動が発生する可能性が極めて高い。さらに、気液界面の運動といった流動現象は音速に比べると非常に遅い現象であるから、その運動を計算するには非現実的な計算時間を要してしまう。また、圧縮性流体の圧力ベース解法であるC-CUP法を用いた場合、界面付近の圧力仕事や状態方程式に対する計
算誤差を許容しなければならないが、この誤差はここで取り上げているような例の場合は、無視しうる程度にならない場合が多い。
今取り上げた例は、気泡のキャビテーション現象についての本質であるから、その数値予測については産業上の重要性は非常に大きい。それにもかかわらず、従来技術では精度の良い計算が難しい。
Pieter Wesseling, Principles of Computational Fluid Dynamics, Springer Gregoire Allire, Sebastein Clerc, Samuel Kokh, A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids, J.Compt.Phiy. 181,577-616 (2002) Fedkiw, R. P., Aslam, T., Merriman, B., Osher, S, A non-oscillatory Eulerian approach to interfaces in multimaterial flows(the ghost fluid method). J. Comput. Phys. 152 (1999) 457-492 矢部孝、内海隆行、尾形陽一、『CIP法』、森北出版株式会社 日本原子力学会熱流動部会編『気液二相流の数値解析』、朝倉書店 秋山守、有冨正憲 監修、『新しい気液二相流数値解析−多次元流動解析−』、コロナ社 A.Minato, T.Nagayoshi, M.Misawa, A.Suzuki, H.Hinokata, S.Koshizuka, Numerical simulation method of complex 3D gas-liquid two-phase flow, Int.Conf.Multiphase Flows (ICMF ’04), Yokohama, Paper No.170, (2004) Toda S., Watanabe K., Utsumi T., Akamatsu M., and Fujii S., Two-phase flow computations using GCUP and level-set method, Proc. of the 2000 annual conference of Japan SIAM, lecture No.OS 1-3 (2000).
本発明の目的は、界面を含む圧縮性二相流の低速流れについて、各流体とその界面の運動を、精度良く、高速に、且つ、界面での数値振動を発生することなく、計算することのできる技術を提供することである。
本発明の第一態様は、二流体の初期条件を含むデータを取得する取得手段と、前記データを用い、支配方程式を離散化された空間及び時間に関して解くことにより、前記二流体の運動を数値計算する計算手段と、前記計算手段の計算結果を出力する出力手段と、を備え、ｋを流体成分を表す添え字、ｉ，ｊを空間成分を表す添え字とし、ｆｋを流体ｋの体積占有率、ρｋを流体ｋの質量密度、ｅｋを流体ｋの内部エネルギー、ｕを流体の速度、ｑｋを流体ｋの熱流速、Ｑｋを流体ｋに与えられる発熱量としたとき、前記支配方程式が下記（式１）〜（式５）を含むことを特徴とする流体の数値計算装置である。
Figure 2011022660
本発明の第二態様は、コンピュータが、二流体の初期条件を含むデータを取得する取得ステップと、コンピュータが、前記データを用い、支配方程式を離散化された空間及び時間に関して解くことにより、前記二流体の運動を数値計算する計算ステップと、コンピュータが、前記計算ステップの計算結果を出力する出力ステップと、を備え、ｋを流体成分を表す添え字、ｉ，ｊを空間成分を表す添え字とし、ｆｋを流体ｋの体積占有率、ρｋを流体ｋの質量密度、ｅｋを流体ｋの内部エネルギー、ｕを流体の速度、ｑｋを流体ｋの熱
流速、Ｑｋを流体ｋに与えられる発熱量としたとき、前記支配方程式が上記（式１）〜（式５）を含むことを特徴とするコンピュータによる流体の数値計算方法である。
本発明の第三態様は、上記数値計算方法の各ステップをコンピュータに実行させるためのプログラムである。
本発明によれば、界面を含む圧縮性二相流の低速流れについて、各流体とその界面の運動を、精度良く、高速に、且つ、界面での数値振動を発生することなく、計算することができる。
本発明の実施形態に係る流体の数値計算装置の構成例である。 本発明の実施形態に係る流体の数値計算方法のフローチャートである。 図２のステップＳ６の処理の詳細を示すフローチャートである。 計算結果の一例（気泡圧力の時間変化）を示したグラフである。
本発明は、二流体（二相流）の低速流れにおける各流体及びその界面の運動をコンピュータにより数値計算（シミュレーション）するための手法である。本発明は、気液二相流、固気二相流、固液二相流、液液二相流などいずれの二相流の数値解析にも適用可能であるが、圧縮性流体の界面運動の計算を精度良く行えることから、液体中の泡の挙動といった、気液二相流の数値解析に特に好ましく利用できる。
まず、本発明において採用した圧縮性二流体のモデルとその支配方程式の導出過程について説明する。その後、本発明の装置による具体的な数値計算手順を説明する。以下の説明では、ＭＫＳ単位系を用いる。空間座標の各成分ｘ，ｙ，ｚは、上付き又は下付きの空間添え字を用いて、それぞれ、ｘ１，ｘ２，ｘ３またはｘ１，ｘ２，ｘ３のように表す。時間はｔで表している。また、数式中の上下に同じ空間添え字が現れた場合は、その添え字について和をとるというアインシュタインの規約（縮約記法）を用いている。
（二流体モデルと支配方程式の導出）
２つの成分（流体の種類）からなる二流体を考える。異なる成分の間の界面が十分に薄く不連続面として捉えられるものと仮定する。その場合、ある１成分のみが空間の一点を占めることになるので、時間・空間に対する微分方程式として、各成分が従う方程式は単成分の方程式と変わらない。ただし、界面において、両成分の解を接続するための境界条件と、界面自身の運動方程式が必要となる。
流体の成分を下付き添え字ｋ（＝１，２）で表すと、各成分が従う方程式は以下の良く知られた流体の支配方程式（質量保存則、運動量保存則、エネルギー保存則）である。
Figure 2011022660
ここで、ρｋは流体ｋの質量密度[Kg/m3]、ｕｋは流体ｋの速度[m/sec]、ｐは流体の圧力[Pa]、Ｆｋは粘性などを含めた圧力勾配以外の体積力[N/m3]、ｅｋは流体ｋの内部エネルギー[J/Kg]、ｑｋは熱流速[W/m2]、Ｑｋは外部の発熱源によるエネルギー生成項（発熱量）[W/m3]である。空間成分を上付き又は下付きの添え字ｉ，ｊ，ｌ（＝１，２，３）で表している。圧力は、流体で区別無く、単一の圧力場で記述できるものとした。１圧力の仮定は、考えている流動の特徴的な時間スケールが流体ごとに異なる圧力の緩和時間より十分長い場合、有効である。
さらに、流体の熱力学的特性を表現する状態方程式 Ｆeos,kと、熱流速、内部エネルギーと温度の関係を与える法則が必要である。熱流速を与える法則としては、フーリエ則
Figure 2011022660
が一般的である。ここで、κは熱伝導率[W/mK]である。また、本説明では、内部エネルギーと温度Ｔ[K]の関係式として、
Figure 2011022660
を用いた。ここで、Ｃは流体の単位質量当たりの熱容量[J/KgK]である。
二相流を構成する各流体成分の空間分布を表すために以下のような関数を定義する。
Figure 2011022660
この関数ｆを体積占有率とよぶ。体積占有率ｆは次式を満たす。
Figure 2011022660
界面の運動方程式を以下で与える。
Figure 2011022660
ここで、ｕＩは界面の速度である。
さらに、界面において両流体の解を接続する境界条件が必要である。界面での境界条件は、界面における質量保存、運動保存、エネルギー保存であり、以下のようなものである。
Figure 2011022660
ここで、ｎＩ[1]は界面において成分ｋの内側から外側に向かう法線ベクトルである。ｍ
Ｉ[Kg/m2sec]は界面での質量生成、ＪＩ[Kg/m sec2]は界面での運動量生成、εＩ[J/m2sec]は界面でのエネルギー生成項である。これらは、相変化や界面張力（界面エネルギー）を考える場合に必要となる。
相変化が無視できる場合には、ｍＩ＝０より、
Figure 2011022660
である。さらに、
Figure 2011022660
なので、ｕＩが消去できて、
Figure 2011022660
が境界条件となる。これは、相変化がない場合、各成分の速度の界面方向成分が連続であることを要求している。
本発明では、体積平均による二流体モデルを使用するため、空間平均について整理しておく。体積平均を次式で定義する。
Figure 2011022660
ここで、Ｃは空間平均をとる領域を表す。
今後使用する式について、いくつか先に計算を済ませておく。
Figure 2011022660
時間微分は、空間平均の外に出すことができる。
Figure 2011022660
空間微分を平均の外に出す場合は、注意が必要である。簡単のために１次元で考えると、
Figure 2011022660
である。空間平均を考えているスケールＬが、計算対象としている領域のスケールに対して十分小さいとき、微分を外に出して考えても十分良い近似となる。本発明では、以後の支配方程式の説明では微分方程式の形で説明を行なうが、具体的な計算では積分方程式のままの計算を行なうので、空間微分についての微妙な議論を気にする必要は無い。
単独の流体の支配方程式に流体の空間分布を表す体積占有率ｆｋをかけて積分することで、空間平均した方程式を導く。
空間平均された質量保存則は以下のとおり。
Figure 2011022660
空間平均された運動量保存則は以下のとおり。
Figure 2011022660
空間平均されたエネルギー保存則は以下のとおり。
ただし、Ｅを全エネルギー密度[J/Kg]
Figure 2011022660
として記述している。
Figure 2011022660
エネルギー方程式は、後で内部エネルギー形式にするため、内部エネルギーが露わに見える形に書くと、以下のとおりである。
Figure 2011022660
全エネルギー保存則は、質量保存則と運動量保存則を考慮することで、内部エネルギーに対する保存則としても書くことができる。
運動量保存則（式２６）にｕｋを掛けることで運動エネルギーの時間変化を求めることができる。
Figure 2011022660
これに質量保存則（式２５）を代入すると、
Figure 2011022660
（式２５）、（式３１）を（式２９）に代入すれば、
Figure 2011022660
これが内部エネルギー表示のエネルギー保存則である。
空間平均した支配方程式は積分方程式である。解析したい領域を有限個の領域に区切り、その一つの領域をセルと呼ぶことにする。今、セルが両流体の界面を解像する程度に小さくとられた場合、多くのセルには単一の流体成分しか存在しない。また、相変化がない場合、界面における境界条件から、両流体の速度の界面に垂直な成分は連続であることが要求される。相変化がある場合も、その相変化による界面移動がセルの大きさに比べて十分に小さいと考えられる場合は、界面に垂直な速度成分の不連続性も十分小さいと考えてよい。これらの理由から、セルの大きさが十分に小さな場合、セルに定義される速度を両流体で連続にとる、すなわち、単一の速度場を用いることは良い近似となりうる。
今までに得た空間平均された支配方程式に、単一の速度場で表すという条件を入れると、界面における速度の不連続性はゼロであるから支配方程式を簡略化することができる。すなわち、
Figure 2011022660
とする。
この単一速度仮定によって、相変化がない場合の界面での境界条件：（式１４）、（式１５）、（式１６）を自然に満たすことができる。
さらに、運動量保存則（式２６）は、流体成分ｋについて和を取り、質量保存則とあわせることで、次のように変形できる。
Figure 2011022660
以上の議論により、体積平均された支配方程式は以下のように書ける。
Figure 2011022660
そもそも、流体の方程式自体が微小体積に関する積分方程式であるから、積分領域が十分に小さいと考えられるスケールで考えた場合は、特に積分記号を明記せずともよい。したがって、支配方程式は、
Figure 2011022660
であると考えてよい。
（式１）〜（式５）に加えて、物質の熱力学的性質を与える構成式を与える。本発明は、流体ごとに任意の状態方程式を与えることが可能である。ここでは具体的な説明のためにStiffened gas状態方程式を考え、以下の手順の説明を行なう。
Figure 2011022660
ここで、γｋ[1]は流体ｋの比熱比であり、Πｋ[N/m2]は物質ごとに違う定数である。ま
た、Ｃｋは流体ｋの定積熱容量であり、Ｔｋは流体ｋの温度、κｋは流体ｋの熱伝導度である。
（数値計算手順の説明）
続いて、上記方程式をコンピュータで数値計算するための具体的な計算手順を説明する
。ここでは簡単のため、１次元空間で説明を行なうが、計算手法の多次元空間への拡張は自明である。また、以下の説明では外力項を落としているが、必要ならば付け加えることができることも自明である。また、以下の説明では空間平均を表す記号＜＞を省略した。
支配方程式を時間方向について半陰的に離散化した方程式は次のようになる。
Figure 2011022660
ここで、上付き添え字ｎはｎ回目の時間ステップを表し、ｎ＋１はｎ＋１回目の時間ステップを表す。Δｔは、ｎステップ目とｎ＋１ステップ目の間の時間間隔である。
我々は方程式をスタッガード格子上で空間離散化を行なうものとする。下付き添え字の整数ｉにてスタッガード格子のセルの位置を表す。半整数はセルの表面位置を表す。また、Δｘは格子間隔である。
運動方程式（式３９）を空間方向に離散化して次式を得る。
Figure 2011022660
ここで、
Figure 2011022660
と置いた。
（式４２）で定義されるｕ＊の計算には、風上を考慮した移流スキームを用いるものとする。
質量保存則を空間方向に離散化して、それに（式４１）を代入することで次式を得る。
Figure 2011022660
ここで、
Figure 2011022660
とした。二重括弧 [[ ]] ではさまれている項については、風上を考慮したスキームで計算を行なうものとする。
エネルギー保存則を空間方向に離散化して、それに（式４１）を代入することで次式を得る。
Figure 2011022660
ここで、
Figure 2011022660
とした。
状態方程式を空間方向に離散化して次の式を得る。
Figure 2011022660
我々が解くべき式をまとめると、以下のとおりである。ここで、Ｍ１＝０、Ｍ２＝０は質量保存則であり、Ｅ１＝０、Ｅ２＝０はエネルギー保存則であり、Ｆ１＝０、Ｆ２＝０は状態方程式である。上付き添え字の１，２は、流体成分ｋを表す。
Figure 2011022660
各計算格子点について（式５１）の連立方程式を解く必要がある。未知数は、次時間ステップにおける圧力、体積占有率、各流体の温度、密度である。（式５１）はこれらの未知数について非線形であるから非線形連立方程式となる。
（式５１）の非線形連立方程式は反復解法により解く。本実施例では、以下に述べるように、非線形連立方程式の解法としてニュートン法を用いる。（式５１）を未知数の１次微小変分で局所線形化した方程式は、次のように書ける。
Figure 2011022660
ここで、ＪはJacobi行列であり次式で計算される。
Figure 2011022660
（式５２）の第５行目、６行目を用いて、１〜４行目のδρ1、δρ2を消去して得られる式
Figure 2011022660
を考える。
ここで、Ｊ１は以下のとおりである。
Figure 2011022660
（式５３）の左からＪ１の逆行列をかけることで次式を得る。実際の数値計算では、（式５２）を同値変形することで得られた（式５４）の方程式を解く。
Figure 2011022660
（式５４）の第１行目は、未知数として圧力変分δｐのみを含む一次方程式であるため、この式をすべてのセルで連立させた圧力方程式を数値的に解くことが可能である。この際には、例えばＢｉＣＧ（Bi-Conjugate Gradient）法などを用いて解けばよい。
圧力がすべてのセルにおいて求まれば、（式５４）の第２行以降から体積占有率変分δｆ、温度変分δＴが求まり、さらに（式５２）の第５、６行目から密度増分δρを求めることができる。
また、速度ｕは（式４１）で更新する。
上記プロセスを（式５１）の左辺Ｍ１、Ｍ２、Ｅ１、Ｅ２、Ｆ１、Ｆ２が十分にゼロに近づくまで繰り返し行なうことで、（式５１）の解、つまり、求める次時間ステップの解を得ることができる。
（装置構成）
次に、本発明に係る流体数値計算装置の一実施例について説明する。図１は、数値計算装置の構成例を示すブロック図である。この装置は、概略、入力部１１、計算部１２、出力部１３を備えるコンピュータにより構成される。入力部１１は、二流体の初期条件などのデータを取得するための取得手段である。計算部１２は、上述した計算処理を実行し、二流体の運動を数値計算するための計算手段である。出力部１３は、計算部１２の計算結
果を出力するための出力手段である。入力部１１はキーボード、マウス等で構成される。計算部１２はＣＰＵ１７およびメモリ１６で構成され、磁気ディスクやフラッシュメモリなどの記録媒体１５に記録された記録内容を読み込むことができる。出力部１３は、計算結果を表示するディスプレイや印刷による出力が可能な印刷部で構成されている。また、磁気ディスクやフラッシュメモリなどの計算結果保存部１４に計算結果を保存しておくことも可能である。記録媒体１５には、流体の数値計算プログラムが記録されている。計算部１２が、数値計算プログラムを読み込み、実行することにより、流体数値計算装置としての機能が実現される。
（数値計算処理のフロー）
図２及び図３のフローチャートに沿って、数値計算装置における数値計算処理のフローを説明する。以下の処理は、特に断りのない限り、計算部１２がプログラムに従って実行する処理である。
（１）ステップＳ１
まず、計算部１２は、数値計算に必要なデータを入力部１１から取得する。ここで入力されるデータには、圧力・体積占有率・温度・速度の初期条件、繰り返し計算の終了判定に用いる収束判定値、各流体の熱物性値（熱容量、熱伝導率、流体の状態方程式に現れる定数）などがある。なお、これらのデータをユーザに都度入力させることもできるし、予め記憶している設定値を読み込んでもよい。
（２）ステップＳ２
計算部１２が、流体内部の熱伝導方程式を流体成分ごとの温度場について計算することにより、熱拡散（熱伝導）による温度の予測値を求める。つまり、計算部１２は、各流体成分について、
Figure 2011022660
を解いて、各セル（格子点）での温度増分δＴ＃を求める。計算部１２によるステップＳ２の機能が本発明の第１計算部に対応する。
（３）ステップＳ３
計算部１２が、下記式により、各流体成分について温度増分δＴ＃による内部エネルギーの予測値Ｅ＃を求める。
Figure 2011022660
（４）ステップＳ４
計算部１２が、保存型の移流方程式により、速度の予測値ｕ＊を求める。
Figure 2011022660
この計算には、風上スキームを用いる。ただし、風上スキームとしてはどのような手法を用いても良い。
（５）ステップＳ５
計算部１２が、ステップＳ３、Ｓ４で得られた予測値Ｅ＃、ｕ＊を用いて、保存型の移流方程式により、体積占有率の予測値Ｄ＊、内部エネルギーの予測値Ｅ＊、質量の予測値Ｍ＊を求める。
Figure 2011022660
この計算には、風上スキームを用いる。ただし、風上スキームとしてはどのような手法を用いても良い。なお、計算部１２によるステップＳ４、Ｓ５の機能が本発明の第２計算部に対応する。
（６）ステップＳ６
次に、計算部１２が、ステップＳ４、Ｓ５で得られた予測値を用いて、連立方程式（式５１）（式４１）を解くことにより、次時間ステップにおける数値解を求める。計算部１２によるステップＳ６の機能が本発明の第３計算部に対応する。
図３は、ステップＳ６の処理の詳細を示している。計算部１２は、ステップＳ４、Ｓ５で計算された速度、質量、エネルギー、体積占有率それぞれの予測値（ｕ＊、Ｍ＊、Ｅ＊、Ｄ＊）を取得し（ステップＳ６０）、これらの予測値から開始して、次の繰り返し計算を実行する。
まず、計算部１２は、（式５４）の第１行目の式を全てのセル（空間格子点）について連立させた連立方程式をＢｉＣＧ法などを用いて解くことにより、全てのセルの圧力変分δｐの値を求める（ステップＳ６１）。次に、計算部１２は、（式５４）の第２〜４行目の式にδｐの値を代入することで、各セルの体積占有率変分δｆ、温度変分δＴ１、δＴ２を求める（ステップＳ６２）。さらに、計算部１２は、（式５２）の第５、６行目の式を用いて、各セルの密度変分δρ１、δρ２を求める（ステップＳ６３）。そして、計算部１２は、求めた各変分の値を、圧力ｐ、体積占有率ｆ、温度Ｔ１、Ｔ２、密度ρ１、ρ２の現在の値に加算することで、これらの未知数を更新する（ステップＳ６４）。また、計算部１２は、更新した圧力ｐの値を（式４１）に代入することで、速度ｕを更新する（ステップＳ６５）。なお、ステップＳ６１が本発明の第１の処理に、ステップＳ６２及びＳ６３が第２の処理に、ステップＳ６４が第３の処理に、ステップＳ６５が第４の処理に、それぞれ対応する。
計算部１２は、更新した未知数の値を（式４４）、（式４７）、（式５０）に代入することで、連立方程式（式５１）の左辺のＭ１、Ｍ２、Ｅ１、Ｅ２、Ｆ１、Ｆ２の値を計算する（ステップＳ６６）。そして、計算部１２は、Ｍ１、Ｍ２、Ｅ１、Ｅ２、Ｆ１、Ｆ２の値が所定の条件を満足するまで（ステップＳ６７）、ステップＳ６１〜Ｓ６６の処理を繰り返す。ステップＳ６７における収束判定はどのように行ってもよい。典型的には、左辺の値が、その方程式に含まれる定数項の値に対して十分に小さくなった場合に、十分にゼロに近づいた（収束した）とみなす。例えば、左辺の値をＶ、定数項の値をＣとしたときに、｜Ｖ｜／｜Ｃ｜＜1.0E-5を収束の条件とすることができる。この場合、1.0E-5が収束判定値である。
（７）ステップＳ７、Ｓ８
計算部１２は、時間ステップをインクリメントして（ステップＳ８）、ステップＳ２〜Ｓ６の処理を繰り返すことにより、所望の時間までの数値解を計算する。
（計算結果の例）
図４は、本実施例の数値計算装置によって計算した気泡成長における気泡中心の圧力履歴を示している（present method）。比較のため、非特許文献４に記載のＣＵＰ法に、非特許文献２で示されている気液混合のStiffened-gas状態方程式（Isobaric closure）を
使用して同じ問題を計算した場合の気泡中心の圧力履歴も示す（CUP+Five-Eqn. Model）
。また、Homogeneous modelとは、気泡内の物理量に空間分布がない１次元の断熱気泡モ
デルによる圧力の計算値である。
本発明による圧力変動の計算結果はHomogeneous modelの圧力履歴と良く一致している
が、従来手法を用いた計算結果はHomogeneous modelと大きく異なっている。これは、気
泡の成長収縮問題では、界面におけるエンタルピー収支を精度良く計算できることが正しいシミュレーション結果を与える条件であるから、本発明による計算が従来手法より優れていることを示すものである。
なお、上述した実施例は本発明の一具体例を示したものにすぎない。本発明は、上記実施例の構成に限られることはなく、その技術思想の範囲内で適宜変形が可能である。例えば、上記実施例では、スタッガード格子を用いているが、空間の離散化に他の手法を用いてもよい。また、上記実施例では、半陰解法を用いているが、陽解法や陰解法などの数値解法を用いてもよい。また、上記実施例では、ニュートン法を用いているが、他の反復解法を用いてもよい。
１１：入力部
１２：計算部
１３：出力部

Claims (8)

  1. 二流体の初期条件を含むデータを取得する取得手段と、
    前記データを用い、支配方程式を離散化された空間及び時間に関して解くことにより、前記二流体の運動を数値計算する計算手段と、
    前記計算手段の計算結果を出力する出力手段と、を備え、
    ｋを流体成分を表す添え字、ｉ、ｊを空間成分を表す添え字とし、
    ｆｋを流体ｋの体積占有率、ρｋを流体ｋの質量密度、ｅｋを流体ｋの内部エネルギー、ｕを流体の速度、ｑｋを流体ｋの熱流速、Ｑｋを流体ｋに与えられる発熱量としたとき、前記支配方程式が下記（式１）〜（式５）を含むことを特徴とする流体の数値計算装置。
    Figure 2011022660
  2. 前記計算手段は、
    流体内部の熱伝導方程式を流体成分ごとの温度場について計算することにより、温度の予測値を求める第１計算部と、
    前記第１計算部で得られた予測値を用いて、保存型の移流方程式を風上スキームによって計算することにより、速度、質量、エネルギー、及び、体積占有率の予測値を求める第２計算部と、
    前記第２計算部で得られた予測値を用いて、前記支配方程式及び状態方程式からなる連立方程式を解くことにより、次時間ステップにおける圧力、体積占有率、各流体成分の温度、各流体成分の密度、速度を求める第３計算部と、を含む
    ことを特徴とする請求項１に記載の流体の数値計算装置。
  3. 前記連立方程式は、圧力、体積占有率、各流体成分の温度、各流体成分の密度を未知数として含む非線形連立方程式であり、
    前記第３計算部は、前記第２計算部で得られた予測値から開始して、反復解法により前記連立方程式を解く
    ことを特徴とする請求項２に記載の流体の数値計算装置。
  4. 前記非線形連立方程式は、
    各流体成分の質量保存則：Ｍ１＝０、Ｍ２＝０、
    各流体成分のエネルギー保存則：Ｅ１＝０、Ｅ２＝０、
    各流体成分の状態方程式：Ｆ１＝０、Ｆ２＝０
    からなり（上付き添え字の１、２は、流体成分を表す。）、
    前記第３計算部は、
    前記非線形連立方程式を未知数の１次微小変分で線形化した連立方程式を同値変形することで得られる、圧力変分のみを未知数として含む方程式を数値的に解く第１の処理と、
    第１の処理で求めた圧力変分から、残りの未知数である体積占有率変分、温度変分、密
    度変分を求める第２の処理と、
    第１及び第２の処理で求めた各変分から、対応する未知数の値を更新する第３の処理と、
    第３の処理で求めた圧力の値を用いて、流体の速度の値を更新する第４の処理とを、
    第３及び第４の処理で更新された圧力、体積占有率、温度、密度、速度によるＭ１、Ｍ２、Ｅ１、Ｅ２、Ｆ１、Ｆ２の値が、所定の条件を満足するまで、繰り返し実行する
    ことを特徴とする請求項３に記載の流体の数値計算装置。
  5. 前記第３計算部は、
    状態方程式（Ｆ１＝０、Ｆ２＝０）を局所線形化した方程式を密度変分δρ１、δρ２に関して解いた式と、
    その密度変分δρ１、δρ２に関する式を用いて、エネルギー保存則（Ｅ１＝０、Ｅ２＝０）、及び状態方程式（Ｆ１＝０、Ｆ２＝０）を局所線形化した方程式からδρ１、δρ２を消去することによって得られた、圧力変分δｐ、体積占有率変分δｆ、温度変分δＴ１、δＴ２に関する式と、を用いるものであり、
    前記第１の処理では、前記圧力変分δｐに関する式を全ての空間格子点について連立させた連立方程式を解くことによって、全ての空間格子点のδｐの値を求め、
    前記第２の処理では、前記体積占有率変分δｆ及び温度変分δＴ１、δＴ２に関する式を用いて、前記第１の処理で求めたδｐの値からδｆ、δＴ１、δＴ２の値を求めた後、さらに、前記密度変分δρ１、δρ２に関する式を用いて、δｐ、δＴ１、δＴ２の値からδρ１、δρ２の値を求める
    ことを特徴とする請求項４に記載の流体の数値計算装置。
  6. 前記流体の状態方程式として、Stiffened gas状態方程式を用いる
    ことを特徴とする請求項２〜５のいずれかに記載の流体の数値計算装置。
  7. コンピュータが、二流体の初期条件を含むデータを取得する取得ステップと、
    コンピュータが、前記データを用い、支配方程式を離散化された空間及び時間に関して解くことにより、前記二流体の運動を数値計算する計算ステップと、
    コンピュータが、前記計算ステップの計算結果を出力する出力ステップと、を備え、
    ｋを流体成分を表す添え字、ｉ、ｊを空間成分を表す添え字とし、
    ｆｋを流体ｋの体積占有率、ρｋを流体ｋの質量密度、ｅｋを流体ｋの内部エネルギー、ｕを流体の速度、ｑｋを流体ｋの熱流速、Ｑｋを流体ｋに与えられる発熱量としたとき、前記支配方程式が下記（式１）〜（式５）を含むことを特徴とするコンピュータによる流体の数値計算方法。
    Figure 2011022660
  8. 請求項７に記載された流体の数値計算方法の各ステップをコンピュータに実行させるためのプログラム。
JP2009164887A 2009-07-13 2009-07-13 流体の数値計算装置 Withdrawn JP2011022660A (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2009164887A JP2011022660A (ja) 2009-07-13 2009-07-13 流体の数値計算装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2009164887A JP2011022660A (ja) 2009-07-13 2009-07-13 流体の数値計算装置

Publications (1)

Publication Number Publication Date
JP2011022660A true JP2011022660A (ja) 2011-02-03

Family

ID=43632703

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2009164887A Withdrawn JP2011022660A (ja) 2009-07-13 2009-07-13 流体の数値計算装置

Country Status (1)

Country Link
JP (1) JP2011022660A (ja)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102649159A (zh) * 2011-02-25 2012-08-29 北京科技大学 一种粉末注射成形坯密度在线预测系统及其方法
JP2012181248A (ja) * 2011-02-28 2012-09-20 Canon Inc 粉体の混合比計算方法及び装置
JP2012181247A (ja) * 2011-02-28 2012-09-20 Canon Inc 粉体の流動状態計算方法及び装置
KR101727022B1 (ko) * 2015-12-21 2017-04-14 한국과학기술정보연구원 전산유체역학에서 최적 2차 변수 계산 경로 제안 방법 및 장치
CN110309541A (zh) * 2019-05-31 2019-10-08 中国航天空气动力技术研究院 一种变比热气体的不同介质多组份流场界面条件构建方法
CN111783365A (zh) * 2020-06-04 2020-10-16 三多(杭州)科技有限公司 应用于低压界面处理的虚拟介质方法、装置和设备

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102649159A (zh) * 2011-02-25 2012-08-29 北京科技大学 一种粉末注射成形坯密度在线预测系统及其方法
JP2012181248A (ja) * 2011-02-28 2012-09-20 Canon Inc 粉体の混合比計算方法及び装置
JP2012181247A (ja) * 2011-02-28 2012-09-20 Canon Inc 粉体の流動状態計算方法及び装置
KR101727022B1 (ko) * 2015-12-21 2017-04-14 한국과학기술정보연구원 전산유체역학에서 최적 2차 변수 계산 경로 제안 방법 및 장치
CN110309541A (zh) * 2019-05-31 2019-10-08 中国航天空气动力技术研究院 一种变比热气体的不同介质多组份流场界面条件构建方法
CN110309541B (zh) * 2019-05-31 2023-04-07 中国航天空气动力技术研究院 一种变比热气体的不同介质多组份流场界面条件构建方法
CN111783365A (zh) * 2020-06-04 2020-10-16 三多(杭州)科技有限公司 应用于低压界面处理的虚拟介质方法、装置和设备
CN111783365B (zh) * 2020-06-04 2023-09-22 三多(杭州)科技有限公司 应用于低压界面处理的虚拟介质方法、装置和设备

Similar Documents

Publication Publication Date Title
Guo et al. An efficient lattice Boltzmann method for compressible aerodynamics on D3Q19 lattice
Yu et al. Direct numerical simulation of viscoelastic drag-reducing flow: a faithful finite difference method
Ceze et al. Anisotropic hp-adaptation framework for functional prediction
De Ridder et al. Modal characteristics of a flexible cylinder in turbulent axial flow from numerical simulations
Nichita et al. A level set method coupled with a volume of fluid method for modeling of gas-liquid interface in bubbly flow
Geng et al. Assessment of RANS turbulence models and Zwart cavitation model empirical coefficients for the simulation of unsteady cloud cavitation
Senocak et al. Interfacial dynamics‐based modelling of turbulent cavitating flows, Part‐2: Time‐dependent computations
Brouwer et al. Enriched piston theory for expedient aeroelastic loads prediction in the presence of shock impingements
Kollmannsberger et al. Fixed‐grid fluid–structure interaction in two dimensions based on a partitioned Lattice Boltzmann and p‐FEM approach
JP2011022660A (ja) 流体の数値計算装置
Margolin et al. Implicit turbulence modeling for high Reynolds number flows
Chen et al. Three-dimensional simplified and unconditionally stable lattice Boltzmann method for incompressible isothermal and thermal flows
Lu et al. Multi-parametric space-time computational vademecum for parametric studies: Application to real time welding simulations
Ceze et al. High-order output-based adaptive simulations of turbulent flow in two dimensions
CN113569450B (zh) 一种预估与控制液滴悬浮与驻留的方法
Yapici et al. Numerical analysis of viscoelastic fluids in steady pressure-driven channel flow
Sawadogo et al. Time domain simulation of the vibration of a steam generator tube subjected to fluidelastic forces induced by two-phase cross-flow
Takahira et al. An improved three-dimensional level set method for gas-liquid two-phase flows
Hosder Stochastic response surfaces based on non-intrusive polynomial chaos for uncertainty quantification
Jarosch Icetools: A full Stokes finite element model for glaciers
Wang et al. A new hybrid lattice-Boltzmann method for thermal flow simulations in low-Mach number approximation
Andersson et al. Numerical simulation of Stirling engines using an unsteady quasi-one-dimensional approach
Badwaik et al. Single-step arbitrary Lagrangian–Eulerian discontinuous Galerkin method for 1-D Euler equations
Codina A nodal‐based implementation of a stabilized finite element method for incompressible flow problems
Ryckelynck et al. Hyper-reduction framework for model calibration in plasticity-induced fatigue

Legal Events

Date Code Title Description
A300 Application deemed to be withdrawn because no request for examination was validly filed

Free format text: JAPANESE INTERMEDIATE CODE: A300

Effective date: 20121002