(57)【要約】【課題】非分岐かつ非直交の構造格子を用いて反復計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、数値計算が不安定とならず、精度の高い解が得られる様に設定する。【解決手段】境界に垂直な方向の微分係数が０でない箇所の計算セルの境界条件を、境界を軸として前記計算セルと対称位置にあるセルをｖ計算領域外に仮想し、仮想した計算セルの反変速度をξη座標系とデカルト座標系との座標変換係数を用いて定義し、反変速度を用いて与えようとする対称境界条件を定義し、運動方程式を反復計算法により解く際に、計算領域の座標変換係数を用いて計算領域外に仮想した計算セルについての座標変換係数を表現し、前回の反復計算で境界条件を反映して得られたデカルト座標の物理速度と境界条件を用いて今回の反復計算における未知数を算出する。【選択図】図２
2011065360.tif 000002
【特許請求の範囲】
【請求項１】
非分岐かつ非直交の構造格子（ξη座標系）を用い、初期条件の設定後反復計算法を用いて流れの数値解析を行う際に、境界に垂直な方向の微分係数が０でない箇所の計算セルの境界条件を設定し、反復計算を用いる流れの数値解析方法であって、
境界に垂直な方向の微分係数が０でない箇所の計算セルの境界条件を、境界を軸として前記計算セルと対称位置にあるセルを計算領域外に仮想し、該仮想した計算セルの反変速度を前記ξη座標系とデカルト座標系との座標変換係数を用いて定義し、さらに定義した反変速度を用いて与えようとする対称境界条件を定義し、
運動方程式を反復計算法により解く際に、計算領域の座標変換係数を用いて計算領域外に仮想した計算セルについての座標変換係数を表現する様にし、前回の反復計算で境界条件を反映して得られたデカルト座標（ｘｙ座標系）の物理速度と境界条件を用いて今回の反復計算におけるデカルト座標の物理速度未知数を算出すること特徴とする非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法。
【請求項２】
前記非分岐かつ非直交の構造格子（ξη系座標）とデカルト座標（ｘｙ系座標）との変換係数をξｘ、ξｙ、ηｘ、ηｙとし、デカルト座標におけるｘ方向をｕ、ｙ方向の速度をｖとして、前記反変速度ＵとＶを、以下の（式６）と（式７）で定義し、前記境界に垂直な方向の微分係数が０でない箇所の計算セルにおける計算セルのＵについての境界（ξ＝Ｎξ）条件
【数１】
2011065360.tif 000049
を、当該計算セルと対象の位置にあると仮想した計算セルのＵを用いて、（式８）で定義し、前記仮想した計算セルに対して以下の（式９）で示される座標変換係数を用い、
これらにより前記仮想した計算セルのｘ方向の速度とｙ方向の速度を、前回の繰返し計算で得た値を用いて解くことを繰り返すことを可能としていること特徴とする請求項１に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法。
【数２】
2011065360.tif 000050
【数３】
2011065360.tif 000051
【数４】
2011065360.tif 000052
【数５】
2011065360.tif 000053
【請求項３】
前記仮想した計算セルのｘ方向の速度とｙ方向の速度を、前回の繰返し計算で得た値を用いて解くことを繰り返すとは、
（式９）に示す関係を（式８）に代入して、
【数６】
2011065360.tif 000054
について解くことにより以下の（式１０）の結果を得、
同じく（式７）に示す関係を（式８）に代入して、
【数７】
2011065360.tif 000055
について解くことにより（式１１）の結果を得、
（式１０）および（式１１）を反復しながら対称境界条件を実現するものであることを特徴とする請求項２に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法。
【数８】
2011065360.tif 000056
【数９】
2011065360.tif 000057
【請求項４】
前記反復計算法は、
コロケート格子法であることを特徴とする請求項２または請求項３に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法の方法に関する。
【背景技術】
【０００２】
複雑な形状の領域（物体の周り）の流れの数値解析は、解析対象の領域をメッシュ分割し、与えられた境界条件、初期値等を基に繰返し計算を行うことによりなされるが、境界条件の１としての流路形状が複雑あるいは大きく歪んでいる場合、分割をなす格子線が碁盤の目の様に直交する座標系（デカルト座標系）では、計算格子を流路の境界（以下、単に「境界」と記す）に沿って上手く設定することが困難である。
【０００３】
このため、例えば、当該領域の外形や解析対象物の外形の内で重要あるいは主な線や面に沿っている等の座標系、一般曲線座標系といわれている座標系の構造格子（以下、単に「構造格子」と記す）を用いて解析することがなされている。図１に、かかる構造格子の一例を示す。図１において、１１は下側の境界であり、１２は右側の境界であり、１３は上側の境界であり、１４は左側の境界である。
【０００４】
そして、構造格子はξ軸とη軸で示す。これら２つの軸は共に数値解析の為に導入した構造格子用の軸であり、各々図１の左下側に示す直交座標（デカルト座標）系の左右（水平）方向のｘ軸、上下（垂直）方向のｙ軸に替えて用いる。そして、各格子線はξ座標値あるいはη座標値の何れか一方が同じ値となる線であり、さらにξ座標およびη座標で測った各升目の大きさ（対向する辺の座標値の差）は等間隔としてある。
【０００５】
図１においては、下側の境界１１はｘ軸に沿った直線状であり、右側の境界１２はｙ軸に沿った直線状であり、このためこれら２つの境界は直交しているが、上側の境界１３は直線状ではあるがｘ軸とｙ軸に対して傾斜しており、左側の境界１４は曲線状である。このため、ξ軸系では、格子線は単に傾斜しているだけでなく、η座標値によってｘ軸、ｙ軸となす傾斜角が相違している。また、η軸系では、格子線は曲線状になっているだけでなく、ξ座標値が０から大きくなるにつれて、即ち左側から右側に行くにつれて、左側の境界１４に沿った曲線から右側の境界１２に沿った直線に近くなっており、右側端では直線となっている。
【０００６】
図１の構造格子は不等間隔格子と呼ばれるものであり、前記したように、ξ軸上のη軸との交点のξ座標値は等間隔だが、ｘ座標およびｙ座標で表した場合、ξ軸に沿って徐々に大きくなる。このようなξ軸上のη軸との交点のｘ、ｙ座標は、ξ＝１、２、・・・Ｎξ、η＝１、２、・・・Ｎηとして、以下の（式１）および（式２）により求められる。
【０００７】
【数１】
2011065360.tif 000003
【０００８】
【数２】
2011065360.tif 000004
【０００９】
上記（式１）および（式２）において、Ｒは、図１の左側境界１４を現す円弧の半径であり、中かっこ内第一項により、曲線状の左側境界１４から直線状の右側境界１２に滑らかに変化するη軸系格子線が表現される。また、θ（η）とχ（ξ）は、以下の（式３）および（式４）により定義される。
【００１０】
【数３】
2011065360.tif 000005
【００１１】
【数４】
2011065360.tif 000006
【００１２】
なお、（式４）におけるΔは、不等間隔格子の配置を調節するパラメータであり、Δを大きくすると最大格子間隔と最小格子間隔の比が大きくなる。図１ではΔ＝０．８程度としている。
【００１３】
数値計算の内容は、計算対象の領域に対して構造格子を設定し、与えられた境界条件、初期値の下で、各計算セルの内部で、あるいは各計算セルとそれに隣接する計算式との間で、連続の式とナビヤ・ストークスの式が成立するものとして計算式をたて、離散化し、近似式を繰返し計算で解く、即ち数値解を得るものである。なお、離散化の方法としては、有限差分法、有限体積法等があり、近似式の作成方法としてはテーラー展開、多項式近似、コンパクト差分法等がある。
【００１４】
具体的な計算の内容としては、図１の全てのセルにおいて、速度ｕ、ｖおよび圧力ｐの初期値を仮定し、これらが連続の式とナビヤ・ストークスの式、および後述する境界条件を満たすように、速度ｕ、ｖおよび圧力ｐの順に、セル（ξ＝１，η＝１）からセル（ξ＝Ｎξ，η＝Ｎη）までの全てのセルにおける値が順次更新されることがなされる。
【００１５】
ここでは、有限体積法に基づく離散化がなされ、あるセルにおける速度ｕ、ｖおよび圧力ｐの値を、当該セルおよび隣接セルにおける速度ｕ、ｖおよび圧力ｐによって表す関係式が得られる。
【００１６】
この関係式は、速度ｕ、ｖおよび圧力ｐの各々に対して得られる。例えば、あるセルにおける速度ｕは、当該セルおよび隣接セルにおける速度ｕ、ｖおよび圧力ｐによって表され、速度ｖおよび圧力ｐについても同様である。
【００１７】
そして、反復計算では、これら各セルにおける速度ｕ、ｖおよび圧力ｐの値が、当該セルおよび隣接セルにおける速度ｕ、ｖおよび圧力ｐの前回反復値を用いて更新され、前回反復値と更新値との差が許容し得るまでに小さくなった状態に達すれば、全てのセルにおける速度ｕ、ｖおよび圧力ｐの値が、連続の式とナビヤ・ストークスの式、および境界条件を満たしたと判断し、反復計算を終了する（特許文献１、非特許文献１）。
【先行技術文献】
【特許文献】
【００１８】
【特許文献１】特開２０００−２８５１０２号公報
【非特許文献】
【００１９】
【非特許文献１】荒川忠一著，「数値流体工学」，東京大学出版会，１９９４年１月
【発明の概要】
【発明が解決しようとする課題】
【００２０】
しかしながら、流路形状が複雑あるいは大きく歪んでいる場合には、構造格子では計算に用いる格子線は境界線（壁）等と必ずしも直交しなくなるため、境界条件を正確には表現できない場合がある。この様な非直交格子（同じ座標値を示す座標線が境界壁に直交しない格子）が問題になるのは、境界線に垂直な方向の微分係数を与えるタイプの境界条件の場合であり、計算格子上に配置される速度ｕ、ｖおよび圧力ｐの値を用いてこの微分係数を正確に表現することができない。
【００２１】
なお、この微分係数が０である場合には、これらの値が境界線に対して傾きをもたず、境界付近で平坦な分布であることを示し，対称境界条件が表現されたことを表す。
【００２２】
この例を、図１と図２を用いて説明する。図２は、図１の右上の境界部分に在る４個のセルと、当該４個のセルのη軸に線対称な位置に仮想した４個のセルを示す図である。図１と図２に示す様に、数値計算の対象となる図形の右側境界と上側境界においては、ξ座標とη座標で示されるセルの辺（格子線）と境界が直交していない。特に図１では右上側の境界ほどその様になっている。
【００２３】
対称境界条件を、例えば速度ｕ、ｖおよび圧力ｐ等の物理量φに与える場合、一般には次の（式５）の様になる。
【００２４】
【数５】
2011065360.tif 000007
【００２５】
なお、（式５）において、
【００２６】
【数６】
2011065360.tif 000008
は、Ｎξ番の計算セルの中心に配置される物理量φを示す。
【００２７】
（式５）は、ξ軸方向（あるいはη軸方向）の微分係数が０という境界条件、即ちξ軸（あるいはη軸）が境界と直交する場合には対称という条件を与えるが、図２に示す様に直交しない場合には微分係数が０、そして領域の境界壁では流入、流出がないため境界線に関して対称であるという条件を表していない。なお、これはξ軸方向が境界線に直交しない場合であるが、η軸方向についても同様である。
【００２８】
その結果、境界で見かけの流入、流出が生じることとなり、質量等の保存則が成立せず、前回反復値と更新値との差が反復に伴って拡大したり、あるいはその差が増減を繰返すことにより十分に小さくならない、いわゆる数値計算の不安定な状態となり、数値解を得られなくなることがある。
【００２９】
図３に、図１に示す領域を対象として、境界で見かけの流入、流出が生じていることとなった計算例を示す。これは、境界条件として右側境界と上側境界においては対称境界条件であるため流れの出入りが無いはずの問題を、コロケート格子法を用い、離散化し、繰返し計算を行って一応の解を求めた途中結果である。即ち、この計算は完全に収束させることはできず、これ以上繰返し計算を行うと解が発散してしまう為、発散する前の結果を示したものである。
【００３０】
そこで、図４のように、ξ軸方向あるいはη軸方向の格子線を境界線に対して直交させるよう無理な修正を行なうと、計算格子が分岐しないという構造格子の条件により、直交しない右側の計算セルが非常に大きくなり、計算領域内での計算セルの最大／最小面積比が大きくなり、数値計算の不安定な状態が生じ、計算精度が低下する要因となる。
【００３１】
その対策として、図５のように、分岐を許容した非構造格子を用いて計算セルの最大／最小面積比を小さくすることもできるが、プログラムが複雑となり、計算に必要なメモリも増大し、計算速度も遅くなる。
【００３２】
このため、非分岐かつ非直交の構造格子を用いて計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、数値計算が不安定とならず、精度の高い解が得られる様に設定し、反復計算を用いる流れの数値解析方法の開発が望まれていた。
【課題を解決するための手段】
【００３３】
本発明は、以上の課題を解決することを目的としてなされたものであり、非分岐かつ非直交の構造格子を用いて計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、境界を軸として対称位置にあるセルを計算領域外に仮想し、反変速度を用いて境界条件を定義し、反復計算法を用いることにより１の境界条件から２つの未知数を算出可能とし、またこのために計算領域内の座標変換係数を用いて、計算領域外に仮想した計算セルについての座標変換係数を表現する様にして、反復計算を用いる流れの数値解析方法としたものである。以下、各請求項の発明を説明する。
【００３４】
請求項１に記載の発明は、
非分岐かつ非直交の構造格子（ξη座標系）を用い、初期条件の設定後反復計算法を用いて流れの数値解析を行う際に、境界に垂直な方向の微分係数が０でない箇所の計算セルの境界条件を設定し、反復計算を用いる流れの数値解析方法であって、
境界に垂直な方向の微分係数が０でない箇所の計算セルの境界条件を、境界を軸として前記計算セルと対称位置にあるセルを計算領域外に仮想し、該仮想した計算セルの反変速度を前記ξη座標系とデカルト座標系との座標変換係数を用いて定義し、さらに定義した反変速度を用いて与えようとする対称境界条件を定義し、
運動方程式を反復計算法により解く際に、計算領域の座標変換係数を用いて計算領域外に仮想した計算セルについての座標変換係数を表現する様にし、前回の反復計算で境界条件を反映して得られたデカルト座標（ｘｙ座標系）の物理速度と境界条件を用いて今回の反復計算におけるデカルト座標の物理速度未知数を算出すること特徴とする非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法である。
【００３５】
本請求項の発明においては、非分岐かつ非直交の構造格子を用いて計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、境界を軸として対称位置にあるセルを計算領域外に仮想し、反変速度を用いて境界条件を定義し、反復計算法を用いることにより１の境界条件から２つの未知数が算出可能となるため、数値計算が不安定とならず、精度の高い解が得られる。
【００３６】
なおここに、「非分岐」とは、升目を構成する格子線が途中で途切れないことを指す。
また、「非直交」とは、格子線と境界とが直交しないことを指す。
また、「構造格子」とは、一般曲線座標系を指す。
また、「座標変換係数」とは、デカルト座標と一般曲線座標との変換係数を指す。
また、「反変速度」とは、一般曲線座標の格子線に垂直な方向の速度を指す。
【００３７】
請求項２に記載の発明は、
前記非分岐かつ非直交の構造格子（ξη系座標）とデカルト座標（ｘｙ系座標）との変換係数をξｘ、ξｙ、ηｘ、ηｙとし、デカルト座標におけるｘ方向をｕ、ｙ方向の速度をｖとして、前記反変速度ＵとＶを、以下の（式６）と（式７）で定義し、前記境界に垂直な方向の微分係数が０でない箇所の計算セルにおける計算セルのＵについての境界（ξ＝Ｎξ）条件
【００３８】
【数７】
2011065360.tif 000009
を、当該計算セルと対象の位置にあると仮想した計算セルのＵを用いて、（式８）で定義し、前記仮想した計算セルに対して以下の（式９）で示される座標変換係数を用い、
これらにより前記仮想した計算セルのｘ方向の速度とｙ方向の速度を、前回の繰返し計算で得た値を用いて解くことを繰り返すことを可能としていること特徴とする請求項１に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法である。
【００３９】
【数８】
2011065360.tif 000010
【００４０】
【数９】
2011065360.tif 000011
【００４１】
【数１０】
2011065360.tif 000012
【００４２】
【数１１】
2011065360.tif 000013
【００４３】
本請求項の発明においては、非分岐かつ非直交、２次元の構造格子（ξη系座標）を用いて計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、境界を軸として対称位置にあるセルを仮想し、（式２）と（式３）で定義される反変速度を用いて境界条件を定義し、反復計算法を用いることにより１の境界条件からデカルト座標における２つの未知数を算出可能とすることにより、数値計算が不安定とならず、精度の高い解が得られる。
【００４４】
請求項３に記載の発明は、
前記仮想した計算セルのｘ方向の速度とｙ方向の速度を、前回の繰返し計算で得た値を用いて解くことを繰り返すとは、
（式９）に示す関係を（式８）に代入して、
【００４５】
【数１２】
2011065360.tif 000014
について解くことにより以下の（式１０）の結果を得、
同じく（式７）に示す関係を（式８）に代入して、
【００４６】
【数１３】
2011065360.tif 000015
について解くことにより（式１１）の結果を得、
（式１０）および（式１１）を反復しながら対称境界条件を実現するものであることを特徴とする請求項２に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法である。
【００４７】
【数１４】
2011065360.tif 000016
【００４８】
【数１５】
2011065360.tif 000017
【００４９】
本請求項の発明においては、以下の（１）から（３）の要点に基づいているため、請求項２の発明の効果がより良く発揮できるだけでなく、プログラム作成が容易となり、計算メモリや計算速度の点でも有利となる。
【００５０】
（１）反変速度が座標軸に垂直な方向を有することに着目して、非対称境界条件を、前記（式８）で表す。
（２）前回反復値に既に境界条件が反映されていることに着目して、（式８）に含まれる２つの未知数を求めるために反復計算法とする。
（３）境界に関する対称性を利用して、計算領域外の座標変換係数を、計算領域内の座標変換係数で表す。
【００５１】
請求項４に記載の発明は、
前記反復計算法は、
コロケート格子法であることを特徴とする請求項２または請求項３に記載の非分岐かつ非直交の構造格子の境界条件を設定し、反復計算を用いる流れの数値解析方法である。
【００５２】
本請求項の発明においては、コロケート格子法で反復計算を行うため、請求項２または請求項３の発明を流れの数値解析に用いる際に便利である。
【発明の効果】
【００５３】
本発明は、非分岐かつ非直交の構造格子を用いて計算を行う際に、境界に垂直な方向の微分係数が０でない場合の境界条件を、境界を軸として対称位置にあるセルを仮想し、反変速度を用いて境界条件を定義し、反復計算法を用いることにより１の境界条件から２つの未知数を算出可能としているため、数値計算が不安定とならず、精度の高い解が得られる。
【図面の簡単な説明】
【００５４】
【図１】非直交の構造格子の１例を示す図である。
【図２】図１の右上の境界部分と当該部分のη軸に線対称な仮想上の部分を示す図である。
【図３】図１に示す領域を対象として、境界で見かけの流入、流出が生じていることとなった計算例を示す図である。
【図４】直交構造格子の１例を示す図である。
【図５】直交非構造格子の１例を示す図である。
【図６】本発明を用いた数値計算の結果を示す図である。
【発明を実施するための形態】
【００５５】
以下、本発明をその実施の形態に基づいて説明する。なお、本発明は、以下の実施の形態に限定されるものではない。本発明と同一および均等の範囲内において、以下の実施の形態に対して種々の変更を加えることが可能である。
【００５６】
図１に示す領域を対象に、同じ構造格子を用いて、即ち同じ格子分割、方法、計算式に与えるのではなく境界壁が本来的に有する境界条件で流れの数値解析を行うものとする。
【００５７】
以下、計算方法について、概略説明する。
計算方法としては、本実施の形態ではコロケート格子法を用いる。この方法は、物理速度ｕとｖは計算セル中心に、反変速度（格子線に垂直な方向の速度）ＵとＶは計算セルの辺の中心に配置するものである。
【００５８】
ｕをｘ軸方向の物理速度、ｖをｙ軸方向の物理速度、Ｕはη軸と垂直な向き（直交する）の反変速度であり、Ｖをξ軸と垂直な向き（直交する）の反変速度とする。反変速度は、以下の（式６）、（式７）にて定義される。なお（式６）、（式７）において、ξｘ、ξｙ、ηｘ、ηｙは、ξ＝ξ（ｘ，ｙ）、η＝η（ｘ，ｙ）とｘ＝ｘ（η，ξ）、ｙ＝ｙ（η，ξ）における座標変換係数である。
【００５９】
【数１６】
2011065360.tif 000018
【００６０】
【数１７】
2011065360.tif 000019
【００６１】
これを用いて、右側境界（ξ＝Ｎξ）における
【００６２】
【数１８】
2011065360.tif 000020
で示す境界条件を、Ｎξ番とＮξ＋１番の計算セル中心に配置される物理速度ｕ及びｖで表すと、次の（式８）の様になる。
【００６３】
【数１９】
2011065360.tif 000021
【００６４】
ここに、
【００６５】
【数２０】
2011065360.tif 000022
とは反変速度が０であることを、即ち格子線に垂直な方向の速度が０であり、流入、流出がないことを指す。
【００６６】
さらに、（式８）を
【００６７】
【数２１】
2011065360.tif 000023
について解くと（式１０）を得る。
【００６８】
【数２２】
2011065360.tif 000024
【００６９】
また、（式８）を
【００７０】
【数２３】
2011065360.tif 000025
について解くと（式１１）を得る。
【００７１】
【数２４】
2011065360.tif 000026
【００７２】
なお、（式１０）および（式１１）における上添え字ｍは、反復回数を示している。また、（式８）から（式１１）における下添え字Ｎξは右側の境界の１つ左のξ座標値であり、Ｎξ＋１は右側の境界線を基準にしてその対称位置におけるξ座標値を表す。そして、鍵かっこ「［    ］」は、複数の下添え字の重複を避けるために用いている。
【００７３】
（式８）の
【００７４】
【数２５】
2011065360.tif 000027
と
【００７５】
【数２６】
2011065360.tif 000028
についての具体的な解き方は、以下の点に着目して行う。
【００７６】
第１に、反復計算を用いることにより１つの境界条件から２つの未知数を算出する。即ち、（式８）には未知数が２つあるため、前回の反復値を用いて（式１０）、（式１１）を得る。
【００７７】
【数２７】
2011065360.tif 000029
【００７８】
【数２８】
2011065360.tif 000030
【００７９】
第２に、計算領域外の座標変換係数を、（式９）に示す様に計算領域内の座標変換係数を用いて算出する。
【００８０】
【数２９】
2011065360.tif 000031
【００８１】
具体的な計算の内容として、図１の全てのセルにおいて速度ｕ、ｖおよび圧力ｐの初期値を仮定し、これらが連続の式とナビヤ・ストークスの式、および後述する境界条件を満たすように、速度ｕ、ｖおよび圧力ｐの順に、セルξ＝１、η＝１からセルξ＝Ｎξ、η＝Ｎηまでの全てのセルにおける値が順次更新されることがなされる。ここでは、有限体積法に基づく離散化がなされ、あるセルにおける速度ｕ、ｖおよび圧力ｐの値が、当該セルおよび隣接セルにおける速度ｕ、ｖおよび圧力ｐによって表される関係式が得られる。
【００８２】
この関係式は、セルＰにおける速度ｕ、ｖおよび圧力ｐをφＰ、それに隣接する４個のセルＥ、Ｗ、Ｎ、Ｓのそれらの値をφＥ、φＷ、φＮ、φＳとして以下の（式１２）で表される。
【００８３】
【数３０】
2011065360.tif 000032
【００８４】
ここで、ａＥ、ａＷ、ａＮ、ａＳは、各々（式１３）から（式１６）で示される。
【００８５】
【数３１】
2011065360.tif 000033
【００８６】
【数３２】
2011065360.tif 000034
【００８７】
【数３３】
2011065360.tif 000035
【００８８】
【数３４】
2011065360.tif 000036
【００８９】
そして、ｂφは格子の非直交性に関係する項である。また、Ｊは計算セルにおける座標変換係数によって計算されるヤコビアンであり、Ｓは流れを駆動するソース項である。
【００９０】
反復計算では、前回反復時のφＥ、φＷ、φＮ、φＳ、および前回反復時のφＰを用いて評価したＳを用いて右辺を計算し、φＰを新しい値に更新する。φＰの前回反復値と更新値との差が許容し得るまでに小さくなった状態に達すれば、全てのセルにおける速度ｕ、ｖおよび圧力ｐの値が、連続の式とナビヤ・ストークスの式、および境界条件を満たしたと判断し、反復計算を終了する。
【００９１】
このφＰの値は、セルξ＝１、η＝１からセルξ＝Ｎξ、η＝Ｎηまでの全てのセルに対して評価されるので、例えば、セルξ＝Ｎξあるいはη＝Ｎηのような境界点において、セルξ＝Ｎξ＋１あるいはη＝Ｎη＋１のような境界外点における値が必要になる。
【００９２】
本発明は、境界線において反変速度の１成分（ＵあるいはＶ、境界ξ＝Ｎξでは前者、η＝Ｎη＋１では後者となる）に対して成立する境界条件に基づく反復解法を用いることにより、これら各境界外点において必要となる速度の２成分（ｕおよびｖ）を与えることを可能とするものである。
【００９３】
すなわち、境界ξ＝Ｎξを例にすると、対称境界条件は、境界線に垂直な速度成分、すなわち反変速度Ｕを用いて厳密にＵ＝０と表せる。これを離散式で表すと（式８）を得る。
【００９４】
【数３５】
2011065360.tif 000037
【００９５】
この（式８）をｕについて解くと（式１０）を得る。
【００９６】
【数３６】
2011065360.tif 000038
【００９７】
また、ｖについて解くと（式１１）を得る。
【００９８】
【数３７】
2011065360.tif 000039
【００９９】
これらの式の右辺におけるｕ、ｖとして、前回反復値を用いて、上記の反復計算を行なえば、対称境界条件を満たすよう、速度ｕ、ｖが更新されることになる。
【０１００】
ところで、上式には計算領域外の座標変換係数ξｘ、ξｙ、ηｘ、ηｙが含まれる。これらを、計算領域内の座標変換係数を用いて算出することにより，上記反復計算を計算領域内で実施することを可能にした。
【０１０１】
（式１０）、（式１１）に含まれる座標変換係数は，座標が既知である計算領域内の計算セルに対し計算可能である。そこで，第２図に示すように，非直交な境界付近において、与えようとしている対称境界条件を満たす計算セルを仮想的に配置することにより、座標の分からない計算領域外の座標変換係数を，計算領域内の座標変換係数により表現する。すなわち、以下の（式９）に示す関係を用いる。
【０１０２】
【数３８】
2011065360.tif 000040
【０１０３】
（式１０）、（式１１）は，この関係を（式８）に代入した結果得られたものである。
【０１０４】
一方、圧力ｐは、速度ｕ、ｖが連続の式を満たすように決める必要があるが、その方法としてＳＩＭＰＬＥ法がある。ＳＩＭＰＬＥ法に基づく反復計算を行なえば、速度ｕ、ｖに対して与えられた上記の境界条件を維持しつつ、連続の式を満たすように速度ｕ、ｖが更新される。そして、全てのセルにおける速度ｕ、ｖおよび圧力ｐの値が、連続の式、ナビヤ・ストークスの式および境界条件Ｕ＝０を満たされたと判断し、反復計算を終了する。
【０１０５】
順序は逆になるが、本実施の形態の方法を適用する運動方程式について説明する。図１に示す領域内で、ｕとｖに関する非圧縮性流体の定常流の運動方程式は、（式１７）、（式１８）である。
【０１０６】
【数３９】
2011065360.tif 000041
【０１０７】
【数４０】
2011065360.tif 000042
【０１０８】
なお、（式１７）、（式１８）に含まれるｐは、計算領域内の質量が保存されるよう、ＳＩＭＰＬＥ解法を用いた。
【０１０９】
また、（式１７）、（式１８）における各記号の意味は、以下の（式１９）から（式２２）に示す通りである。また、ρは密度（ｋｇ／ｍ３）であり、νは動粘度（ｍ２／ｓ）であり、ｐは圧力（Ｐａ）である。
【０１１０】
【数４１】
2011065360.tif 000043
【０１１１】
【数４２】
2011065360.tif 000044
【０１１２】
【数４３】
2011065360.tif 000045
【０１１３】
【数４４】
2011065360.tif 000046
【０１１４】
なお、この（式１７）と（式１８）は、等間隔で分割したそして分岐しないデカルト座標系においては、各々（式２３）と（式２４）の様になる。
【０１１５】
【数４５】
2011065360.tif 000047
【０１１６】
【数４６】
2011065360.tif 000048
【０１１７】
図６に、数値計算の結果を示す。図３と異なり、上側及び右側の境界において流れの対称条件が充たされており、流入も流出も生じていない水路内を循環する流れが示されている。
【産業上の利用可能性】
【０１１８】
本発明は、対象とする領域が複雑な形状であるだけでなく、座標線と境界が直交せず境界に垂直方向の微分係数を与えるタイプの計算が多い原子炉等の数値解析に利用される。
【符号の説明】
【０１１９】
１１      下側の境界
１２      右側の境界
１３      上側の境界
１４      左側の境界
【図１】
2011065360.tif 000058
【図２】
2011065360.tif 000059
【図３】
2011065360.tif 000060
【図４】
2011065360.tif 000061
【図５】
2011065360.tif 000062
【図６】
2011065360.tif 000063