燃焼解析方法、燃焼解析装置、およびコンピュータプログラム
【課題】化学平衡法を用いつつも、非平衡な燃焼状態における計算精度の悪化を防止する。
【解決手段】燃焼を数値解析する燃焼解析方法であって、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、を含む。
【解決手段】燃焼を数値解析する燃焼解析方法であって、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、を含む。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、燃焼についての数値解析を行う方法、燃焼解析装置、およびコンピュータプログラムに関するものである。
【背景技術】
【０００２】
燃焼現象を、高速で簡潔に計算するために化学平衡法に基づくアルゴリズムが存在する（例えば、非特許文献１参照）。
燃焼においては、酸素と燃料が乱流や拡散によって混合され、これが燃焼反応の原料となる。化学平衡法のうち、例えば、ギブス自由エネルギー極小化法では、これらの原料中に含まれる元素（一般的な燃焼場では、Ｃ，Ｈ，Ｏ，Ｎの４元素）の数が、燃焼の前後で保存されるという条件と、系の持つギブスの自由エネルギーが最小になるという条件から導かれるｎ＋１個（ｎは原料中に含まれる元素数）多元連立非線形方程式を解くことで、燃焼生成物が算出される。
【０００３】
【非特許文献１】小川将希，大上芳文，離散渦法による燃焼解析，日本機械学会関西支部第82期定時総会講演会講演論文集，No.074-1 (2007), p.3-23
【発明の開示】
【発明が解決しようとする課題】
【０００４】
化学平衡法では、原料中に含まれる元素（例えば、Ｈ，Ｃ，Ｏ，Ｎ）の量が解れば、化学反応式を必要とすることなく、どんな生成物も予測できるというメリットがある。
しかし、化学平衡法に基づくアルゴリズムは、平衡でない現象に対しては計算精度の悪化が生じる。つまり、化学平衡法に基づくアルゴリズムでは、反応速度が無限大、すなわち反応時間が無限小、であることを前提としており、反応速度を制御できない。したがって、反応速度が遅い反応生成物では計算精度が低下する。
【０００５】
そこで、本発明は、化学平衡法を用いつつも、非平衡な燃焼状態における計算精度の悪化を防止することを目的とする。
【課題を解決するための手段】
【０００６】
本発明は、燃焼を数値解析する燃焼解析方法であって、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、を含むことを特徴とする燃焼解析方法である。
【０００７】
他の観点からみた本発明は、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する手段と、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する手段と、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する手段と、を備えていることを特徴とする燃焼解析装置である。
【０００８】
さらに他の観点からみた本発明は、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、をコンピュータに実行させるためのコンピュータプログラムである。
【０００９】
上記各発明によれば、燃焼に係わる化学種のうちの一部の化学種（第１化学種）については、化学平衡法に基づくアルゴリズムによって高速で簡潔に生成量を求めることができる。
一方、他の一部の化学種（第２化学種）については、反応速度論に基づくアルゴリズムによって反応速度を求めて、その反応速度から生成量を算出することで、反応速度を制御しつつ、生成量を算出できる。この結果、化学平衡法で求めると計算精度が悪化するような反応速度の遅い化学種を第２化学種とすることで、精度良く算出することができる。
【００１０】
したがって、例えば、燃焼に係わる化学種のうちの主要な化学種（第１化学種）については化学平衡法によって高速で簡潔に生成量を求めつつも、反応速度が遅く化学平衡法では計算精度が低下してしまうが計算精度が必要な化学種（第２化学種）については、部分的に反応速度論を用いて計算精度低下を防止することが可能である。
しかも、全ての化学種について計算量の多い反応速度論を用いるわけではないので、高速な演算が可能である。
【発明の効果】
【００１１】
よって、本発明によれば、化学平衡法を用いつつも、非平衡な燃焼状態における計算精度の悪化を防止することができる。
【発明を実施するための最良の形態】
【００１２】
以下、本発明の実施形態を図面に基づいて説明する。
【００１３】
［１．燃焼解析装置の全体構成］
図１は、実施形態に係る燃焼解析装置による燃焼解析処理の手順を示すフローチャートである。この燃焼解析装置は、図１のフローチャートの各ステップを実行するためのコンピュータプログラムが搭載されたコンピュータによって構成されている。
また、図２及び図３は、上記コンピュータプログラムによって実現される燃焼解析装置の機能ブロックを示している。
【００１４】
この燃焼解析装置は、解析対象の燃焼（乱流拡散燃焼）に係わる複数の化学種それぞれの生成量を算出することができる。燃焼解析装置は、解析対象の燃焼に係わる複数の化学種のうち、１又は複数の主要な化学種については、化学平衡法に基づくアルゴリズムに従って演算を行う。
一方、残りの一又は複数の化学種については、反応速度論（化学動力学）に基づくアルゴリズムに従って演算を行う。
【００１５】
このため、燃焼解析装置は、図２に示すように、解析対象の燃焼に係わる複数の化学種のうち、化学平衡法に従って演算が行われる１又は複数の化学種を「第１化学種」として設定するための第１化学種設定部１１と、反応速度論に従って演算が行われる１又は複数の化学種を「第２化学種」として設定するための第２化学種設定部２１と、を備えている。
したがって、燃焼解析装置のユーザは、設定部２１によって、化学平衡法ではなく、反応速度論で演算をしたい化学種（第２化学種）を自由に設定することができる。
【００１６】
また、燃焼解析装置は、設定（選択）された第１化学種それぞれについて、化学平衡法に基づくアルゴリズムに従って、第１化学種の質量分率（燃焼生成量）を算出する第１化学種の質量分率算出部１２を備えている。本実施形態において第１化学種の質量分率算出部１２は、化学平衡論（ＣＥ；Ｃｈｅｍｉｃａｌ Ｅｑｕｉｌｉｂｒｉｕｍ）の一種であるギブス自由エネルギー極小化法によって燃焼生成物を求める演算を行う。
【００１７】
ギブス自由エネルギー極小化法は、燃焼生成物を求める際に、平衡状態ではギブス自由エネルギーが極小値をとるという条件から元素の化学ポテンシャルを直接求めるものである。
化学平衡論で燃料生成物を求める方法としては、ギブス自由エネルギー極小化法以外に、圧力保存と原子数保存を拘束条件として平衡定数から各成分の濃度比を決める方法もあるが、この方法では、例えば、１００成分からなる燃焼場を解く場合、１００元の多元連立方程式を解く必要がある。
一方、ギブス自由エネルギー極小化法では、考慮する元素の数＋１の多元連立方程式を解けばよく、一般的な燃焼場では、Ｈ，Ｃ，Ｏ，Ｎの４元素を考慮すればよい。
【００１８】
さらに、燃焼解析装置は、設定（選択）された第２化学種それぞれについて、反応速度論に基づくアルゴリズムに従って、第２化学種の反応速度を算出する第２化学種それぞれの反応速度算出ｂ２２を備えている。本実施形態において第２化学種の反応速度算出部２２は、反応速度論に基づく燃焼モデルである渦消散コンセプト（ＥＤＣ；Ｅｄｄｙ Ｄｉｓｓｉｐａｔｉｏｎ Ｃｏｎｃｅｐｔ）モデルに従って反応速度を求める演算を行う。
なお、渦消散コンセプトモデル（ＥＤＣモデル）は、一段又は二段の総括反応の反応速度を算出することが可能な渦消散モデルを、多段階の化学反応式が解けるようにしたものである。
【００１９】
また、燃焼解析装置は、反応速度算出部２２によって求めた第２化学種それぞれの反応速度に基づいて、第２化学種それぞれの質量分率（燃焼生成量）を算出する第２化学種の質量分率算出部２３を備えている。
この第２化学種の質量分率算出部２３では、第２化学種の保存式に、前記反応速度を代入することで、燃焼生成物を求める。
【００２０】
また、渦消散コンセプトモデル（ＥＤＣモデル）は、ファイルスケールの体積分率算出部２４、タイムスケール算出部２５、原子の生成／消滅速度算出部２６も備えており、各値を算出することができる。
【００２１】
上記のように、本実施形態では、燃焼モデルが、化学平衡論（ＣＥ）と渦消散コンセプトモデル（ＥＤＣモデル）とによって構成されている。以下では、本実施形態の燃焼モデルを「ＥＤＣ／ＣＥモデル」とよぶ。
【００２２】
さらに、燃焼解析装置は、図３に示すように、燃焼解析のためのその他の演算を行う制御部３０を備えている。この制御部３０には、算出された第１化学種及び第２化学種の質量分率が与えられ、燃焼解析のための各種パラメータが算出される。制御部３０によって算出された各種パラメータは、第１化学種及び第２化学種の質量分率等の算出に用いられる。
この制御部３０は、処理ループを終了させるための収束判定部３１、速度と圧力を求める速度・圧力算出部３２、温度を算出する温度算出部３３、レイノルズ応力・乱流粘性係数を求めるレイノルズ応力・乱流粘性係数算出部３４、原子Ｉの質量分率を求める原子質量分率算出部３５、密度を求める密度算出部３６を備えている。
【００２３】
なお、図１のフローチャート及び図２及び図３ブロック図に関連して後に説明する数式は、上記コンピュータの記憶部に格納されたデータ（数式中の各項の値（パラメータ）に対応）に基づき、上記コンピュータの演算部によって演算され、その演算結果は前記記憶部に格納される。
【００２４】
［２．燃焼解析方法］
本実施形態に係る燃焼解析装置は、ＡＮＳＹＳ社の汎用流体解析ソフトウェアＦＬＵＥＮＴ（バージョン６．３．２６）に対して、本発明の機能をユーザ定義関数（User Defined Function）として組み込むことで実現した。ただし、本発明が適用されるソフトウェアがＦＬＵＥＮＴに限定されるものではない。
【００２５】
［２．１ 支配方程式］
解析対象の燃焼を支配する方程式としては、運動方程式、レイノルズ応力モデル、化学種Ｊの保存式、エネルギー方程式、気体の状態方程式、温度の算出式があり、具体的には下記の通りである。下記の支配方程式は、燃焼の数値解析の際に演算される。
【００２６】
［２．１．１ 運動方程式］
運動方程式は、下記の通りである。前記ＦＬＵＥＮＴでは、以下の連続の式とナビエストークスの式をSIMPLE 法［Suhas V Patanka. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, 1980 参照］により解く
【数１】
【数２】
【００２７】
［２．１．２ レイノルズ応力モデル］
本実施形態での解析対象の燃焼は、乱流であるため、乱流モデルとして、以下の式で表されるレイノルズ応力モデルを用いた。
【数３】
【００２８】
【数４】
【数５】
【数６】
【数７】
ここで、σkの推奨値は０．８２である。
【００２９】
θij,は、以下の式で与えられる。
【数８】
【００３０】
ここで，θij,1，θij,2 はそれぞれスローターム，ラピッドタームと呼ばれ，それぞれ以下のようにモデル化される。
【数９】
【数１０】
ここで、Ｃ1とＣ2はモデル定数であり、推奨値はそれぞれ１．８，０．６０であり、Ｐ＝（１／２）×Ｐkk，Ｃ＝（１／２）×Ｃkkである。
【００３１】
εijは、以下の式で与えられる。
【数１１】
【数１２】
【数１３】
ここで、Ｃμの推奨値は０．０９である。
なお、μは、粘性係数であり、流体が持つ固有の物性値である。一方、μtは、乱流によって生じる拡散の強さを表す指標（乱流による拡散率）であり、乱流モデルによって求められる。
【００３２】
［２．１．３ 化学種Ｊの保存式］
燃焼に係わる化学種Ｊの質量分率の式はＦＬＵＥＮＴにより以下の式で計算される。
【数１４】
【数１５】
【数１６】
【００３３】
［２．１．４ エネルギー方程式］
エンタルピーは、ＦＬＵＥＮＴのユーザ定義スカラー（ＵＤＳ）によって計算される。エネルギー方程式は下記の通りである。
【数１７】
【数１８】
【数１９】
【００３４】
［２．１．５ 気体の状態方程式］
気体の密度は、次式によって計算される。
【数２０】
【００３５】
［２．１．６ 温度の求め方］
温度は、次式で、ＦＬＵＥＮＴのＵＤＳにより計算される。
【数２１】
ここで、Ｔ0は参照温度で、２９８．１５Ｋとし、ＣpmJは化学種Ｊの平均定圧比熱、ｈJ0は化学種Ｊの標準生成エンタルピーを指す。
【００３６】
［２．２ ＥＤＣ／ＣＥモデルの説明］
本実施形態のＥＤＣ／ＣＥモデル（渦消散コンセプト・平衡法モデル）は、下記の通りである。なお、以下では、解析対象をＨ2−ａｉｒ乱流拡散火炎とし、Ｈ2をＥＤＣモデルにて計算した場合の例に基づいて説明する。
【００３７】
ａ．燃料（ここでは、Ｈ２とＮ２の混合気）の反応速度をＥＤＣデルで計算する。
ｂ．反応後の燃料は化学平衡に達していると仮定し、これを反応物として平衡計算をする。
【００３８】
前記反応速度は、ＥＤＣ（渦消散コンセプト）モデルの式で計算する。ＥＤＣモデルでは、燃焼は、乱流の微細構造、すなわちファインスケールで起こるものと仮定する。ファインスケールの体積分率は、次式で定義される（I R Gran and B F Magnussen. A numerical study of a bluff-body stabilized diffusion flame .part2. influence of combustion modeling and finite-rate chemistry. Combustion Science Technology, Vol.119, pp. 191_217, 1998. 参照）。
【数２２】
【００３９】
さらに、反応は次式で表されるタイムスケールで進行するものとする。
【数２３】
【００４０】
さらに、τ*のタイムスケールで反応が定圧反応容器内で起こるものと仮定し、格子内の現在の化学種と温度を初期条件とし、次の常微分方程式をｔ＝０〜τ*で数値積分する。なお、上記の「格子内」とは、数値解析の対象空間を格子状に分割した小さな空間をいう。数値流体力学では、その分割したそれぞれの空間（格子）に対して、運動方程式や例のずる応力などの方程式を解いていく。この格子は、セルとも呼ばれる。
【数２４】
【数２５】
【数２６】
ここで、Ｙ*J ，Ｒ*J はそれぞれファインスケール内の化学種Ｊの質量分率と反応速度を指し、ｖ'J K とν''J Kは、それぞれ化学種Ｊの化学反応式Ｋにおける正と負の量論係数、ＸJは化学種Ｊのモル濃度、ｋKは化学反応式Ｋの反応速度定数，ＫMAXは考慮する反応速度式の数である。Ａは頻度因子，Ｅは活性化エネルギー、βはそれぞれ温度の指数で実験パラメータである．また，式（２６）はアレニウスの式と呼ばれ，この式を含む常微分方程式は一般にstiff であることが知られている．この常微分方程式を数値積分するために，本実施形態ではＣＶＯＤＥ[S D Cohen and A C Hindmarsh. CVODE User Guide, llnl report edition, 1994. UCRL-MA-118618. 参照] を用いた。
【００４１】
実際の化学種Ｊの保存式（１４）に現れる反応速度は、次式でモデル化される。
【数２７】
【００４２】
また、燃料中のＮ2の反応速度は、Ｈ2に比例すると仮定し、燃料中のＮ2の質量分率及びＨ2の質量分率に基づき、次式で計算される。つまり、本実施形態における第２化学種の反応速度算出部２２は、燃料中のＮ2については、式（２８）で計算し、燃料中のＮ2以外の化学種については、式（２７）で反応速度を計算する。
【００４３】
【数２８】
【００４４】
上記の式（２７）は、詳細化学反応式を取り入れる場合である。しかし、化学反応式を取り入れず、反応は乱流混合が律速すると仮定すると、化学反応式は総括反応で表される。Ｈ2−ａｉｒを例にとると、Ｈ2＋１／２Ｏ2−＞Ｈ2Ｏとなり、このときの燃料であるＨ2の反応速度は以下の渦消散モデルで計算することができる。
【数２９】
ここで、Ｙfuel、Ｙox、Ｙpはそれぞれ燃料、酸化剤、燃焼生成物の質量分率を指し、ＡmixとＢmixはそれぞれ渦消散モデルの実験パラメータでそれぞれ４．０と０．５である。stは量論酸素燃料質量比を指す。
【００４５】
原子Ｉの生成または消滅速度は次式で計算される。
【数３０】
さらに、ａIJは化学種Ｊに対する原子Ｉの数の行列を示す．例えば，化学種Ｈ２に対する原子Ｈの数はａH,H2＝２であり、化学種ＯＨに対する原子Ｈの数はａH,OH＝１、原子Ｏの数は、ａO,OH＝１となる。
【００４６】
原子Ｉの保存式は次式で表される。
【数３１】
【００４７】
このときのエンタルピーは次式で近似される。
【数３２】
ＩMAXは、燃焼に関与する原子の数である。通常の燃焼場では、Ｈ，Ｃ，Ｏ，Ｎだけでよい。つまり、ＩMAX＝４となる。
【００４８】
本実施形態のＥＤＣ／ＣＥモデルでは、さらにＮを調節するパラメータを有する。
【数３３】
このＮＣＴＲＬの値を変化させることで、平衡計算による計算結果を調整できる。
平衡計算のためのプログラムは、ギブス自由エネルギー極小化法［JSME, editor. JSME Combustion Handbook. JSME, 2002. 参照］の理論に基づいたＮＡＳＡのＣＥＡ［S Gordon and B Mc Bride. Computer program for calculations of complex equilibrium compositions and applications. NASA Reference Publication, No. 1311, 1994. 参照］をベースに作成し、これをＵＤＦとしてＦＬＵＥＮＴに組み込んだ。
【００４９】
［２．３ 処理手順］
図１に戻り、燃焼解析の処理手順を説明する。ここでは、後述の図６の「Ｃａｓｅ３」の例に則して処理手順を説明する。なお、図１のフローチャートにおいて、ステップａ，ｂ，ｄ，ｅ，及びｇの部分が、ＥＤＣ／ＣＥモデルである。また、ステップａがＥＤＣモデルの部分、ｂがＣＥの部分、ｄ，ｅ，ｇがそれ以外で化学種の保存式や温度計算などを行っている部分である。その他のステップｃ，ｆ，ｈ，ｉは、ＦＬＵＥＮＴから提供される機能を利用した。
【００５０】
このＥＤＣ／ＣＥモデルにおいて考慮する化学種は、Ｈ2，Ｏ2，Ｈ2Ｏ，ＯＨ，Ｈ，Ｏ，ＨＯ2，Ｎ，Ｎ2である。そして、第１化学種設定部１１では、これらの化学種のうち、第１化学種として、Ｈ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2が設定（選択）されているものとする。また、第２化学種設定部２１では、残りの化学種であるＨ2，Ｈ，Ｏ，ＯＨが第２化学種として設定（選択）されているものとする。
【００５１】
ステップａ：
ステップａでは、ＥＤＣモデルの反応速度算出部２２によって、第２化学種であるＨ2，Ｈ，Ｏ，ＯＨ（Ｊ＝Ｈ2，Ｈ，Ｏ，ＯＨ）の反応速度ＲJを求める。この演算は、式（２４）〜（２９）に基づいて行われる。
また、ＥＤＣモデルの原子の生成／消滅速度算出部２６では、上記の反応速度から、式（３０）に基づき、原子Ｉの生成または消滅速度を求める。
また、ＥＤＣモデルでは、ファインスケールの体積分率算出部２４が、式（２２）に基づき体積分率を求め、タイムスケール算出部２５が、式（２３）に基づきタイムスケールを求める。
【００５２】
ステップｂ：
ステップｂでは、ＣＥ（ギブス自由エネルギー極小化法）の質量分率算出部１２によって、第１化学種であるＨ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2についての平衡計算を行う。この平衡計算では、原子Ｉ（Ｉ＝Ｃ，Ｈ，Ｏ，Ｎ）の質量分率ｄIを燃焼反応物とし、燃焼前後で式（３２）のエンタルピーは変化しないことを条件に、第１化学種であるＨ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2の質量分率（生成量）ＹH2O，ＹO2，ＹHO2，ＹN，ＹN2を求める。
【００５３】
ステップｃ：
ステップｃでは、速度・圧力算出部３２が、ＳＩＭＰＬＥ法（式（１）〜（２））に基づき、速度と圧力を求める。
【００５４】
ステップｄ：
ステップｄでは、温度算出部３３が、エネルギー方程式（式（１７））を解き、これから求まるエンタルピーから、式（２１）に基づき、温度を求める。
【００５５】
ステップｅ：
ステップｅでは、第２化学種の質量分率算出部２３が、ステップａで求めた第２化学種の反応速度ＲJそれぞれを、対応する第２化学種Ｊの保存式（式（１４））のＲJに代入し、第２化学種であるＨ2，Ｈ，Ｏ，ＯＨの質量分率（生成量）ＹH2，ＹH，ＹO，ＹOHを求める。
【００５６】
ステップｆ：
ステップｆでは、レイノルズ応力・乱流粘性係数算出部３４が、レイノルズ応力モデルに基づき、レイノルズ応力と乱流粘性係数を求める（式（３）〜（１３））。なお、レイノルズ応力は、運動方程式に現れる。また、乱流粘性係数は、エネルギー方程式、化学種Ｊの保存式、原子Ｉの保存式に現れる。
【００５７】
ステップｇ：
ステップｇでは、原子質量分率算出部３５が、ステップａで求めた原子Ｉ（Ｉ＝Ｃ，Ｈ，Ｏ，Ｎ）の生成／消滅速度ＲIそれぞれを、対応する原子Ｉの保存式（式（３１））に代入し、原子Ｉの質量分率ｄIを求める。
【００５８】
ステップｈ：
ステップｈでは、密度算出部３６が、化学種Ｊの質量分率ＹJと圧力から、式（２０）に基づき、密度を求める。
【００５９】
ステップｉ：
ステップｉでは、収束判定部３１が、図１のステップａ〜ｈの処理ループを終了させるために、収束判定を行う。
収束判定は、運動方程式、レイノルズ応力モデル、エネルギー方程式、化学種Ｊの保存式、原子Ｉの質量分率の式が正しい解に到達したかを判定する誤差判定である。
この収束判定は、ＦＬＵＥＮＴが本来有する機能によって自動的に実行される。ＦＬＵＥＮＴが採用している有限体積法をベースとしているＳＩＭＰＬＥ法では、離散化後、セルＰにおける任意の変数の保存式は、次式のように表される。
【数３４】
ここで、ａpは中心係数、ａnbは隣接セルの各影響係数、ｂは生成項Ｓ＝Ｓc＋ＳPφの式中の定数部分Ｓcを指す。式（３４）では次式が成り立つ。
【００６０】
【数３５】
【００６１】
さらに、このときのφの残差は次式で表される。
【数３６】
ＦＬＵＥＮＴでは、各変数φについて、この残差が一定の値以下（通例で、１．０×１０-3〜１．０×１０-6）になるまで、計算を続ける。
【００６２】
［３．実施例と比較例の計算条件］
本実施形態のＥＤＣ／ＣＥモデルによる計算精度を確かめるために、水素噴流火炎の計算（実施例と比較例）を行った。図４に計算領域と境界条件を示す。計算条件は高城・古藤の水素噴流火炎の実験[Toshimi Takagi and Satoru Kotou. A prediction of flow and combustion in turbulent diffusion flame in japanese. Transactions of the Japan Society of Mechanical Engineers, Series B, Vol. 48, No. 436, 1981. 参照] に準ずるものとした。
半径が２．４５ｍｍの円筒形ノズルから水素と窒素のモル比が０．４：０．６の燃料が５５．４ｍ／ｓで噴出されるものとし、周囲流は空気とした。空気については初めから平衡状態であるとし、モル比でｄＯ：ｄＮ＝０．２１：０．７９を入り口条件として与える。
【００６３】
さらに，燃料と周囲流の温度はともに３００Ｋとし、計算結果を比較する場所は、中心軸上とノズルから６０ｍｍと１６０ｍｍの地点の断面上とした。全ての方程式に二次精度の風上差分を適用し、図４の中心線を軸とした軸対象流れとし、乱流モデルとしてレイノルズ応力モデルを使用した。
【００６４】
図５に今回実行した計算（実施例と比較例）の代表的なシミュレーション条件（Ｃａｓｅ１，３，４）を示す。なお、Ｃａｓｅ２は欠番である。
また、図６に、上記計算で使用した化学反応式を示す。
このＥＤＣ／ＣＥモデルにおいて考慮する化学種は、Ｈ2，Ｏ2，Ｈ2Ｏ，ＯＨ，Ｈ，Ｏ，ＨＯ2，Ｎ，Ｎ2である。
【００６５】
図５のＣａｓｅ１（比較例）では、燃焼モデルとして従来のＥＤＣモデルを採用した。
このＣａｓｅ１のＥＤモデルでは、図６に示す化学反応式を用いた。
Ｃａｓｅ３（実施例）では、Ｈ2，Ｈ，Ｏ，ＯＨの化学種（第２化学種）のみ反応速度を考慮し、残りの化学種は第１化学種とし、平衡法で計算をした。
Ｃａｓｅ４（実施例）では、ＯＨが関わる図６の１〜８，１１，１３の化学反応式を考慮し、さらに、Ｈ2については、混合律速を仮定し、式（２９）で評価した。
【００６６】
［４．実施例及び比較例の計算結果と実験値との比較］
［４．０．１ 乱流エネルギーｋ、平均速度Ｕ、主要化学種の中心軸上での比較］
図７に、図４の中心線（中心軸）上における乱流エネルギーｋ、平均流速Ｕの分布を示し、図８に化学種と温度Ｔの分布を示す。
図７において、乱流エネルギーｋを見ると、実験値（ＥＸＰ）ではＬ＝１８０あたりまで緩慢に増えるのに対して、Ｃａｓｅ１，３では，Ｌ＝１００付近に向けて急激に増えて、ここを境に減少している。過去の計算例[Toshimi Takagi and Satoru Kotou. A prediction of flow and combustion in turbulent diffusion flame in japanese. Transactions of the Japan Society of Mechanical Engineers, Series B, Vol. 48, No. 436, 1981.][ Kimio Iino, Shinji Murakami, and Hoshitaka Kikkawa. Numerical simulation of high pressure methane-oxygen coaxial turbulent nonpremixed flames. Journal of High Temperature Society, Vol. 31, pp. 112_121, 2005.] においても、平均流速Ｕの実験値との整合性は良いが、ｋについてはシミュレーション値が過大評価されることが報告されている．この乱流拡散火炎のｋの過大評価を改善するためには，より高精度な乱流モデルのＬＥＳモデルを選択することや，レイノルズ応力モデルを改善する必要がある。
【００６７】
ここでは、全てのケースで中心軸上で温度の最大値の位置が一致するように図５のＲＳＭのパラメータを選んだ。実験ではノズルから１５０ｍｍ付近で温度は最大値の約１６００Ｋをとる。Ｃａｓｅ３のＨ2Ｏと温度の最大値は、それぞれ１８．９％、１５５６Ｋとなっており，実験値（ＥＸＰ）でＨ2Ｏと温度の最大値の約１８．３５％、約１６００Ｋとなっている。
【００６８】
［４．０．２ 乱流エネルギー、平均速度、主要化学種，温度のＬ＝６０での比較］
図９に、Ｌ＝６０での半径方向のｋ，Ｕ、図１０に、化学種と温度Ｔの分布を示す。図９のＣａｓｅ３、Ｃａｓｅ１のＵは、実験値（ＥＸＰ）でのＵと良く一致している。しかし、ｋについては中心軸上の分布と同じく実験値（ＥＸＰ）と差異が見られる。図１０ を見ると、Ｒ＝６ｍｍ付近の温度、Ｈ2ＯについてＣａｓｅ３、Ｃａｓｅ１ともにＥＸＰとの誤差が大きい。Ｒ＝６ｍｍ付近ではＴ，Ｈ2ＯのＥＸＰとの誤差は大きいがそれ以外の領域ではＥＸＰとの一致が良好なためである．
また、Ｃａｓｅ１，３ともにＲ＝６ｍｍ付近の温度を正しく予測できておらず，両者ともに似た予測値を示していることから，この付近のＴの誤差は燃焼モデルを起因とする誤差ではないと推測できる。
【００６９】
［４．０．３ 乱流エネルギー、平均速度、主要化学種、温度のＬ＝１６０での比較］
図１１にＬ＝１６０での半径方向のｋ，Ｕ、図１２に化学種と温度Ｔの分布を示す。
図１１を見ると、Ｒが大きくなるにつれて、Ｃａｓｅ１，３のＵの予測値と実験でのＵとの間の差異が大きくなる。この差異が存在する部分で拡散が大きく評価されれば、Ｕは拡散するので図１１中のＣａｓｅ１，３によるＵは小さくなると予測できる．このことから、この部分での差異は乱流拡散を実験よりも小さく見積もっていることが原因と推測される。さらにこの差異を改善するためには，より高精度な乱流モデルのＬＥＳモデルを選択することや、レイノルズ応力モデルを改善すればよい。また、図8１２を見ると、Ｃａｓｅ１，３の温度をＥＰＸと比較するとＲが大きくなるにつれて（火炎の外側へいくにつれて）Ｃａｓｅ１，３のＴの予測値と実験でのＴとの間の差異が大きくなる。これは、上記のＵ と同じく拡散を小さく見積もっていることと、さらにいえば、Ｕが大きいためＨ2Ｏが流されてしまったことが原因と考える。
【００７０】
［４．１ ラジカル類のEDC モデルとの比較］
［４．１．１ 中心軸上］
図１３，１４に中心軸上でのＯＨとＨのＥＤＣ／ＣＥモデルとＥＤＣモデルによる予測値を表す。図５にまとめたとおり、図５のＣａｓｅ１はＥＤＣモデル（比較例）、Ｃａｓｅ３，４はＥＤＣ／ＣＥモデル（実施例）での計算値を表す。Ｃａｓｅ３，４の違いはＥＤＣモデルで考慮する化学種の数の違いである。Ｃａｓｅ３ はＨ2，Ｈ，Ｏ，ＯＨの反応速度をＥＤＣモデルで考慮しているのに対し、Ｃａｓｅ４はＯＨのみを考慮している。
図１３に示すように、ＯＨの予測値を見ると、Ｃａｓｅ１，３，４との計算結果を比較するとＬ＝５０付近とＬ＝２００付近で差異が見られるが，それ以外の領域では整合性は良い．
これは、Ｃａｓｅ３，４ともにＯＨが式（２２）〜（２７）で示されるＥＤＣモデルの理論で計算されているからである。
【００７１】
図１４に示すように、Ｈの予測値ではＬ＝１５０まではＣａｓｅ１とＣａｓｅ３が良く一致している。しかし、Ｌ＝１５０以降はＣａｓｅ３のＨの予測値がＣａｓｅ１による予測値より小さく見積もるという差異が見られる。この差異は、Ｈ2Ｏ，Ｏ2，ＨＯ2の反応速度をＥＤＣモデルで考慮しなかったことが原因と考えられる。Ｃａｓｅ１とＣａｓｅ４ を比較すると、Case４ではＨをＣａｓｅ１より低く見積もっている。
【００７２】
［４．１．２ Ｌ＝６０］
図１６，１５にＬ６０でのＯＨとＨのＥＤＣ／ＣＥモデルとＣａｓｅ１のＥＤＣモデルによる予測値を表す。
図１６に示すＯＨの予測値を見ると、Ｃａｓｅ１とＣａｓｅ３，４との計算結果を比較するとＲ＝７付近まで良く一致しているが、それ以降ではＥＤＣ／ＣＥモデルとＥＤＣモデルの予測値に若干の差異が見られる。概ねＣａｓｅ１，３，４ともに近い予測値を示しているといえる。これは、Ｃａｓｅ３，４ともにＯＨについて反応速度を考慮したからである。
図１５に示すＨの予測値では、Ｃａｓｅ１とＣａｓｅ３を比較すると、Ｒ＝０付近とＲ＝６以降で若干の誤差が見られる。しかし、概ね両者の一致は良好である。
Ｃａｓｅ４のＨの予測値はＣａｓｅ１に比べ、Ｒ＝０付近では多く見積もり、Ｒ＝８付近を堺にしてＣａｓｅ１の予測値より低く見積もるようになる。
【００７３】
［４．１．３ Ｌ＝１６０］
図１８，１７にＬ＝１６０でのＯＨとＨのＥＤＣ／ＣＥモデルとＥＤＣモデルによる予測値を示す。
図１８のＯＨの予測値で、Ｃａｓｅ１，３，４との計算結果を比較すると、３者は似たグラフの傾きを示しているといえる。
図１７のＨの予測値では、Ｃａｓｅ１，３は良く一致している。Ｃａｓｅ４のＨの予測値はＣａｓｅ１に比べ、Ｒ＝０付近では多く見積もり、Ｒ＝１３付近を堺にしてＣａｓｅ１の予測値より低く見積もるようになる。ＨはＲ＝１７付近で完全に消費されている。
【００７４】
［５ 考察］
［５．１ 実施例の燃焼モデルと従来の燃焼モデルとの対比］
従来のＥＤＣモデルは乱流拡散火炎において詳細な化学反応式を考慮することができ、多成分系の反応を取り扱える。一方で、実施例に係るＥＤＣ／ＣＥモデルは化学反応式と化学平衡法の両方を用いて多成分系の反応を取り扱う。ここでは、ＥＤＣモデルと本実施例のモデルの計算値を比較し，両者の特徴について議論する。
【００７５】
ＥＤＣモデルでは，特定の化学種の発生量を調節するためには化学反応式の数を増やしたり、化学反応式中の実験パラメータとＥＤＣモデルのＣτ ，Ｃξ などのパラメータを調節する。しかし、これらのパラメータの組み合わせは多く、化学反応式の数を増やすとパラメータの数も増える。一方、実施例のモデルはＥＤＣモデルのパラメータと任意の化学反応式を考慮するだけで良いので、調節はＥＤＣデルより簡便である。
【００７６】
また、ＥＤＣ／ＣＥモデルでは考慮すべき化学反応式はシミュレーションの実行者の要求する計算精度に従う。
例えば、Ｃａｓｅのように任意の化学種（ここではＯＨ) について計算精度が必要であれば、その化学種について反応速度で考慮する。また、Ｃａｓｅ３のように他の化学種についても計算精度が必要ならば，他の化学種についても化学反応式を考慮する。こうすることで、ＥＤＣ／ＣＥモデルの計算結果はＥＤＣモデルの計算結果に近づけることができる。ＥＤＣ／ＣＥモデルでは、仮に平衡法で化学種を一つも考慮せず，全ての化学種をＥＤＣモデルで考慮したとすると、EDC モデルと同じになる。Ｃａｓｅ４で計算したラジカル類はEDC モデルの計算結果と比較すると、ＯＨについては良い一致を見せた。しかし、Ｈについては定量的な一致は見せなかった。
【００７７】
とはいえ、図１５，１７ を見ると、Ｈについても発生や消滅の傾向は捉えられていると考える。このことから、ＥＤＣ／ＣＥモデルの使い方としては，数値解析の初段階ではＣａｓｅ４のように必要な化学種のみにＥＤＣモデルによる反応速度計算を適用し、その後、他の化学種について精度が必要となれば、Ｃａｓｅ３のようにＥＤＣモデルで考慮する化学種を増やしていけばよい。ＥＤＣ／ＣＥモデルの適用としては，例えば、サーマルＮＯxを次の化学反応式で表される拡大ゼルドヴィッチ機構にて計算する際にＨ，Ｏ，ＯＨ については平衡法で考慮し、ＮＯxについてはＥＤＣモデルで考慮するといった使い方ができる。
【００７８】
【数３７】
【数３８】
【数３９】
【００７９】
他のＥＤＣ／ＣＥモデルの適用としては，乱流拡散火炎の詳細化学反応式を構築するためにも有用である。すなわち、Ｃａｓｅ４のように燃料の反応速度にだけＥＤＣモデルを適用し、他の化学種については平衡法を適用する。その後、燃焼場への寄与が大きいと推定される化学種について化学反応式を追加していく。燃焼場への寄与の大きさについては、平衡法で求まる予測値を参照すれば推測できる。
【００８０】
［５．１．１ 計算時間について］
ＥＤＣ／ＣＥモデルにおいて、平衡法では考慮する元素の数＋１を未知数とした非線形の連立方程式を解く。通常の燃焼場ではＣ，Ｈ，Ｏ，Ｎ＋１で５元となる。一方でＥＤＣモデルは考慮する化学種の数を未知数とした常微分連立方程式を解く。化学反応式が増え、考慮する化学種が増えるにつれ、常微分の連立方程式を解くための計算時間も増える。それに対して、平衡法では常に５元であり、反応物の濃度によって収束までの時間に差異はあるものの、EDC モデルより計算時間の増減は顕著ではない。また、ＥＤＣ／ＣＥモデルでは、平衡法で考慮する化学種については化学種Ｊの保存式（１４）を計算する必要がない。ＥＤＣモデルでは考慮する全ての化学種について化学種Ｊの保存式を解く必要がある。このために、ＥＤＣモデルでは考慮する化学種が増えるにつれ、解くべき化学種の保存式の数も増えるため計算時間は増加する。このような理由から、計算に考慮する化学種の数が多ければ多いほど、ＥＤＣ／ＣＥモデルはＥＤＣモデルより高速であるといえる。
【００８１】
［５．１．２ ＥＤＣ／ＣＥの適用範囲］
ＥＤＣ／ＣＥモデルの適用できる燃焼場は、例えば、乱流拡散火炎である。つまり、流れ場は乱流であり、燃料と酸化剤が別々の場所から流入する。また、反応速度の計算にＥＤＣモデルを使用しているため化学種の予測精度は最も良い精度でもＥＤＣモデルと同じとなる。ＮＯxなどの非平衡性の強い化学種に対して平衡法を適用すると予測制度が低下するので、このような化学種についてはEDC モデルで反応速度を計算しなければならない。
【００８２】
また、最終的に生じる生成物、例えばＨ2ＯやＣＯ2については平衡法で考慮する。燃焼によって主に生成する化学種を定めることで、平衡法の計算を安定化させるためである。言い換えれば、燃焼生成物の化学種の大部分を平衡法で計算し，精度が必要な化学種に限って補助的にＥＤＣモデルを取り入れることが好ましい。また、ＥＤＣ／ＣＥモデルではアレニウスの式を使っており、アレニウスの式は温度に依存するのでしばしば消炎が起こる。この問題については，燃料の反応速度についてアレニウスの式を考慮しない、式（２９）を使うと回避できる。
【００８３】
［６ 結論」
結論として、実施例に係るＥＤＣ／ＣＥモデルの特徴と、その比較・検討についてまとめる。
（１）詳細化学反応式を必要とせず全ての化学種の生成量を予測でき，任意の化学種について化学反応式を取り入れることでＥＤＣモデルと同程度の精度で計算することができる。
（２）Ｈ2Ｏなどの最終生成物だけでなく、ＯＨ，Ｏ，Ｈなどの中間生成物も予測することができ，これとゼルドビィッチ機構などの既存のモデルと組み合わせることでＮＯxなどの生成量の予測などの応用が可能である。
（３）乱流混合・燃焼をＥＤＣモデルにより評価しているので，化学反応式を考慮し化学反応に乱流と温度の影響を取り入れることができる。
【００８４】
さらに、実験値（ＥＸＰ）、従来のＥＤＣモデルによる予測値との比較・検討についてまとめる。
（１）Ｈ2，Ｏ2，Ｈ2Ｏについて、ＥＤＣ／ＣＥモデルによる予測値は、実験値とＥＤＣ モデルによる予測値と一致する。
（２）ラジカル類である、ＯＨ，Ｈについては反応速度を考慮することで、ＥＤＣモデルの予測値と良い一致を見せるようになる。
（３）多くの化学種について反応速度を考慮すると、ＥＤＣモデルによる予測値に近づく。
【００８５】
なお、本発明は、上記実施形態及び実施例に限定されるものではなく、様々な変形が可能である。
【図面の簡単な説明】
【００８６】
【図１】燃焼解析処理のフローチャートである。
【図２】燃焼解析装置の機能ブロック図である。
【図３】燃焼解析装置の機能ブロック図である。
【図４】水素噴流火炎の計算領域と境界条件を示す図である。
【図５】Ｈ２−Ｏ２化学反応メカニズムと関連するレート係数を示す表である。
【図６】シミュレーション条件を示す表である。
【図７】中心線に沿った位置での乱流エネルギーと平均速度を示すグラフである。
【図８】中心線に沿った位置での化学種のモル分率と温度を示すグラフである。
【図９】Ｌ＝６０の位置での乱流エネルギーと平均速度を示すグラフである。
【図１０】Ｌ＝６０の位置での化学種のモル分率と温度を示すグラフである。
【図１１】Ｌ＝１６０の位置での乱流エネルギーと平均速度を示すグラフである。
【図１２】Ｌ＝１６０の位置での化学種のモル分率と温度を示すグラフである。
【図１３】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、中心線に沿った位置でのＯＨのモル分率を示すグラフである。
【図１４】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、中心線に沿った位置でのＨのモル分率を示すグラフである。
【図１５】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝６０の位置でのＨのモル分率を示すグラフである。
【図１６】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝６０の位置でのＯＨのモル分率を示すグラフである。
【図１７】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝１６０の位置でのＨのモル分率を示すグラフである。
【図１８】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝１６０の位置でのＯＨのモル分率を示すグラフである。
【符号の説明】
【００８７】
１１ 第１化学種設定部
１２ 第１化学種の質量分率算出部
２１ 第２化学種設定部
２２ 第２化学種の反応速度算出部
２３ 第２化学種の質量分率算出部
【技術分野】
【０００１】
本発明は、燃焼についての数値解析を行う方法、燃焼解析装置、およびコンピュータプログラムに関するものである。
【背景技術】
【０００２】
燃焼現象を、高速で簡潔に計算するために化学平衡法に基づくアルゴリズムが存在する（例えば、非特許文献１参照）。
燃焼においては、酸素と燃料が乱流や拡散によって混合され、これが燃焼反応の原料となる。化学平衡法のうち、例えば、ギブス自由エネルギー極小化法では、これらの原料中に含まれる元素（一般的な燃焼場では、Ｃ，Ｈ，Ｏ，Ｎの４元素）の数が、燃焼の前後で保存されるという条件と、系の持つギブスの自由エネルギーが最小になるという条件から導かれるｎ＋１個（ｎは原料中に含まれる元素数）多元連立非線形方程式を解くことで、燃焼生成物が算出される。
【０００３】
【非特許文献１】小川将希，大上芳文，離散渦法による燃焼解析，日本機械学会関西支部第82期定時総会講演会講演論文集，No.074-1 (2007), p.3-23
【発明の開示】
【発明が解決しようとする課題】
【０００４】
化学平衡法では、原料中に含まれる元素（例えば、Ｈ，Ｃ，Ｏ，Ｎ）の量が解れば、化学反応式を必要とすることなく、どんな生成物も予測できるというメリットがある。
しかし、化学平衡法に基づくアルゴリズムは、平衡でない現象に対しては計算精度の悪化が生じる。つまり、化学平衡法に基づくアルゴリズムでは、反応速度が無限大、すなわち反応時間が無限小、であることを前提としており、反応速度を制御できない。したがって、反応速度が遅い反応生成物では計算精度が低下する。
【０００５】
そこで、本発明は、化学平衡法を用いつつも、非平衡な燃焼状態における計算精度の悪化を防止することを目的とする。
【課題を解決するための手段】
【０００６】
本発明は、燃焼を数値解析する燃焼解析方法であって、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、を含むことを特徴とする燃焼解析方法である。
【０００７】
他の観点からみた本発明は、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する手段と、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する手段と、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する手段と、を備えていることを特徴とする燃焼解析装置である。
【０００８】
さらに他の観点からみた本発明は、解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、をコンピュータに実行させるためのコンピュータプログラムである。
【０００９】
上記各発明によれば、燃焼に係わる化学種のうちの一部の化学種（第１化学種）については、化学平衡法に基づくアルゴリズムによって高速で簡潔に生成量を求めることができる。
一方、他の一部の化学種（第２化学種）については、反応速度論に基づくアルゴリズムによって反応速度を求めて、その反応速度から生成量を算出することで、反応速度を制御しつつ、生成量を算出できる。この結果、化学平衡法で求めると計算精度が悪化するような反応速度の遅い化学種を第２化学種とすることで、精度良く算出することができる。
【００１０】
したがって、例えば、燃焼に係わる化学種のうちの主要な化学種（第１化学種）については化学平衡法によって高速で簡潔に生成量を求めつつも、反応速度が遅く化学平衡法では計算精度が低下してしまうが計算精度が必要な化学種（第２化学種）については、部分的に反応速度論を用いて計算精度低下を防止することが可能である。
しかも、全ての化学種について計算量の多い反応速度論を用いるわけではないので、高速な演算が可能である。
【発明の効果】
【００１１】
よって、本発明によれば、化学平衡法を用いつつも、非平衡な燃焼状態における計算精度の悪化を防止することができる。
【発明を実施するための最良の形態】
【００１２】
以下、本発明の実施形態を図面に基づいて説明する。
【００１３】
［１．燃焼解析装置の全体構成］
図１は、実施形態に係る燃焼解析装置による燃焼解析処理の手順を示すフローチャートである。この燃焼解析装置は、図１のフローチャートの各ステップを実行するためのコンピュータプログラムが搭載されたコンピュータによって構成されている。
また、図２及び図３は、上記コンピュータプログラムによって実現される燃焼解析装置の機能ブロックを示している。
【００１４】
この燃焼解析装置は、解析対象の燃焼（乱流拡散燃焼）に係わる複数の化学種それぞれの生成量を算出することができる。燃焼解析装置は、解析対象の燃焼に係わる複数の化学種のうち、１又は複数の主要な化学種については、化学平衡法に基づくアルゴリズムに従って演算を行う。
一方、残りの一又は複数の化学種については、反応速度論（化学動力学）に基づくアルゴリズムに従って演算を行う。
【００１５】
このため、燃焼解析装置は、図２に示すように、解析対象の燃焼に係わる複数の化学種のうち、化学平衡法に従って演算が行われる１又は複数の化学種を「第１化学種」として設定するための第１化学種設定部１１と、反応速度論に従って演算が行われる１又は複数の化学種を「第２化学種」として設定するための第２化学種設定部２１と、を備えている。
したがって、燃焼解析装置のユーザは、設定部２１によって、化学平衡法ではなく、反応速度論で演算をしたい化学種（第２化学種）を自由に設定することができる。
【００１６】
また、燃焼解析装置は、設定（選択）された第１化学種それぞれについて、化学平衡法に基づくアルゴリズムに従って、第１化学種の質量分率（燃焼生成量）を算出する第１化学種の質量分率算出部１２を備えている。本実施形態において第１化学種の質量分率算出部１２は、化学平衡論（ＣＥ；Ｃｈｅｍｉｃａｌ Ｅｑｕｉｌｉｂｒｉｕｍ）の一種であるギブス自由エネルギー極小化法によって燃焼生成物を求める演算を行う。
【００１７】
ギブス自由エネルギー極小化法は、燃焼生成物を求める際に、平衡状態ではギブス自由エネルギーが極小値をとるという条件から元素の化学ポテンシャルを直接求めるものである。
化学平衡論で燃料生成物を求める方法としては、ギブス自由エネルギー極小化法以外に、圧力保存と原子数保存を拘束条件として平衡定数から各成分の濃度比を決める方法もあるが、この方法では、例えば、１００成分からなる燃焼場を解く場合、１００元の多元連立方程式を解く必要がある。
一方、ギブス自由エネルギー極小化法では、考慮する元素の数＋１の多元連立方程式を解けばよく、一般的な燃焼場では、Ｈ，Ｃ，Ｏ，Ｎの４元素を考慮すればよい。
【００１８】
さらに、燃焼解析装置は、設定（選択）された第２化学種それぞれについて、反応速度論に基づくアルゴリズムに従って、第２化学種の反応速度を算出する第２化学種それぞれの反応速度算出ｂ２２を備えている。本実施形態において第２化学種の反応速度算出部２２は、反応速度論に基づく燃焼モデルである渦消散コンセプト（ＥＤＣ；Ｅｄｄｙ Ｄｉｓｓｉｐａｔｉｏｎ Ｃｏｎｃｅｐｔ）モデルに従って反応速度を求める演算を行う。
なお、渦消散コンセプトモデル（ＥＤＣモデル）は、一段又は二段の総括反応の反応速度を算出することが可能な渦消散モデルを、多段階の化学反応式が解けるようにしたものである。
【００１９】
また、燃焼解析装置は、反応速度算出部２２によって求めた第２化学種それぞれの反応速度に基づいて、第２化学種それぞれの質量分率（燃焼生成量）を算出する第２化学種の質量分率算出部２３を備えている。
この第２化学種の質量分率算出部２３では、第２化学種の保存式に、前記反応速度を代入することで、燃焼生成物を求める。
【００２０】
また、渦消散コンセプトモデル（ＥＤＣモデル）は、ファイルスケールの体積分率算出部２４、タイムスケール算出部２５、原子の生成／消滅速度算出部２６も備えており、各値を算出することができる。
【００２１】
上記のように、本実施形態では、燃焼モデルが、化学平衡論（ＣＥ）と渦消散コンセプトモデル（ＥＤＣモデル）とによって構成されている。以下では、本実施形態の燃焼モデルを「ＥＤＣ／ＣＥモデル」とよぶ。
【００２２】
さらに、燃焼解析装置は、図３に示すように、燃焼解析のためのその他の演算を行う制御部３０を備えている。この制御部３０には、算出された第１化学種及び第２化学種の質量分率が与えられ、燃焼解析のための各種パラメータが算出される。制御部３０によって算出された各種パラメータは、第１化学種及び第２化学種の質量分率等の算出に用いられる。
この制御部３０は、処理ループを終了させるための収束判定部３１、速度と圧力を求める速度・圧力算出部３２、温度を算出する温度算出部３３、レイノルズ応力・乱流粘性係数を求めるレイノルズ応力・乱流粘性係数算出部３４、原子Ｉの質量分率を求める原子質量分率算出部３５、密度を求める密度算出部３６を備えている。
【００２３】
なお、図１のフローチャート及び図２及び図３ブロック図に関連して後に説明する数式は、上記コンピュータの記憶部に格納されたデータ（数式中の各項の値（パラメータ）に対応）に基づき、上記コンピュータの演算部によって演算され、その演算結果は前記記憶部に格納される。
【００２４】
［２．燃焼解析方法］
本実施形態に係る燃焼解析装置は、ＡＮＳＹＳ社の汎用流体解析ソフトウェアＦＬＵＥＮＴ（バージョン６．３．２６）に対して、本発明の機能をユーザ定義関数（User Defined Function）として組み込むことで実現した。ただし、本発明が適用されるソフトウェアがＦＬＵＥＮＴに限定されるものではない。
【００２５】
［２．１ 支配方程式］
解析対象の燃焼を支配する方程式としては、運動方程式、レイノルズ応力モデル、化学種Ｊの保存式、エネルギー方程式、気体の状態方程式、温度の算出式があり、具体的には下記の通りである。下記の支配方程式は、燃焼の数値解析の際に演算される。
【００２６】
［２．１．１ 運動方程式］
運動方程式は、下記の通りである。前記ＦＬＵＥＮＴでは、以下の連続の式とナビエストークスの式をSIMPLE 法［Suhas V Patanka. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, 1980 参照］により解く
【数１】
【数２】
【００２７】
［２．１．２ レイノルズ応力モデル］
本実施形態での解析対象の燃焼は、乱流であるため、乱流モデルとして、以下の式で表されるレイノルズ応力モデルを用いた。
【数３】
【００２８】
【数４】
【数５】
【数６】
【数７】
ここで、σkの推奨値は０．８２である。
【００２９】
θij,は、以下の式で与えられる。
【数８】
【００３０】
ここで，θij,1，θij,2 はそれぞれスローターム，ラピッドタームと呼ばれ，それぞれ以下のようにモデル化される。
【数９】
【数１０】
ここで、Ｃ1とＣ2はモデル定数であり、推奨値はそれぞれ１．８，０．６０であり、Ｐ＝（１／２）×Ｐkk，Ｃ＝（１／２）×Ｃkkである。
【００３１】
εijは、以下の式で与えられる。
【数１１】
【数１２】
【数１３】
ここで、Ｃμの推奨値は０．０９である。
なお、μは、粘性係数であり、流体が持つ固有の物性値である。一方、μtは、乱流によって生じる拡散の強さを表す指標（乱流による拡散率）であり、乱流モデルによって求められる。
【００３２】
［２．１．３ 化学種Ｊの保存式］
燃焼に係わる化学種Ｊの質量分率の式はＦＬＵＥＮＴにより以下の式で計算される。
【数１４】
【数１５】
【数１６】
【００３３】
［２．１．４ エネルギー方程式］
エンタルピーは、ＦＬＵＥＮＴのユーザ定義スカラー（ＵＤＳ）によって計算される。エネルギー方程式は下記の通りである。
【数１７】
【数１８】
【数１９】
【００３４】
［２．１．５ 気体の状態方程式］
気体の密度は、次式によって計算される。
【数２０】
【００３５】
［２．１．６ 温度の求め方］
温度は、次式で、ＦＬＵＥＮＴのＵＤＳにより計算される。
【数２１】
ここで、Ｔ0は参照温度で、２９８．１５Ｋとし、ＣpmJは化学種Ｊの平均定圧比熱、ｈJ0は化学種Ｊの標準生成エンタルピーを指す。
【００３６】
［２．２ ＥＤＣ／ＣＥモデルの説明］
本実施形態のＥＤＣ／ＣＥモデル（渦消散コンセプト・平衡法モデル）は、下記の通りである。なお、以下では、解析対象をＨ2−ａｉｒ乱流拡散火炎とし、Ｈ2をＥＤＣモデルにて計算した場合の例に基づいて説明する。
【００３７】
ａ．燃料（ここでは、Ｈ２とＮ２の混合気）の反応速度をＥＤＣデルで計算する。
ｂ．反応後の燃料は化学平衡に達していると仮定し、これを反応物として平衡計算をする。
【００３８】
前記反応速度は、ＥＤＣ（渦消散コンセプト）モデルの式で計算する。ＥＤＣモデルでは、燃焼は、乱流の微細構造、すなわちファインスケールで起こるものと仮定する。ファインスケールの体積分率は、次式で定義される（I R Gran and B F Magnussen. A numerical study of a bluff-body stabilized diffusion flame .part2. influence of combustion modeling and finite-rate chemistry. Combustion Science Technology, Vol.119, pp. 191_217, 1998. 参照）。
【数２２】
【００３９】
さらに、反応は次式で表されるタイムスケールで進行するものとする。
【数２３】
【００４０】
さらに、τ*のタイムスケールで反応が定圧反応容器内で起こるものと仮定し、格子内の現在の化学種と温度を初期条件とし、次の常微分方程式をｔ＝０〜τ*で数値積分する。なお、上記の「格子内」とは、数値解析の対象空間を格子状に分割した小さな空間をいう。数値流体力学では、その分割したそれぞれの空間（格子）に対して、運動方程式や例のずる応力などの方程式を解いていく。この格子は、セルとも呼ばれる。
【数２４】
【数２５】
【数２６】
ここで、Ｙ*J ，Ｒ*J はそれぞれファインスケール内の化学種Ｊの質量分率と反応速度を指し、ｖ'J K とν''J Kは、それぞれ化学種Ｊの化学反応式Ｋにおける正と負の量論係数、ＸJは化学種Ｊのモル濃度、ｋKは化学反応式Ｋの反応速度定数，ＫMAXは考慮する反応速度式の数である。Ａは頻度因子，Ｅは活性化エネルギー、βはそれぞれ温度の指数で実験パラメータである．また，式（２６）はアレニウスの式と呼ばれ，この式を含む常微分方程式は一般にstiff であることが知られている．この常微分方程式を数値積分するために，本実施形態ではＣＶＯＤＥ[S D Cohen and A C Hindmarsh. CVODE User Guide, llnl report edition, 1994. UCRL-MA-118618. 参照] を用いた。
【００４１】
実際の化学種Ｊの保存式（１４）に現れる反応速度は、次式でモデル化される。
【数２７】
【００４２】
また、燃料中のＮ2の反応速度は、Ｈ2に比例すると仮定し、燃料中のＮ2の質量分率及びＨ2の質量分率に基づき、次式で計算される。つまり、本実施形態における第２化学種の反応速度算出部２２は、燃料中のＮ2については、式（２８）で計算し、燃料中のＮ2以外の化学種については、式（２７）で反応速度を計算する。
【００４３】
【数２８】
【００４４】
上記の式（２７）は、詳細化学反応式を取り入れる場合である。しかし、化学反応式を取り入れず、反応は乱流混合が律速すると仮定すると、化学反応式は総括反応で表される。Ｈ2−ａｉｒを例にとると、Ｈ2＋１／２Ｏ2−＞Ｈ2Ｏとなり、このときの燃料であるＨ2の反応速度は以下の渦消散モデルで計算することができる。
【数２９】
ここで、Ｙfuel、Ｙox、Ｙpはそれぞれ燃料、酸化剤、燃焼生成物の質量分率を指し、ＡmixとＢmixはそれぞれ渦消散モデルの実験パラメータでそれぞれ４．０と０．５である。stは量論酸素燃料質量比を指す。
【００４５】
原子Ｉの生成または消滅速度は次式で計算される。
【数３０】
さらに、ａIJは化学種Ｊに対する原子Ｉの数の行列を示す．例えば，化学種Ｈ２に対する原子Ｈの数はａH,H2＝２であり、化学種ＯＨに対する原子Ｈの数はａH,OH＝１、原子Ｏの数は、ａO,OH＝１となる。
【００４６】
原子Ｉの保存式は次式で表される。
【数３１】
【００４７】
このときのエンタルピーは次式で近似される。
【数３２】
ＩMAXは、燃焼に関与する原子の数である。通常の燃焼場では、Ｈ，Ｃ，Ｏ，Ｎだけでよい。つまり、ＩMAX＝４となる。
【００４８】
本実施形態のＥＤＣ／ＣＥモデルでは、さらにＮを調節するパラメータを有する。
【数３３】
このＮＣＴＲＬの値を変化させることで、平衡計算による計算結果を調整できる。
平衡計算のためのプログラムは、ギブス自由エネルギー極小化法［JSME, editor. JSME Combustion Handbook. JSME, 2002. 参照］の理論に基づいたＮＡＳＡのＣＥＡ［S Gordon and B Mc Bride. Computer program for calculations of complex equilibrium compositions and applications. NASA Reference Publication, No. 1311, 1994. 参照］をベースに作成し、これをＵＤＦとしてＦＬＵＥＮＴに組み込んだ。
【００４９】
［２．３ 処理手順］
図１に戻り、燃焼解析の処理手順を説明する。ここでは、後述の図６の「Ｃａｓｅ３」の例に則して処理手順を説明する。なお、図１のフローチャートにおいて、ステップａ，ｂ，ｄ，ｅ，及びｇの部分が、ＥＤＣ／ＣＥモデルである。また、ステップａがＥＤＣモデルの部分、ｂがＣＥの部分、ｄ，ｅ，ｇがそれ以外で化学種の保存式や温度計算などを行っている部分である。その他のステップｃ，ｆ，ｈ，ｉは、ＦＬＵＥＮＴから提供される機能を利用した。
【００５０】
このＥＤＣ／ＣＥモデルにおいて考慮する化学種は、Ｈ2，Ｏ2，Ｈ2Ｏ，ＯＨ，Ｈ，Ｏ，ＨＯ2，Ｎ，Ｎ2である。そして、第１化学種設定部１１では、これらの化学種のうち、第１化学種として、Ｈ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2が設定（選択）されているものとする。また、第２化学種設定部２１では、残りの化学種であるＨ2，Ｈ，Ｏ，ＯＨが第２化学種として設定（選択）されているものとする。
【００５１】
ステップａ：
ステップａでは、ＥＤＣモデルの反応速度算出部２２によって、第２化学種であるＨ2，Ｈ，Ｏ，ＯＨ（Ｊ＝Ｈ2，Ｈ，Ｏ，ＯＨ）の反応速度ＲJを求める。この演算は、式（２４）〜（２９）に基づいて行われる。
また、ＥＤＣモデルの原子の生成／消滅速度算出部２６では、上記の反応速度から、式（３０）に基づき、原子Ｉの生成または消滅速度を求める。
また、ＥＤＣモデルでは、ファインスケールの体積分率算出部２４が、式（２２）に基づき体積分率を求め、タイムスケール算出部２５が、式（２３）に基づきタイムスケールを求める。
【００５２】
ステップｂ：
ステップｂでは、ＣＥ（ギブス自由エネルギー極小化法）の質量分率算出部１２によって、第１化学種であるＨ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2についての平衡計算を行う。この平衡計算では、原子Ｉ（Ｉ＝Ｃ，Ｈ，Ｏ，Ｎ）の質量分率ｄIを燃焼反応物とし、燃焼前後で式（３２）のエンタルピーは変化しないことを条件に、第１化学種であるＨ2Ｏ，Ｏ2，ＨＯ2，Ｎ，Ｎ2の質量分率（生成量）ＹH2O，ＹO2，ＹHO2，ＹN，ＹN2を求める。
【００５３】
ステップｃ：
ステップｃでは、速度・圧力算出部３２が、ＳＩＭＰＬＥ法（式（１）〜（２））に基づき、速度と圧力を求める。
【００５４】
ステップｄ：
ステップｄでは、温度算出部３３が、エネルギー方程式（式（１７））を解き、これから求まるエンタルピーから、式（２１）に基づき、温度を求める。
【００５５】
ステップｅ：
ステップｅでは、第２化学種の質量分率算出部２３が、ステップａで求めた第２化学種の反応速度ＲJそれぞれを、対応する第２化学種Ｊの保存式（式（１４））のＲJに代入し、第２化学種であるＨ2，Ｈ，Ｏ，ＯＨの質量分率（生成量）ＹH2，ＹH，ＹO，ＹOHを求める。
【００５６】
ステップｆ：
ステップｆでは、レイノルズ応力・乱流粘性係数算出部３４が、レイノルズ応力モデルに基づき、レイノルズ応力と乱流粘性係数を求める（式（３）〜（１３））。なお、レイノルズ応力は、運動方程式に現れる。また、乱流粘性係数は、エネルギー方程式、化学種Ｊの保存式、原子Ｉの保存式に現れる。
【００５７】
ステップｇ：
ステップｇでは、原子質量分率算出部３５が、ステップａで求めた原子Ｉ（Ｉ＝Ｃ，Ｈ，Ｏ，Ｎ）の生成／消滅速度ＲIそれぞれを、対応する原子Ｉの保存式（式（３１））に代入し、原子Ｉの質量分率ｄIを求める。
【００５８】
ステップｈ：
ステップｈでは、密度算出部３６が、化学種Ｊの質量分率ＹJと圧力から、式（２０）に基づき、密度を求める。
【００５９】
ステップｉ：
ステップｉでは、収束判定部３１が、図１のステップａ〜ｈの処理ループを終了させるために、収束判定を行う。
収束判定は、運動方程式、レイノルズ応力モデル、エネルギー方程式、化学種Ｊの保存式、原子Ｉの質量分率の式が正しい解に到達したかを判定する誤差判定である。
この収束判定は、ＦＬＵＥＮＴが本来有する機能によって自動的に実行される。ＦＬＵＥＮＴが採用している有限体積法をベースとしているＳＩＭＰＬＥ法では、離散化後、セルＰにおける任意の変数の保存式は、次式のように表される。
【数３４】
ここで、ａpは中心係数、ａnbは隣接セルの各影響係数、ｂは生成項Ｓ＝Ｓc＋ＳPφの式中の定数部分Ｓcを指す。式（３４）では次式が成り立つ。
【００６０】
【数３５】
【００６１】
さらに、このときのφの残差は次式で表される。
【数３６】
ＦＬＵＥＮＴでは、各変数φについて、この残差が一定の値以下（通例で、１．０×１０-3〜１．０×１０-6）になるまで、計算を続ける。
【００６２】
［３．実施例と比較例の計算条件］
本実施形態のＥＤＣ／ＣＥモデルによる計算精度を確かめるために、水素噴流火炎の計算（実施例と比較例）を行った。図４に計算領域と境界条件を示す。計算条件は高城・古藤の水素噴流火炎の実験[Toshimi Takagi and Satoru Kotou. A prediction of flow and combustion in turbulent diffusion flame in japanese. Transactions of the Japan Society of Mechanical Engineers, Series B, Vol. 48, No. 436, 1981. 参照] に準ずるものとした。
半径が２．４５ｍｍの円筒形ノズルから水素と窒素のモル比が０．４：０．６の燃料が５５．４ｍ／ｓで噴出されるものとし、周囲流は空気とした。空気については初めから平衡状態であるとし、モル比でｄＯ：ｄＮ＝０．２１：０．７９を入り口条件として与える。
【００６３】
さらに，燃料と周囲流の温度はともに３００Ｋとし、計算結果を比較する場所は、中心軸上とノズルから６０ｍｍと１６０ｍｍの地点の断面上とした。全ての方程式に二次精度の風上差分を適用し、図４の中心線を軸とした軸対象流れとし、乱流モデルとしてレイノルズ応力モデルを使用した。
【００６４】
図５に今回実行した計算（実施例と比較例）の代表的なシミュレーション条件（Ｃａｓｅ１，３，４）を示す。なお、Ｃａｓｅ２は欠番である。
また、図６に、上記計算で使用した化学反応式を示す。
このＥＤＣ／ＣＥモデルにおいて考慮する化学種は、Ｈ2，Ｏ2，Ｈ2Ｏ，ＯＨ，Ｈ，Ｏ，ＨＯ2，Ｎ，Ｎ2である。
【００６５】
図５のＣａｓｅ１（比較例）では、燃焼モデルとして従来のＥＤＣモデルを採用した。
このＣａｓｅ１のＥＤモデルでは、図６に示す化学反応式を用いた。
Ｃａｓｅ３（実施例）では、Ｈ2，Ｈ，Ｏ，ＯＨの化学種（第２化学種）のみ反応速度を考慮し、残りの化学種は第１化学種とし、平衡法で計算をした。
Ｃａｓｅ４（実施例）では、ＯＨが関わる図６の１〜８，１１，１３の化学反応式を考慮し、さらに、Ｈ2については、混合律速を仮定し、式（２９）で評価した。
【００６６】
［４．実施例及び比較例の計算結果と実験値との比較］
［４．０．１ 乱流エネルギーｋ、平均速度Ｕ、主要化学種の中心軸上での比較］
図７に、図４の中心線（中心軸）上における乱流エネルギーｋ、平均流速Ｕの分布を示し、図８に化学種と温度Ｔの分布を示す。
図７において、乱流エネルギーｋを見ると、実験値（ＥＸＰ）ではＬ＝１８０あたりまで緩慢に増えるのに対して、Ｃａｓｅ１，３では，Ｌ＝１００付近に向けて急激に増えて、ここを境に減少している。過去の計算例[Toshimi Takagi and Satoru Kotou. A prediction of flow and combustion in turbulent diffusion flame in japanese. Transactions of the Japan Society of Mechanical Engineers, Series B, Vol. 48, No. 436, 1981.][ Kimio Iino, Shinji Murakami, and Hoshitaka Kikkawa. Numerical simulation of high pressure methane-oxygen coaxial turbulent nonpremixed flames. Journal of High Temperature Society, Vol. 31, pp. 112_121, 2005.] においても、平均流速Ｕの実験値との整合性は良いが、ｋについてはシミュレーション値が過大評価されることが報告されている．この乱流拡散火炎のｋの過大評価を改善するためには，より高精度な乱流モデルのＬＥＳモデルを選択することや，レイノルズ応力モデルを改善する必要がある。
【００６７】
ここでは、全てのケースで中心軸上で温度の最大値の位置が一致するように図５のＲＳＭのパラメータを選んだ。実験ではノズルから１５０ｍｍ付近で温度は最大値の約１６００Ｋをとる。Ｃａｓｅ３のＨ2Ｏと温度の最大値は、それぞれ１８．９％、１５５６Ｋとなっており，実験値（ＥＸＰ）でＨ2Ｏと温度の最大値の約１８．３５％、約１６００Ｋとなっている。
【００６８】
［４．０．２ 乱流エネルギー、平均速度、主要化学種，温度のＬ＝６０での比較］
図９に、Ｌ＝６０での半径方向のｋ，Ｕ、図１０に、化学種と温度Ｔの分布を示す。図９のＣａｓｅ３、Ｃａｓｅ１のＵは、実験値（ＥＸＰ）でのＵと良く一致している。しかし、ｋについては中心軸上の分布と同じく実験値（ＥＸＰ）と差異が見られる。図１０ を見ると、Ｒ＝６ｍｍ付近の温度、Ｈ2ＯについてＣａｓｅ３、Ｃａｓｅ１ともにＥＸＰとの誤差が大きい。Ｒ＝６ｍｍ付近ではＴ，Ｈ2ＯのＥＸＰとの誤差は大きいがそれ以外の領域ではＥＸＰとの一致が良好なためである．
また、Ｃａｓｅ１，３ともにＲ＝６ｍｍ付近の温度を正しく予測できておらず，両者ともに似た予測値を示していることから，この付近のＴの誤差は燃焼モデルを起因とする誤差ではないと推測できる。
【００６９】
［４．０．３ 乱流エネルギー、平均速度、主要化学種、温度のＬ＝１６０での比較］
図１１にＬ＝１６０での半径方向のｋ，Ｕ、図１２に化学種と温度Ｔの分布を示す。
図１１を見ると、Ｒが大きくなるにつれて、Ｃａｓｅ１，３のＵの予測値と実験でのＵとの間の差異が大きくなる。この差異が存在する部分で拡散が大きく評価されれば、Ｕは拡散するので図１１中のＣａｓｅ１，３によるＵは小さくなると予測できる．このことから、この部分での差異は乱流拡散を実験よりも小さく見積もっていることが原因と推測される。さらにこの差異を改善するためには，より高精度な乱流モデルのＬＥＳモデルを選択することや、レイノルズ応力モデルを改善すればよい。また、図8１２を見ると、Ｃａｓｅ１，３の温度をＥＰＸと比較するとＲが大きくなるにつれて（火炎の外側へいくにつれて）Ｃａｓｅ１，３のＴの予測値と実験でのＴとの間の差異が大きくなる。これは、上記のＵ と同じく拡散を小さく見積もっていることと、さらにいえば、Ｕが大きいためＨ2Ｏが流されてしまったことが原因と考える。
【００７０】
［４．１ ラジカル類のEDC モデルとの比較］
［４．１．１ 中心軸上］
図１３，１４に中心軸上でのＯＨとＨのＥＤＣ／ＣＥモデルとＥＤＣモデルによる予測値を表す。図５にまとめたとおり、図５のＣａｓｅ１はＥＤＣモデル（比較例）、Ｃａｓｅ３，４はＥＤＣ／ＣＥモデル（実施例）での計算値を表す。Ｃａｓｅ３，４の違いはＥＤＣモデルで考慮する化学種の数の違いである。Ｃａｓｅ３ はＨ2，Ｈ，Ｏ，ＯＨの反応速度をＥＤＣモデルで考慮しているのに対し、Ｃａｓｅ４はＯＨのみを考慮している。
図１３に示すように、ＯＨの予測値を見ると、Ｃａｓｅ１，３，４との計算結果を比較するとＬ＝５０付近とＬ＝２００付近で差異が見られるが，それ以外の領域では整合性は良い．
これは、Ｃａｓｅ３，４ともにＯＨが式（２２）〜（２７）で示されるＥＤＣモデルの理論で計算されているからである。
【００７１】
図１４に示すように、Ｈの予測値ではＬ＝１５０まではＣａｓｅ１とＣａｓｅ３が良く一致している。しかし、Ｌ＝１５０以降はＣａｓｅ３のＨの予測値がＣａｓｅ１による予測値より小さく見積もるという差異が見られる。この差異は、Ｈ2Ｏ，Ｏ2，ＨＯ2の反応速度をＥＤＣモデルで考慮しなかったことが原因と考えられる。Ｃａｓｅ１とＣａｓｅ４ を比較すると、Case４ではＨをＣａｓｅ１より低く見積もっている。
【００７２】
［４．１．２ Ｌ＝６０］
図１６，１５にＬ６０でのＯＨとＨのＥＤＣ／ＣＥモデルとＣａｓｅ１のＥＤＣモデルによる予測値を表す。
図１６に示すＯＨの予測値を見ると、Ｃａｓｅ１とＣａｓｅ３，４との計算結果を比較するとＲ＝７付近まで良く一致しているが、それ以降ではＥＤＣ／ＣＥモデルとＥＤＣモデルの予測値に若干の差異が見られる。概ねＣａｓｅ１，３，４ともに近い予測値を示しているといえる。これは、Ｃａｓｅ３，４ともにＯＨについて反応速度を考慮したからである。
図１５に示すＨの予測値では、Ｃａｓｅ１とＣａｓｅ３を比較すると、Ｒ＝０付近とＲ＝６以降で若干の誤差が見られる。しかし、概ね両者の一致は良好である。
Ｃａｓｅ４のＨの予測値はＣａｓｅ１に比べ、Ｒ＝０付近では多く見積もり、Ｒ＝８付近を堺にしてＣａｓｅ１の予測値より低く見積もるようになる。
【００７３】
［４．１．３ Ｌ＝１６０］
図１８，１７にＬ＝１６０でのＯＨとＨのＥＤＣ／ＣＥモデルとＥＤＣモデルによる予測値を示す。
図１８のＯＨの予測値で、Ｃａｓｅ１，３，４との計算結果を比較すると、３者は似たグラフの傾きを示しているといえる。
図１７のＨの予測値では、Ｃａｓｅ１，３は良く一致している。Ｃａｓｅ４のＨの予測値はＣａｓｅ１に比べ、Ｒ＝０付近では多く見積もり、Ｒ＝１３付近を堺にしてＣａｓｅ１の予測値より低く見積もるようになる。ＨはＲ＝１７付近で完全に消費されている。
【００７４】
［５ 考察］
［５．１ 実施例の燃焼モデルと従来の燃焼モデルとの対比］
従来のＥＤＣモデルは乱流拡散火炎において詳細な化学反応式を考慮することができ、多成分系の反応を取り扱える。一方で、実施例に係るＥＤＣ／ＣＥモデルは化学反応式と化学平衡法の両方を用いて多成分系の反応を取り扱う。ここでは、ＥＤＣモデルと本実施例のモデルの計算値を比較し，両者の特徴について議論する。
【００７５】
ＥＤＣモデルでは，特定の化学種の発生量を調節するためには化学反応式の数を増やしたり、化学反応式中の実験パラメータとＥＤＣモデルのＣτ ，Ｃξ などのパラメータを調節する。しかし、これらのパラメータの組み合わせは多く、化学反応式の数を増やすとパラメータの数も増える。一方、実施例のモデルはＥＤＣモデルのパラメータと任意の化学反応式を考慮するだけで良いので、調節はＥＤＣデルより簡便である。
【００７６】
また、ＥＤＣ／ＣＥモデルでは考慮すべき化学反応式はシミュレーションの実行者の要求する計算精度に従う。
例えば、Ｃａｓｅのように任意の化学種（ここではＯＨ) について計算精度が必要であれば、その化学種について反応速度で考慮する。また、Ｃａｓｅ３のように他の化学種についても計算精度が必要ならば，他の化学種についても化学反応式を考慮する。こうすることで、ＥＤＣ／ＣＥモデルの計算結果はＥＤＣモデルの計算結果に近づけることができる。ＥＤＣ／ＣＥモデルでは、仮に平衡法で化学種を一つも考慮せず，全ての化学種をＥＤＣモデルで考慮したとすると、EDC モデルと同じになる。Ｃａｓｅ４で計算したラジカル類はEDC モデルの計算結果と比較すると、ＯＨについては良い一致を見せた。しかし、Ｈについては定量的な一致は見せなかった。
【００７７】
とはいえ、図１５，１７ を見ると、Ｈについても発生や消滅の傾向は捉えられていると考える。このことから、ＥＤＣ／ＣＥモデルの使い方としては，数値解析の初段階ではＣａｓｅ４のように必要な化学種のみにＥＤＣモデルによる反応速度計算を適用し、その後、他の化学種について精度が必要となれば、Ｃａｓｅ３のようにＥＤＣモデルで考慮する化学種を増やしていけばよい。ＥＤＣ／ＣＥモデルの適用としては，例えば、サーマルＮＯxを次の化学反応式で表される拡大ゼルドヴィッチ機構にて計算する際にＨ，Ｏ，ＯＨ については平衡法で考慮し、ＮＯxについてはＥＤＣモデルで考慮するといった使い方ができる。
【００７８】
【数３７】
【数３８】
【数３９】
【００７９】
他のＥＤＣ／ＣＥモデルの適用としては，乱流拡散火炎の詳細化学反応式を構築するためにも有用である。すなわち、Ｃａｓｅ４のように燃料の反応速度にだけＥＤＣモデルを適用し、他の化学種については平衡法を適用する。その後、燃焼場への寄与が大きいと推定される化学種について化学反応式を追加していく。燃焼場への寄与の大きさについては、平衡法で求まる予測値を参照すれば推測できる。
【００８０】
［５．１．１ 計算時間について］
ＥＤＣ／ＣＥモデルにおいて、平衡法では考慮する元素の数＋１を未知数とした非線形の連立方程式を解く。通常の燃焼場ではＣ，Ｈ，Ｏ，Ｎ＋１で５元となる。一方でＥＤＣモデルは考慮する化学種の数を未知数とした常微分連立方程式を解く。化学反応式が増え、考慮する化学種が増えるにつれ、常微分の連立方程式を解くための計算時間も増える。それに対して、平衡法では常に５元であり、反応物の濃度によって収束までの時間に差異はあるものの、EDC モデルより計算時間の増減は顕著ではない。また、ＥＤＣ／ＣＥモデルでは、平衡法で考慮する化学種については化学種Ｊの保存式（１４）を計算する必要がない。ＥＤＣモデルでは考慮する全ての化学種について化学種Ｊの保存式を解く必要がある。このために、ＥＤＣモデルでは考慮する化学種が増えるにつれ、解くべき化学種の保存式の数も増えるため計算時間は増加する。このような理由から、計算に考慮する化学種の数が多ければ多いほど、ＥＤＣ／ＣＥモデルはＥＤＣモデルより高速であるといえる。
【００８１】
［５．１．２ ＥＤＣ／ＣＥの適用範囲］
ＥＤＣ／ＣＥモデルの適用できる燃焼場は、例えば、乱流拡散火炎である。つまり、流れ場は乱流であり、燃料と酸化剤が別々の場所から流入する。また、反応速度の計算にＥＤＣモデルを使用しているため化学種の予測精度は最も良い精度でもＥＤＣモデルと同じとなる。ＮＯxなどの非平衡性の強い化学種に対して平衡法を適用すると予測制度が低下するので、このような化学種についてはEDC モデルで反応速度を計算しなければならない。
【００８２】
また、最終的に生じる生成物、例えばＨ2ＯやＣＯ2については平衡法で考慮する。燃焼によって主に生成する化学種を定めることで、平衡法の計算を安定化させるためである。言い換えれば、燃焼生成物の化学種の大部分を平衡法で計算し，精度が必要な化学種に限って補助的にＥＤＣモデルを取り入れることが好ましい。また、ＥＤＣ／ＣＥモデルではアレニウスの式を使っており、アレニウスの式は温度に依存するのでしばしば消炎が起こる。この問題については，燃料の反応速度についてアレニウスの式を考慮しない、式（２９）を使うと回避できる。
【００８３】
［６ 結論」
結論として、実施例に係るＥＤＣ／ＣＥモデルの特徴と、その比較・検討についてまとめる。
（１）詳細化学反応式を必要とせず全ての化学種の生成量を予測でき，任意の化学種について化学反応式を取り入れることでＥＤＣモデルと同程度の精度で計算することができる。
（２）Ｈ2Ｏなどの最終生成物だけでなく、ＯＨ，Ｏ，Ｈなどの中間生成物も予測することができ，これとゼルドビィッチ機構などの既存のモデルと組み合わせることでＮＯxなどの生成量の予測などの応用が可能である。
（３）乱流混合・燃焼をＥＤＣモデルにより評価しているので，化学反応式を考慮し化学反応に乱流と温度の影響を取り入れることができる。
【００８４】
さらに、実験値（ＥＸＰ）、従来のＥＤＣモデルによる予測値との比較・検討についてまとめる。
（１）Ｈ2，Ｏ2，Ｈ2Ｏについて、ＥＤＣ／ＣＥモデルによる予測値は、実験値とＥＤＣ モデルによる予測値と一致する。
（２）ラジカル類である、ＯＨ，Ｈについては反応速度を考慮することで、ＥＤＣモデルの予測値と良い一致を見せるようになる。
（３）多くの化学種について反応速度を考慮すると、ＥＤＣモデルによる予測値に近づく。
【００８５】
なお、本発明は、上記実施形態及び実施例に限定されるものではなく、様々な変形が可能である。
【図面の簡単な説明】
【００８６】
【図１】燃焼解析処理のフローチャートである。
【図２】燃焼解析装置の機能ブロック図である。
【図３】燃焼解析装置の機能ブロック図である。
【図４】水素噴流火炎の計算領域と境界条件を示す図である。
【図５】Ｈ２−Ｏ２化学反応メカニズムと関連するレート係数を示す表である。
【図６】シミュレーション条件を示す表である。
【図７】中心線に沿った位置での乱流エネルギーと平均速度を示すグラフである。
【図８】中心線に沿った位置での化学種のモル分率と温度を示すグラフである。
【図９】Ｌ＝６０の位置での乱流エネルギーと平均速度を示すグラフである。
【図１０】Ｌ＝６０の位置での化学種のモル分率と温度を示すグラフである。
【図１１】Ｌ＝１６０の位置での乱流エネルギーと平均速度を示すグラフである。
【図１２】Ｌ＝１６０の位置での化学種のモル分率と温度を示すグラフである。
【図１３】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、中心線に沿った位置でのＯＨのモル分率を示すグラフである。
【図１４】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、中心線に沿った位置でのＨのモル分率を示すグラフである。
【図１５】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝６０の位置でのＨのモル分率を示すグラフである。
【図１６】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝６０の位置でのＯＨのモル分率を示すグラフである。
【図１７】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝１６０の位置でのＨのモル分率を示すグラフである。
【図１８】ＥＤＣモデルとＥＤＣ／ＣＥモデルにおいて、Ｌ＝１６０の位置でのＯＨのモル分率を示すグラフである。
【符号の説明】
【００８７】
１１ 第１化学種設定部
１２ 第１化学種の質量分率算出部
２１ 第２化学種設定部
２２ 第２化学種の反応速度算出部
２３ 第２化学種の質量分率算出部
【特許請求の範囲】
【請求項１】
燃焼を数値解析する燃焼解析方法であって、
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、
を含むことを特徴とする燃焼解析方法。
【請求項２】
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する手段と、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する手段と、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する手段と、
を備えていることを特徴とする燃焼解析装置。
【請求項３】
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、
をコンピュータに実行させるためのコンピュータプログラム。
【請求項１】
燃焼を数値解析する燃焼解析方法であって、
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、
を含むことを特徴とする燃焼解析方法。
【請求項２】
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する手段と、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する手段と、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する手段と、
を備えていることを特徴とする燃焼解析装置。
【請求項３】
解析対象の燃焼に係わる化学種のうちの一部である１又は複数の第１化学種の生成量を、化学平衡法に基づくアルゴリズムに従って算出する第１ステップと、
前記燃焼に係わる化学種のうちの他の一部である１又は複数の第２化学種の前記燃焼による反応速度を、反応速度論に基づくアルゴリズムに従って算出する第２ステップと、
前記第２ステップによって求めた反応速度を、前記第２化学種の保存式に代入することで、前記第２化学種の生成量を算出する第３ステップと、
をコンピュータに実行させるためのコンピュータプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【公開番号】特開２０１０−３８４１０（Ｐ２０１０−３８４１０Ａ）
【公開日】平成２２年２月１８日（２０１０．２．１８）
【国際特許分類】
【出願番号】特願２００８−２００１２９（Ｐ２００８−２００１２９）
【出願日】平成２０年８月１日（２００８．８．１）
【出願人】（５９３００６６３０）学校法人立命館 (359)
【Ｆターム（参考）】
【公開日】平成２２年２月１８日（２０１０．２．１８）
【国際特許分類】
【出願日】平成２０年８月１日（２００８．８．１）
【出願人】（５９３００６６３０）学校法人立命館 (359)
【Ｆターム（参考）】
[ Back to top ]