(57)【要約】      （修正有）【課題】化学反応を伴う多成分系の速度場において、既存の硬い系専用に開発された１次元解法（Ｇｅｅｒの公式以外の解法も含む）を用いることによって、３次元解を簡単に近似できるようにする。【解決手段】解析領域内における速度場の任意の空間格子点Ｐでの濃度計算方法は、周りの６個の近傍点の各々の速度成分から格子点Ｐに流入する近傍点を見つけ、各成分ごとに既存の１次元解法を用いて、各々の１次元解を求め、これらの解を点Ｐに流入する速度成分によって加重平均し、この値をもって点Ｐでの濃度の近似値とする。また、これら一連の手続きで得られた全格子点での値を、それぞれ次回の全格子点での出発値とするような反復計算を有限回行ない、これらによって得られた最終値をもって、濃度に関する３次元連立１階偏微分方程式の近似解とすることが出来る。【選択図】図１
2005004709.tif 000002
【特許請求の範囲】
【請求項１】
解析領域内における速度場の任意の空間格子点Ｐに対して、周りの６個の近傍点の各々の速度成分を調べ、前記格子点Ｐに流入する近傍点の特定及び前記格子点Ｐに流入する速度成分を算出することを特徴とする計算方法。
【請求項２】
請求項１記載の計算方法において、各成分ごとに既存の１次元解法を用いて、各々の１次元解を求め、これらの解を前記格子点Ｐに流入する前記算出の速度成分によっで加重平均し、この値をもって、濃度に関する３次元連立１階偏微分方程式の、点Ｐにおける近似解とすること特徴とする計算方法。
【請求項３】
請求項２記載の計算方法において、流れに逆流がある場合、一連の手続きで得られた全格子点での値を、それぞれ次回の前記全格子点の出発値とするような反復計算を有限回行い、これらによって得られた前記全格子点での最終値をもって、前記偏微分方程式の近似解とすることができることを特徴とする計算方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】本発明は、化学反応を伴う多成分系の複雑な流体において、濃度分布を効率良く計算する方法に関するものである。
【０００２】
【従来の技術】数値流体力学では、現象を支配する偏微分方程式を解く数値解法の１つの手段として、汎用性の高い差分法を用いることが多い。複数の化学種を含む化学反応解析では、解析特有の硬い系（“ｓｔｉｆｆな系”）の問題が生じる。この問題は１階の連立微分方程式で表され、１次元の場合、現在までいろいろ研究され、これを解くための解法が既に幾つか存在している。この中では、Ｇｅａｒの公式が特に有名である。現在、この公式を用いたプログラムが文献やインターネット上で数多く見うけられる。
【０００３】
【発明が解決しようとする課題】しかしながら、２次元以上、特に３次元の場合、１次元のＧｅａｒの公式に相当する公式は、現在のところ見当たらない。
【０００４】本発明は、与えられた、化学反応を伴う多成分系の速度場において、既存の硬い系専用に開発された１次元解法（Ｇｅａｒの公式以外の解法も含む）を用いることによって、３次元解を簡単に近似できるようにすることを目的とする。
【０００５】
【課題を解決するための手段】解析領域内における速度場の任意の空間格子点Ｐでの濃度計算方法は、周りの６個の近傍点の各々の速度成分から格子点Ｐに流入する近傍点を見つけ、各成分ごとに既存の１次元解法を用いて、各々の１次元解を求め、これらの解を点Ｐに流入する速度成分によっで加重平均し、この値をもって点Ｐでの濃度の近似値とすること特徴とする。また、これら一連の手続きで得られた全格子点での値を、それぞれ次回の全格子点での出発値とするような反復計算を有限回行い、これらによって得られた最終値をもって、濃度に関する３次元連立１階偏微分方程式の近似解とすることができることを特徴とする。
【０００６】
【発明の実施の形態】以下、この発明の実施形態を説明する前に、本発明の概要を説明する。
【０００７】本発明では、解析領域内の速度場に関して、次の２つの仮定を設ける。仮定１：流れは解析領域の入口から出口へ向かい、解析領域内部では、湧き出しは無いものとする。仮定２：解析領域内の隣り合う２つの空間格子点において、互いに向きあうような流れの速度成分がある場合、その速度成分に関しては、２点間で流れは相殺される、即ち、流れは無いものとして取り扱う。
【０００８】更に、流れが計算の進行方向に対して、順方向か逆方向かを明確にするために、次の２つの定義を与える。定義１：速度場の格子点Ｐにおいて、Ｐに流入する速度成分の方向がすべて計算の進行方向と一致する場合、点Ｐで順流という。定義２：前記定義１の点Ｐにおいて、Ｐに流入する速度成分のうちで、計算の進行方向と逆方向の成分がある場合、点Ｐでは逆流成分をもつという。
【０００９】化学反応を伴う多成分系の速度場では、流れは、レイノルズ数が十分高い場合、化学種ｉの濃度Ｕｉは数１の３次元連立１階偏微分方程式で計算される。
【００１０】
【数１】
2005004709.tif 000003
【００１１】ここで、Ｖｘ，Ｖｙ，Ｖｚは其々ｘ，ｙ，ｚ方向の速度成分で、ｆｉはＵｉの非線形関数である。また、Ｎは化学種の総数である。今後、特に断らない限り添字ｉは数１と同じ意味に用いる。
【００１２】しかしながら、数１の離散化に、直接３次元差分法を適用すると、硬い方程式のため、空間メッシュを相当細かくしても計算途中で解が発散してしまう。そこで、この難点を克服するために、１次元連立１階常微分方程式を解くための、例えば、Ｇｅａｒの公式のような硬い方程式専用の、精度の良い、安定した解法を用いて、３次元差分方程式の解を近似する。以下において、この近似法の概要を説明する。
【００１３】ｑ位のＧｅａｒの公式は、多段階後退差分法を用いて、数２の形で与えられる。
【００１４】
【数２】
2005004709.tif 000004
2005004709.tif 000005
す。またｈはメッシュ幅である。
【００１６】数１のｘ，ｙ，ｚの各方向にＧｅａｒの公式（数２）を適用すると、数３が得られる。
【００１７】
【数３】
2005004709.tif 000006
2005004709.tif 000007
【００１８】ここでｈｑは、解析領域をｘ，ｙ，ｚ方向其々に等間隔ｈに分割したとき、ｈをさらにｑ分割したメッシュ幅である。ｊ±ｍ／ｑは、ｊとｊ±１の間のｊからｍ番目の格子点を表し、符号±は、格子点ｊに流入する流れが逆流のとき＋、順流のとき−をとる。ｋ，ｌについても同様である。
【００１９】Ｕｉ（ｊ±ｍ／ｑ，ｋ，ｌ），Ｕｉ（ｊ，ｋ±ｍ／ｑ，ｌ），Ｕｉ（ｊ，ｋ，ｌ±ｍ／ｑ）は既知の値であるから、数３は未知変数Ｕｉ（ｊ，ｋ，ｌ）に関する非線形関数となる。そこで、数３を満たす解を数４で定義する。
【００２０】
【数４】
2005004709.tif 000008
【００２１】一方、数１と同じ速度成分及び非線型関数ｆｉを用いて、ｘ，ｙ，ｚの各方向に対して、数５のような１次元連立１階常微分方程式を考える。
【００２２】
【数５】
2005004709.tif 000009
【００２３】ここで、数５の各式にＧｅａｒの公式（数２）を適用すると、数６を得る。但し、記号の意味は、数３と同じである。
【００２４】
【数６】
2005004709.tif 000010
2005004709.tif 000011
【００２５】数６のＵｉ（ｊ，ｋ，ｌ）に関して、ｘ，ｙ，ｚの各方向について解いたときの解を数７で定義する。
【００２６】
【数７】
2005004709.tif 000012
【００２７】ここで、数７の解と格子点（ｊ，ｋ，ｌ）に流入する速度成分Ｖｘ，Ｖｙ，Ｖｚを用いて加重平均、数８を定義する。但し、これらの速度成分はすべて正にとる。
【００２８】
【数８】
2005004709.tif 000013
【００２９】ここで、符号±は、数３及び数６と同じ意味で、格子点（ｊ，ｋ，ｌ）に流入する速度成分が、順流ならば−、逆流ならば＋をとることを示す。従って、数４は、数３，６，７，８を用いると数９のような形で表される。
【００３０】
【数９】
2005004709.tif 000014
【００３１】ここで、Ｖｓｕｍは数８の分母を表し、Ｆｉ（ｊ，ｋ，ｌ）は数１０で与えられる。
【００３２】
【数１０】
2005004709.tif 000015
【００３３】ここで、数１０の右辺の第１項は数１１のようにテイラー展開できる。
【００３４】
【数１１】
2005004709.tif 000016
【００３５】数１０と数１１から数１２が得られる。
【００３６】
【数１２】
Ｆｉ（ｊ，ｋ，ｌ）＝Ａ＋Ｃ＋Ｂｈ＋Ｏ（ｈ２）
【００３７】従って、数１２を数９に代入すると、ｈｑ＝ｈ／ｑであるから最終的には、数１３のような式が得られる。
【００３８】
【数１３】
2005004709.tif 000017
【００３９】ここで、Ｏ（ｈ）は誤差評価におけるｈのオーダーを表す。またＡ，Ｂ，Ｃは、すべて格子点（ｊ，ｋ，ｌ）よりも過去の点、即ち計算の進行方向に対して上流側で計算される値である。
【００４０】数１３は、数１を３次元後退差分によって離散化したときの解は、１次元解を合成することによって、ｈのオーダーで近似できることを示している。
【００４１】次に、逆流がある場合の速度場では、全格子点での濃度は、有限回の反復計算で近似値できること説明する。いま、解析領域内の、逆流がある任意の格子点をＰ１とすると、Ｐ１に流入する速度成分には、計算の進行方向とは逆向きの成分が少なくとも１つが存在することになる。この逆流成分を逆向きに辿っていくと、仮定１により、必ず入口境界に辿り着く。どの辿り方を選んでも同じことがいえる。何故なら、もしそうでないとすると、途中で湧き出しが存在することになる。
【００４２】従って、Ｐ１から入口境界の点ＰＢまでの格子点の道筋（Ｐ１，…，Ｐｎ−１，Ｐｎ，…，ＰＢ）には、Ｐｎ−１に流入するＰｎからの逆流成分があり、Ｐｎに流入するすべての道筋の流れが順流であるような隣合う２点Ｐｎ−１、Ｐｎが存在することは明らかである。Ｐｎ−１に流入する逆流成分方向の濃度は、Ｐｎにおける濃度の近似値とＰｎからＰｎ−１に流入する速度成分から計算される。しかし、計算の進行方向はＰｎ−１→Ｐｎであるから、Ｐｎ−１に流入する逆流成分方向の濃度は、Ｐｎでの仮の値を用いた計算値となる。ところが、計算処理をＰｎに移すと、Ｐｎに流入する流れは、どの道筋をとってもすべて順流なので、Ｐｎでの濃度の近似値は、入口境界値から順次計算される近似値だけを用いて計算される。このことは、全格子点での計算を２回反復すれば、Ｐｎ−１に流入する逆流成分方向の濃度は、仮の値を用いないで、すべて順次計算される過去の近似値だけで計算できることを示している。また、Ｐｎ−１に流入する他の逆流成分がある場合ついては、同様の議論を順次適用すればよい。ところで、このようなＰｎ−１につながる逆流成分をもつ点は有限個しかないことは明らかである。従って、Ｐｎ−１における濃度は、有限回の反復計算を行えば、入口の境界値から順次計算される過去の近似値だけを用いて計算できることになる。
【００４３】上述の議論は、順次Ｐ１にまで適用することができる。ここで、Ｐ１は逆流成分をもつ任意の格子点であったから、全格子点での濃度は、高々有限回の反復計算で近似できることになる。
【００４４】以下、この発明の実施の形態を図を参照して説明する。図１は本発明の硬い系における３次元解の合成解による近似計算法の概要を示すフローチャートである。図中、ステップ１１においては、数１の右辺ｆｉを表し、これは化学種の生成消滅の項と外力項の和である。ステップ１１２では、右手座標系を設定する。但し、座標原点は入口境界面の４つの頂点のうちの任意の１つとし、またｚ方向の向きは、入口境界面から見て出口側の方向を正にとる。ステップ１２では、３次元速度データ及び数１の数値計算に必要な諸データを設定する。ステップ１３は、解析領域において流れが入口から出口へ向うときの入口境界での境界値と、流れが逆流を含むような場合に必要な反復計算の出発値及び計算条件データ等を設定する。
【００４５】ステップ１４では、計算格子点Ｐ（ｊ，ｋ，ｌ）で化学種ｉの濃度Ｕｉを計算するためにＰ（ｊ，ｋ，ｌ）に流入する近傍点と流入速度成分を決める。これについては、図２で詳細に説明する。
【００４６】ステップ１５では、数５の１次元連立１階常微分方程式を解く硬い方程式専用の、例えば、Ｇｅａｒの公式を利用した解法（数６を解く既存のプログラム）を用いて、ｘ，ｙ，ｚの各方向について流れ成分があれば１次元解、即ち濃度を計算する。但し、１次元解を計算するとき、流入速度はすべて正の値を用いる。
【００４７】ステップ１６は、ステップ１４で決められたＰ（ｊ，ｋ，ｌ）に流入する速度成分とステップ１５で計算した各方向からの濃度を数８のように加重平均する。この加重平均をもって、計算格子点Ｐ（ｊ，ｋ，ｌ）における化学種ｉの濃度Ｕｉの近似値とする。
【００４８】ステップ１７では、解析領域の全格子点（入口境界は除く）について、ステップ１４〜１６までの一連の処理を繰り返す（ステップ２０）。
【００４９】ステップ１８では、計算された全格子点での濃度の近似値に対して、前回の値と今回の値を比較し、同じでなければ、前回の値を今回の値で更新して反復計算を行う（ステップ２１）。全格子点での濃度の近似値が前回の値と同じになれば、次の処理へ進む（ステップ１９）。
【００５０】図２〜４は、流入近傍点と流入速度成分を決めるときの手順を示している。但し、添字が増加する方向を計算の進行方向としている。図２の３１は、格子点Ｐ（ｊ，ｋ，ｌ）のｘ方向の速度成分Ｖが正の場合、流れがＰ（ｊ，ｋ，ｌ）に流入するのは点Ｐ（ｊ−１，ｋ，ｌ）からだけであることを示している。従って、点Ｐ（ｊ−１，ｋ，ｌ）の速度成分Ｖｂを調べて、もしＶｂが正ならば、流入近傍点は点Ｐ（ｊ−１，ｋ，ｌ）で流入速度成分はＶｂとなる。ここでＶａが負のときは、２点Ｐ（ｊ＋１，ｋ，ｌ）とＰ（ｊ，ｋ，ｌ）の間では、仮定２により、ｘ方向に関する流れは無いものとして取り扱う。
【００５１】図３の３２では、図２の３１において、Ｖを負にしてＶｂとＶａを入れ替えた場合に相当する。このときの流れの流入近傍点と流入速度成分は其々Ｐ（ｊ＋１，ｋ，ｌ）とＶａになる。ここでＶｂが正のときは、２点Ｐ（ｊ，ｋ，ｌ）とＰ（ｊ−１，ｋ，ｌ）の間では、図２の３１と同様、ｘ方向に関する流れは無いものとして取り扱う。
【００５２】図４の３３は、Ｖが０で、ＶａとＶｂの向きが反対の場合である。ここでＶａｂを数１４のように定義する。
【００５３】
【数１４】
Ｖａｂ≡Ｖａ＋Ｖｂ
【００５４】このとき、Ｖａｂの値によってＰ（ｊ，ｋ，ｌ）への流れは、図４のａ）〜ｃ）のように３つのケースに分けられる。図４の３４では、Ｖａｂが正の場合で、Ｐ（ｊ，ｋ，ｌ）に流入するのは、Ｖａｂの相対速度をもつＰ（ｊ−１，ｋ，ｌ）からだけで、Ｐ（ｊ＋１，ｋ，ｌ）からの流れは無いものとして取り扱う。
【００５５】図４の３５はＶａｂが負の場合で、３４と全く逆の流れの場合に相当する。このときの流入近傍点と流入速度成分は其々Ｐ（ｊ＋１，ｋ，ｌ）とＶａｂであり、Ｐ（ｊ−１，ｋ，ｌ）からの流れはない。
【００５６】図４の３６は、ＶａとＶｂの向きが反対で大きさが同じ場合である。このとき、Ｐ（ｊ，ｋ，ｌ）では、流れが平衡状態と考え、Ｐ（ｊ−１，ｋ，ｌ）とＰ（ｊ＋１，ｋ，ｌ）からの流れは無いものとして取り扱う。
【００５７】以上、図２〜４において、ｘ方向についての流入近傍点と流入速度成分の特定手順を示してきたが、ｙ方向とｚ方向についても全く同様の手順が成り立つ。
【００５８】
【発明の効果】化学反応を伴う速度場の濃度解析では、濃度を３次元差分方程式から直接計算することは、膨大な計算量、計算の安定性、更に反復計算の収束性等の問題が生じるため、非常に難しい。しかしながら、本発明によれば、３次元差分解は、既存の硬い方程式専用の１次元解法で求めた解を、流入する速度成分で加重平均した合成解によって、近似できるという効果がある。しかも、計算は１次元計算と加重平均だけなので非常に簡単であるという利点がある。
【００５９】更に、流れに逆流があるような場合、有限回の反復計算だけで近似解が得られるという利点がある。
【図面の簡単な説明】
【図１】本発明の１実施例による硬い系における３次元解の合成解による近似計算法のフローチャート。
【図２】計算格子点でのｘ方向の速度成分が正のときの流入近傍点と流入速度成分の特定例を示す図である。
【図３】計算格子点でのｘ方向の速度成分が負のときの流入近傍点と流入速度成分の特定例を示す図である。
【図４】計算格子点でのｘ方向の速度成分が零のときの流入近傍点と流入速度成分の特定例を示す図である。
【図１】
2005004709.tif 000018
【図２】
2005004709.tif 000019
【図３】
2005004709.tif 000020
【図４】
2005004709.tif 000021