構造格子を用いたシミュレーション方法
【課題】流体解析等に好適に使用できる陰解法を用いたシミュレーション解法に関する。
【解決手段】二次元又は三次元空間内に、構造格子のモデルを設定する工程と、この構造格子のモデルに関連付けられた物理量及び条件に基づいてポアソン方程式を設定する工程と、該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、前記計算工程は、予め定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記ポアソン方程式を近似的に解くことを特徴とする。
【解決手段】二次元又は三次元空間内に、構造格子のモデルを設定する工程と、この構造格子のモデルに関連付けられた物理量及び条件に基づいてポアソン方程式を設定する工程と、該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、前記計算工程は、予め定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記ポアソン方程式を近似的に解くことを特徴とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、例えばコンピュータを用いた流体解析等に好適に使用できる構造格子を用いたシミュレーション方法に関する。
【背景技術】
【０００２】
任意の空間を流れる流体を解析する場合、コンピュータを用いたシミュレーションが行われる。このシミュレーションでは、任意の三次元空間を、構造格子を用いてメッシュ化し、各格子点での流体の圧力、温度及び／又は流速等が計算される。そして、この中の圧力に関する方程式は、流体を非圧縮であると定義できる場合、通常、ポアソン方程式を解くことにより行われる。
【０００３】
前記ポアソン方程式の計算方法として、例えばヤコビ法、ガウス・ザイデル法、ＳＯＲ法などの反復法が知られている。反復法は、初期値を与え（手順１）、反復計算を行い（手順２）、収束判定条件を満たす解になるまで手順２を繰り返す（手順３）、ことで行われる。しかしながら、上記方法は、解くべきマトリクスサイズが大きくなると、極端に収束性が悪化し、多くの時間を必要とする。
【０００４】
このため、近年では、反復数がメッシュ規模に依存しないよう、ＡＭＧ（Algebraic Multigrid Methods）法が多用される傾向がある。周知のように、ＡＭＧ法は、直交格子だけでなく非構造格子も解くことができるように、幾何マルチグリッド法を拡張したものである。ＡＭＧ法は、マルチグリッド法と同様、反復数がメッシュサイズに依存しないという特徴があるため、現在、最も高速な連立一次方程式の解法と言え、市販されている汎用ソルバーにも、多く採用されている（例えば、ＡＮＳＹＳ社の「Fluent」、ＣＤ−adapco社の「ＳＴＡＲ−ＣＤ」等）。
【０００５】
しかしながら、ＡＭＧ法は、大規模マトリクスを解くときの収束性に優れる一方、コンピュータに大規模なメモリーサイズを要求する。このため、大きなメモリーを搭載したコンピュータでなければ、ＡＭＧ法で大規模な計算ができないという問題があった。
【０００６】
一方、必要なメモリーサイズが少なくても大規模マトリクスの計算が可能で収束性の良い方法として、直接法の一つであるブロックサイクリックリダクション法が知られている。しかしながら、ブロックサイクリックリダクション法は、ポアソン方程式（例えば請求項１の式１）で、kが不変量(constant value)である場合を前提としており、kが一定値ではなく、マトリクスの中で一カ所でも違う値をとってしまい、マトリクスの外に出せない様なケースでは、ブロックサイクリックリダクション法）を用いて解く事はできない。
【０００７】
従って、ブロックサイクリックリダクション法を正しく適用するためには、シミュレーションで用いられるkの値を一定にすること、つまり、陽解法を用いる（クーラン条件（CFL条件）を守る）必要があり、計算時の時間刻み幅を小さくとる必要があるので、計算コストが大きくなってしまうという欠点がある。
【０００８】
なお、CFL条件（Courant-Friedrichs-Lewy Condition）とは、コンピュータシミュレーションの計算（数値解析）において、「情報が伝播する速さ」を「実際の現象で波や物理量が伝播する速さ」よりも早くしなければならないという必要条件のことである。例えば、離散格子系において波動を扱う場合に、その運動方程式の数値解を求める際に用いる時間ステップ（Δt）の値は、実際の波動が隣り合う格子に伝達するまでの時間よりも小さくなければならない。もし（Δt）の値がCFL条件の上限を超えると、数値発散が生じてしまい、物理的に意味の無い解を得てしまう。また、格子の間隔が小さくなると、時間ステップの上限値も減少する。そして、実際の波(特性速度)の速さをC、計算上の時間ステップ幅をΔt、格子幅をΔxとすれば、 移流方程式等の場合、 Δx/Δt>C (情報が伝播する速さ>実際の波等の速さ)がCFL条件となる。このCFL条件は、陽解法の時間進展を行う際に用いられる条件であり、この条件を回避するために、しばしば陰解法が用いられる場合がある。 陰解法を用いることでCFL条件を回避しかつ緩和できる理由としては様々な説明が存在するが、最も簡潔に説明すると、陽解法は1ステップ前の自分の周りのごくわずかな格子点のみから情報を得て 次の時間の値を決めるのに対して、陰解法は1ステップ前の(ほぼ)全ての格子点の情報を処理して 次の時間の値を決めるためであり、CFL条件におけるΔxが実質的に巨大になるためである。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開平５−１３５０９１号公報
【特許文献２】特開２０００−２９３５０８号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
本発明は、以上のような問題点に鑑み案出なされたもので、二次元又は三次元空間内に、均一又は不均一構造格子を設定し、この構造格子に関連付けられた物理量及び条件に基づいて特定形式のポアソン方程式を設定し、このポアソン方程式を、ブロックサイクリックリダクション法を近似的に計算することを基本として、不均一構造格子を前提とするポアソン方程式についても小さなメモリーサイズのコンピュータで迅速に計算することが可能な不均一構造格子を用いたシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【００１１】
本発明のうち請求項１記載の発明は、二次元又は三次元空間内に、構造格子のモデルを設定する工程と、
この構造格子のモデルに関連付けられた物理量及び条件に基づいて式（１）のポアソン方程式を設定する工程と、
該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、
前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、
この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、
前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする。
【数２】
【００１２】
また請求項２記載の発明は、構造格子は、二次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法である。
【００１３】
また請求項３記載の発明は、前記構造格子は、三次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法である。
【００１４】
また請求項４記載の発明では、前記計算工程は、予め前記式（２）の方程式をフーリエ変換することによって、次元数を二次元に落とすことを特徴とする請求項３記載の構造格子を用いたシミュレーション方法である。
【発明の効果】
【００１５】
本発明によれば、二次元又は三次元空間内に、構造格子で空間を離散化したモデルを設定する工程と、この構造格子のモデルに関連付けられた物理量及び条件に基づいて特定のポアソン方程式を設定する工程と、該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする。
【００１６】
即ち、ブロックサイクリックリダクション法は直接法であり、本来では、１回の計算で誤差のほとんど無い解が得られるが、k≠constant型のポアソン方程式については、ブロックサイクリックリダクション法の適用条件を満たさず、そのような計算はできない。しかし、本発明では、誤差の概念を導入し、この誤差を、ブロックサイクリックリダクション法を用いて反復計算し、補正によって少しずつ小さくすることを特徴の一つとする。つまり、結果的に、k≠constant型のポアソン方程式を、前記ブロックサイクリックリダクション法を利用して少ない反復回数で解くことを可能にしている。
【００１７】
従って、k≠constantを前提とするポアソン方程式についても、ブロックサイクリックリダクション法を適用できるため、小さなメモリーサイズのコンピュータで計算することができ、かつ、陰解法を適用できる様になるため、大きな時間刻み幅で非定常計算を実施できる様になる。このため、本発明方法では、計算コストの削減にもつながる。これにより、例えば複雑な乱流が生じる空気流体解析等を、高い空間解像度、時間解像度で解くことができ、ゴルフボールのディンプル開発等において、効率化やコスト削減に役立つ。また、ブロックサイクリックリダクション法では、反復回数が少ないため、計算結果に桁落ち等による計算誤差が積み上がる可能性も低く、計算精度の向上にも繋がる。
【図面の簡単な説明】
【００１８】
【図１】本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。
【図２】（ａ）は均一構造格子、（ｂ）は不均一構造格子をそれぞれ説明する斜視図である。
【図３】（ａ）は三次元の構造格子の一つのセルの斜視図、（ｂ）はその正面図である。
【図４】（ａ）は本発明を説明するための二次元の構造格子の例であり、（ｂ）は均一構造格子、（ｃ）は不均一構造格子の例を示す。
【図５】構造格子の互いのマトリクス要素の関係を示す。
【図６】図５の構造格子の関係を示すパターンをマトリクスで示す。
【図７】図６のマトリクスをブロックごとに示す。
【図８】図７のマトリクスのブロックに番号を付したものを示す。
【図９】図８のマトリクスを３×３のブロックで区分し５×５の三重対角マトリクスで示す。
【図１０】（ａ）、（ｂ）は前記マトリクスを用いた方程式の一例を示す。
【図１１】有限差分法を説明する線図である。
【図１２】有限差分法の一般式を示す。
【図１３】図４（ａ）に示したｋ＝１の第１のケースについての有限差分法に基づいた方程式を示す。
【図１４】図１３の方程式を書き直したものである。
【図１５】ポアソン方程式を、サイクリックリダクション法を用いて計算する処理手順の一例を示すフローチャートである。
【図１６】図１４の方程式から誤差ｒを求める形に書き直したものである。
【図１７】誤差を小さくするために修正すべき量φ’corrを求める方程式である。
【図１８】φの計算結果を示す。
【図１９】２回目の反復計算での誤差ｒの計算結果を示す。
【図２０】２回目の反復計算での誤差を小さくするために修正すべき量φ’corrの計算結果を示す。
【図２１】２回目の反復計算でのφの計算結果を示す。
【図２２】３回目の反復計算での誤差ｒの計算結果を示す。
【図２３】実施例、比較例の計算に用いられた不均一構造格子の三次元空間モデルの斜視図である。
【図２４】格子密度５００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２５】格子密度１０００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２６】格子密度１５００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２７】実施例及び比較例の計算時間を示すグラフである。
【図２８】格子密度１０００万個のケースで、反復回数を増やしたときの誤差の変化を示すグラフである。
【図２９】実施例及び比較例の計算時にアロケートするメモリー容量の比較を示すグラフである。
【図３０】（ａ）、（ｂ）は反復計算の回数を変えたときのポアソン方程式の解を比較した結果を示す。
【図３１】１億個の格子で分割されかつ構造格子の不等分割度も高いケースでのテスト結果を示す。
【図３２】図３１のモデルの収束性評価結果を示す。
【発明を実施するための形態】
【００１９】
以下、本発明の実施の一形態が図面に基づき説明される。
図１は、本実施形態のシミュレーション方法の概略を示すフローチャートである。本実施形態のシミュレーション方法は、言うまでもなくコンピュータ装置を用いて行われ、二次元又は三次元空間（図１では三次元空間）内に、均一又は不均一な格子で構成された不均一構造格子のモデルを設定する工程（Ｓ１）と、この不均一構造格子のモデルに関連付けられた物理量及び条件に基づいて特定式で表されるポアソン方程式を設定する工程（Ｓ２）と、該ポアソン方程式を解くことにより前記物理量を計算する計算工程とを含む（Ｓ３）。
【００２０】
前記二次元又は三次元空間は、例えば直交座標系、円筒座標系又は極座標系が使用される。前記座標系のいずれのかの二次元又は三次元空間内で、構造格子（オイラーメッシュ）のモデルが設定され、それらの情報がコンピュータ装置（図示省略）に記憶される。このモデルは、格子で分割された空間と、そこに与えられた境界条件や変数などを含み、例えば流体シミュレーションでは流体の流れ場を表す。そして、このモデルから、例えば各格子点での圧力や流速等の物理量が計算される。
【００２１】
本明細書において、「均一構造格子」とは、図２（ａ）に示されるように、前記三次元空間Ａに配置されたメッシュの最小単位である各構造格子ａ１…（この格子が区画する最小空間を「セル」と言うことがある。）が、各次元の軸（三次元の場合、ｘ、ｙ及びｚ軸）について、一定の長さで分割されたものとする。この直交座標系の例では、結果的に、各セルは同一形状かつ同一容積を持つ。一方、「不均一構造格子」は、図２（ｂ）に示されるように、三次元空間Ａに配置された構造格子が、各次元の軸（三次元の場合、ｘ、ｙ及びｚ軸）のいずれかについて、非一定の長さで分割されたものとする。この直交座標系の例では、各セルは、結果的に容積又は形状が異なる少なくとも２種類の構造格子ａ２、ａ３を含んで分割される。
【００２２】
不均一構造格子のモデルでは、例えば、図２（ｂ）のように、特定の空間Ａ１をより詳細に解析したい場合、その空間だけを相対的に小さな構造格子ａ２で分割する一方、詳細な解析が特に必要ない空間Ａ２については、相対的に大きな構造格子ａ３で分割することができる。従って、不均一構造格子のモデルは、必要な計算精度を維持しつつ計算コストを減じるのに役立つ。
本発明方法は、均一又は不均一構造格子のいずれについても適用可能である。即ち、k≠constantのモデルについても適用することが可能である。
【００２３】
次に、工程Ｓ２では、前記コンピュータ装置により、下記式（１）のポアソン方程式が設定される。コンピュータ装置は、作業者によって入力された構造格子のモデルの仕様、境界条件及び必要なパラメータに基づいて前記ポアソン方程式を設定する。
【数３】
【００２４】
式（１）において、係数ｋは各格子に定められた定数であり、CFL条件を満たさない、陰解法を用いた場合、非一定である。また、φは計算される変数であり、流体シミュレーションの場合、圧力Ｐや速度三成分（ｕ、ｖ、ｗ）等になるが、これらを代表して一つのφという変数で示される。さらに、"src"はソースタームの略号であり、流体シミュレーションの場合、質量流束（mass flux）の項となり、通常、上記変数φの関数ｆ（φ）となる。
【００２５】
また、上記ｋは、各構造格子の容積VolをＡｐで除して次のように表すこともできる。
ｋ＝Vol／Ａｐ
ここで、Ａｐは、運動量の式（momentum equation）の対角項成分である。また、Ａｐは質量流束の関数であり、この値は一次元の流体シミュレーションの場合、下式のように表すことができ（例えば、「Computational Methods for Fluid Dynamics 3rd Edition」 J.H.Ferziger、M.Peric著）、セル毎に異なる値となる。
Ａｐ＝−（ＡＥ＋ＡＷ）＋ρ／Δｔ
ここで、
ＡＥ＝（ρｕ／２Δｘ）−（Γ／Δｘ2）
Ａｗ＝−（ρｕ／２Δｘ）−（Γ／Δｘ2）
ρ：流体の密度（音速以下では一定）、
ｕ：流体の速度
Γ：流体の粘性
上記式において、Ａpが流体の速度(ｕ)の関数になっている。また、速度ｕは、一様流でない限り、空間で分布を持っている。従って、k≠constant となることがわかる。他方、陽解法の場合、 ｋ＝ρ／Ａｐ、かつ、Ａｐ＝ρ／Δｔ となるため、ｋ＝Δｔとなる。そして、陽解法の場合、通常、時間刻みΔｔは常に一定であるため、k=constant となり得る。
【００２６】
また、構造格子が三次元空間に定義された流体シミュレーションを陰解法で解く場合、Ａｐは次のように求められる。先ず、各セルは、図３（ａ）に示されるように、一つのセル隣り合う６つのセル、即ち、ノース（Ｎ）、サウス（Ｓ）、ウエスト（Ｗ）、イースト（Ｅ）、フロント（Ｆ）及びバック（Ｂ）を持つ。そして、今、各々の隣接するセルの中心（黒丸部分）間の距離を図３（ｂ）に示されるように、「Dist_」に上記括弧内の符号を付して表す。また、セルの上記６つの面Ｎ、Ｓ、Ｗ、Ｅ、Ｆ及びＢは、それぞれの表面積及びこの表面積中心での流体の速度を有し、これらは、「Area_」及び「Velocity_」に上記括弧内の符号を付して表す。そして、各セルの表面での質量流量（「Flux_」に上記括弧内の符号を付して表される。）は、次の式で与えられ、これらは流速が異なるため、質量流量は同一値にはならない。
Flux_Ｎ＝ρ×Area_Ｎ×Velocity_Ｎ
Flux_Ｓ＝ρ×Area_Ｓ×Velocity_Ｓ
Flux_Ｗ＝ρ×Area_Ｗ×Velocity_Ｗ
Flux_Ｅ＝ρ×Area_Ｅ×Velocity_Ｅ
Flux_Ｆ＝ρ×Area_Ｆ×Velocity_Ｆ
Flux_Ｂ＝ρ×Area_Ｂ×Velocity_Ｂ
【００２７】
また、流体の粘性をΓで表すと、非対角項成分（非対角要素：off diagonal element）Ａ＿＊については、以下の式が成り立つ。
A_Ｎ＝−Γ×（Area_Ｎ／Dist_Ｎ）＋Min（Flux_Ｎ，0）
A_Ｓ＝−Γ×（Area_Ｓ／Dist_Ｓ）＋Min（Flux_Ｓ，0）
A_Ｗ＝−Γ×（Area_Ｗ／Dist_Ｗ）＋Min（Flux_Ｗ，0）
A_Ｅ＝−Γ×（Area_Ｅ／Dist_Ｅ）＋Min（Flux_Ｅ，0）
A_Ｆ＝−Γ×（Area_Ｆ／Dist_Ｆ）＋Min（Flux_Ｆ，0）
A_Ｂ＝−Γ×（Area_Ｂ／Dist_Ｂ）＋Min（Flux_Ｂ，0）
なお、上記式において、「Min（Flux_＊,０）」は、引数Flux_＊がマイナス値ならFlux_＊を、プラスの値であれば０をそれぞれ戻り値として返す関数である。
【００２８】
上記の変数を用いれば、構造格子が三次元の場合には、Ａｐは、次の式で表すことができる。
Ａｐ＝−１×（A_Ｎ＋A_Ｓ＋A_Ｗ＋A_Ｅ＋A_Ｆ＋A_Ｂ）
ここで、各面の質量流束は同一ではないため、各々のセルのＡｐも異なったものになり、その結果、k≠constant になる。
【００２９】
陽解法で解く場合、場合、ｋ＝ρ／Ａｐ、 Ａｐ＝ρ／Δｔ となり、k＝Δｔ＝constant となるため、下記式（３）のように、ポアソン方程式の外に出せる。従って、この場合のポアソン方程式は、通常、下式のようにｋの無い形式になる。一方、陰解法を用いて解く場合、ｋ＝Vol／Ａｐ Ａｐ＝−（ＡＥ＋ＡＷ）＋ρ／Δｔとなり、Ａpは流速(ｕ)の関数となる。通常、場所によって速度は変化しているので、ｋは一定にならない。このため、ｋをポアソン方程式の外に出せなくなる。
【数４】
【００３０】
そこで、本発明では、式（１）のようなポアソン方程式を、新たな手法で短時間かつ近似的に計算する。具体的には、ブロックサイクリックリダクション法を反復法のように用いて前記ポアソン方程式を近似的に解くことを特徴とする（計算工程Ｓ３）。
【００３１】
上で述べたように、ブロックサイクリックリダクション法は、直接法であり、１回の計算で殆ど誤差のない解を得ることができる。しかしながら、このブロックサイクリックリダクション法は、本来、k＝constantのポアソン方程式を解く事は可能であるが、k≠constantのポアソン方程式を直接解くことはできない。そこで、本実施形態では、誤差の概念を導入し、この誤差を、ブロックサイクリックリダクション法を用いて反復計算し、かつ、補正によって少しずつ小さくする。つまり、結果的に、k≠constantのポアソン方程式を、前記サイクリックリダクション法を利用して少ない回数の反復計算で解くことを可能にしている。
【００３２】
以下、陰解法では、何故、ブロックサイクリックリダクション法がそのままでは適用できないのかの説明を含めて、具体的な例を用いて説明する。
【００３３】
図４（ａ）には、二次元空間の例を示す。この例では、二次元空間が横３×縦５の計１５の構造格子で分割され、各セルには図のように番号が割り当てられている。なお、この二次元空間の領域は、縦１ｍ×横１ｍの正方形とする。
【００３４】
先ず、２つのケースについてポアソン方程式が設定される。第一のケースは均一構造格子のモデル（図４（ｂ））、第二のケースは不均一構造格子のモデル（図４（ｃ））とする。
【００３５】
ここでの目的の一つは、k≠constantのモデル場合、この例題がブロックサイクリックリダクション法を使って解くことができないことを示すことである。なお、この例題において、構造格子の右端の境界Ｗは、固定された値（＝１０．０）の条件が定義される。また、残りの３つの境界には０の勾配が与えられる。さらに、ソースタームとして、セル１について−１００が、その他のセル２乃至１５については０がそれぞれ定義される。なお、この例題は、単なるマトリクス演算であり、実際の物理現象を再現しているものではない。
【００３６】
また、上記２つのケースについてｋの値がそれぞれ定義される。即ち、図４（ｂ）に示されるように、第１のケースの均一構造格子については全てｋ＝１で一定とする。また、図４（ｃ）に示されるように、第２の不均一構造格子のｋは、中央のセルは全て０．７５、その両側のセルは全て１．０として非一定とする。
【００３７】
また、図５に示される各々の矢印は、これから解くことになるポアソン方程式のマトリクス要素の関係を示している。さらに、図６には、図５の構造格子の関係を示すパターンをマトリクスで示している。縦軸及び横軸の数字は、各セルの番号を意味している。図６から明らかなように、このマトリクスの中において、ゼロでない要素（＃で表されている）は、前記の各々の矢印に対応している。さらに、図７から明らかなように、図６のマトリクスの中には、いくつかの３×３のサブマトリクスがあることが分かる。これらのサブマトリクスの一つを、「ブロック」と呼ぶ。また、これらのブロックには、図８に示されるように、名前を付けることとする。
【００３８】
さらに、図８のマトリクスを３×３のブロックで区分すると、図９に示されるような５×５の三重対角マトリクスと等価であることがわかる。つまり、図１０（ａ）に示される方程式は、図１０（ｂ）の方程式と等価と言える。
【００３９】
次に、図１０（ｂ）の方程式を解くために、本実施形態では、有限差分法（Finite Difference Method）が用いられる。有限差分法では、例えば図１１に示されるように、各セルについて、隣り合う全てのセルとの関係が考慮される。即ち、例えば、セル８について見ると、セル５、７、９及び１１との関係が考慮される。また、有限差分法では、「Computational methods for fluid dynamics」Ferziger J H peric 著の５０頁等で示されるように、図１２で示される一般式が用いられる。
【００４０】
ここで、本例題の場合、Δｘ＝１／３≒０．３３３３３３ｍ、Δｙ＝１／５＝０．２ｍとなる。よって、これらの値と先に述べた境界条件とを用いることにより、図４（ｂ）に示したｋ＝１の第１のケースについては、図１３の方程式に書き表すことができる。同様に、図４（ｃ）に示した第２のケースについては、図１４の方程式に書き表すことができる。
【００４１】
図１３及び図１４から明らかなように、第１のケースでは、図８のブロックＡ２乃至Ａ５に相当する対角のブロックは、Ａi＝ａiＩの形式で表すことができる（なお、ａiはスカラー，Ｉは単位マトリクスである。）。しかしながら、図１４の第２のケースでは、対角項のブロックＡ２乃至Ａ５は、Ａi＝ａiＩの形式で表すことができない。
【００４２】
一方、方程式に（ブロック）サイクリックリダクション法を適用するためには、その方程式が下記の３条件、即ち、Ａ乃至Ｃのブロックが次の形式であることが必要である（この点については、例えば「A Direct Method for the Discrete Solution of Separable Elliptic Equations」 Paul N Swarztrauber著を参照のこと）。
（条件１）Ａi＝ａi×Ｉ
（条件２）Ｂi＝Ｂ＋ｂi×Ｉ
（条件３）Ｃi＝ｃi×Ｉ
ここで、ａi、ｂi及びｃiはスカラー，Ｉは単位マトリクスである。
【００４３】
第１及び第２のケースが、上記３条件を満たすか否かについて検討する。第１のケースでは、条件１乃至３を全て満たしている。従って、均一構造格子（ｋ＝１で一定）のポアソン方程式については、サイクリックリダクション法を適用し、直接法として計算することができる。しかしながら、第２のケース（ｋ≠一定）では、ブロックＡ及びＣは、スカラー×単位マトリクスの形式にはならないため、上記条件１及び３を満たさない。従って、不均一構造格子（ｋ≠一定）については、そのままでは、（ブロック）サイクリックリダクション法を直接法として適用することはできない。
【００４４】
本発明では、上述のような不均一構造格子（ｋ≠一定）のポアソン方程式を、サイクリックリダクション法を反復法的に用いて高速に計算するものであり、その処理手順の一例が図１５に示されている。
【００４５】
図１５に示されるように、本実施形態では、先ず、予め定義された誤差ｒが計算される。該誤差ｒは、本来ゼロとなるべきものであるから、下式（４）のように定義されかつ計算される（ステップＳ３１）。なお、この式（４）は、請求項１の式（２）に対応している。
誤差ｒ＝src−Ａ・φ …式（４）
ここで、srcは前記ソースターム、符号Ａは図１４で示したｋ≠一定のときのポアソン方程式のマトリクス、φは計算される変数であるが初回計算時にはゼロとし、二回目以降は下式（５）で与えられるものとする（この式（５）は、請求項１のφi+1＝φi＋φ'i+1に相当している。）。
φ＝φold＋φ’corr …（５）
また、φoldはφの前回計算値、φ’corrは、前記誤差を小さくするために修正すべき量を示す修正パラメータ（請求項１ではφ'i+1で表されている。）であり、下記式（６）で定義される（なお、請求項１の形式では、式（７）で定義される）。
φ’corr＝Ａuni-1・［ｒ／ｋ_bar］
【数５】
【００４６】
上記式において、Ａuni-1は、図１３に示した第１のケースの均一構造格子のマトリクスＡ（すなわち、第２のケースとは境界条件、格子の分割数等は同一であるが、ｋ＝一定である点のみが異なる第１のケースの均一構造格子のマトリクス）の逆行列であり、ｋ_barは、各セルのｋの平均値である。このφ’corrの計算式において、前記マトリックスＡは、サイクリックリダクション法の適用条件１乃至３の全てを満たすので、その逆行列であるマトリックスＡuni-1も、適用条件１乃至３の全てを満たす。従って、前記修正パラメータφ’corrの計算には、サイクリックリダクション法を直接法として適用でき、１回の計算で精度良く修正パラメータの解を得ることができる。
【００４７】
図１４、図１６の方程式のように、１回目の計算において、誤差ｒはソースタームと等しくなる。この例題では、誤差ｒの総和Σ｜ｒi｜は１００＋９０×５＝５５０である。
【００４８】
次に、本実施形態では、誤差ｒの総和Σ｜ｒi｜が予め定めた許容範囲か否かが判断され（ステップＳ３２）。そして、許容範囲内であれば処理を終える（ステップＳ３２でＹ）。
【００４９】
他方、誤差ｒの総和Σ｜ｒi｜が許容範囲内でない場合（ステップＳ３２でＮ）、上記式に基づいてφ’corr、及び新たなφが計算され、ステップＳ３１に戻って誤差ｒが再度計算される。図１７には、このような方程式が示される。また、図１８には、変数φの計算結果が示されている。
【００５０】
以上のステップが、誤差ｒの総和が許容範囲内になるまで反復計算される。例えば、図１９乃至図２１には、２回目のイタレーションを実行した結果が示される。このときの誤差ｒの総和は６８．７５０９である。
【００５１】
さらに、図２２には、３回目の反復計算を実行したときの誤差ｒの計算式が示されている。このときの誤差ｒの総和は８．８６６９である。このように、わずか３回のイタレーションを実行することで、誤差ｒの総和を顕著に減らすことができる。そして、この誤差ｒの総和が予め定めた許容範囲内になったときの変数φ（請求項１の式（２）ではφi+1）が求めようとする解である。この許容範囲については、シミュレーションの内容や目的に応じて適宜決定することができる。
【００５２】
以上の説明から明らかなように、ポアソン方程式を、陰解法を用いて解くときには、k≠constantとなり、直接法であるサイクリックリダクション法でそのまま解くことはできない。しかしながら、前記した誤差の概念を導入し、この誤差の計算に、不均一構造格子モデル（k≠constant）の境界条件と同じ条件の均一構造格子モデル（即ち、ｋ＝constant）のマトリックスの逆行列を用いることで、サイクリックリダクション法を適用可能とし、かつ、この誤差を補正によって少しずつ小さくすることで、結果的に、k≠constantのポアソン方程式を、サイクリックリダクション法を利用して計算可能としている。従って、ＡＭＧソルバーなどを用いた場合に比して、少ないメモリーで計算時間及び計算コストを大幅に低減することができる。
【００５３】
上記実施形態では、前記構造格子は、二次元のものを示したが、三次元空間として設定されてもよい。この場合、前記計算工程では、予め前記式（１）の方程式をフーリエ変換する。これにより、式（１）のマトリックスの次元数を二次元に落とすことができ、以降は、上で説明した手順がそのまま適用できる。
【００５４】
具体的に説明すると、上記式（７）は、行列式で表すと、式（８）のようになる。なお、式（８）の各要素Ｔiは、例えば図１４のマトリクスＡを意味している。
【数６】
【００５５】
上記式（８）のマトリクスを、フーリエ変換すると、式（９）のマトリクスに変換できる。
【数７】
【００５６】
式（９）の一つ一つの要素を抜き出すと、式（１０）のように表せ、この要素は、ブロックサイクリックリダクション法を用いて上記手順にて計算することができる。なお、
上記式（８）のΦｉは、φが集まって、一回り大きくなったマトリクスを、ブロックとして捉えて表現されたものである。
【００５７】
また、
式（８）で使われている符号ｓｉは、ソースタームであり、［ｒｉ／ｋ_bar］に相当する。さらに、式（８）の４×４の行列が、前記Ａuni-1 に相当する。なお、式（９）及び（１０）のチルダは、フーリエ空間に行ったことを示す。
【実施例】
【００５８】
以下、本発明をより具現化した実施例について説明する。
＜実施例＞
本発明方法を適用して不均一構造格子モデルのポアソン方程式を計算する。
＜比較例１＞
Fluent(ver.6.3.26)に搭載されているＡＭＧ法ソルバー（Additive corrective Multi-Grid, Coarsen by 2 cells, Pre-smoothing iteration=0, post smoothing iterations =2, smoother = ガウス・ザイデル, V-cycle ）を用いて、実施例と同一の不均一構造格子モデルのポアソン方程式を計算する。
【００５９】
＜評価方法＞
実施例及び比較例のテストに使用した計算機の仕様は次の通りである。
ハードウエア：HP社のデスクトップマシン（HP xw9300 Workstation /CPU AMD Dual-core Opteron275 (2.2GHz))
ＯＳ：Linux(x86_64-redhat-linux)
メモリー搭載量：３２GB
実施例のプログラムソフトはC++で作成されている。また、上記ハードウエアのＣＰＵは２つのコアが搭載されているが、計算は一つのコアだけで実行させている。
【００６０】
＜問題設定＞
図２３には、実施例、比較例の計算に用いられた不均一構造格子の三次元空間モデルが示される。該モデルは円筒状をなし、ポアソン方程式の係数ｋは、空間における格子の座標の関数として、ｋ＝√（ｘ2＋ｙ2＋ｚ2）として定義されている。従って、実施例、比較例では、下記の式（１１）のポアソン方程式を計算することになる。
【数８】
【００６１】
また、空間の格子密度は、５００万個、１０００万個及び１５００万個の３つのケースとした。そして、誤差が最初の計算値の１／１０００になるまでに要する反復回数及び計算時間が評価された。
【００６２】
＜テスト結果＞
収束性評価結果（最初のステップでの誤差の値を１で表している。）は図２４乃至図２６に示される通りである。各セルの誤差は先に述べた通り、例えばsrc−Ａ・φで表され、 誤差の総和はΣ｜ｒi｜ で表される。
【００６３】
テストの結果、FluentのＡＭＧ法ソルバーを使用した比較例では、１．０Ｅ-3まで誤差を減らすのに、最低４回の反復計算が必要になった。一方、実施例では、全てのモデルにおいて３回で１．０Ｅ-3よりも小さな誤差になっていることが確認できる。また、計算時間の比較を図２７に示すが、実施例は、比較例の約半分以下の短い時間で収束解が得られていることが分かる。
【００６４】
次に、格子密度１０００万個のモデルについて、さらに反復計算の回数を増やすと、誤差の値がどうなるかを調べた。結果は、図２８に示される通りである。実施例では、リニアに誤差を小さくできるが、比較例では途中からほとんど誤差が減らなくなることが分かる。これも、本発明の利点の一つである。
【００６５】
さらに、図２９には、各実施例での計算にアロケートされるメモリー容量の比較が示される。格子数をＮとすると、比較例の必要なArrayは、７×Ｎ×２＝１４Ｎとなる。一方、 実施例の必要なArrayは、９×Ｎである。従って、実施例は比較例の約６５％のメモリー使用量で計算できることが分かる。
【００６６】
次に、それぞれのポアソン方程式の解を比較した結果が、図３０（ａ）及び（ｂ）に示される。理論的には全ての領域における値が１００になっているのが正解である。図３０のから明らかなように、実施例では３回反復計算させると、既に収束解の１００が得られている。しかし、比較例では２０００回反復させても１００未満であり、誤差が小さくなっていないことがわかる。
【００６７】
さらに、大規模（１億個）の格子で分割されかつ構造格子の不等分割度も高いケースでのテスト結果を示す。不均一構造格子のモデルは図３１に示される通りである。また、ポアソン方程式は、先の例と同じ形である。収束性評価結果（最初のステップでの誤差値を１としている）は図３２に示される通りである。
【００６８】
実施例では、この様に大規模なマトリクスにおいても、たったの５回の反復計算で、誤差を１×Ｅ-8まで減少させることができる。この大規模問題で、この様な速さで収束解を得るのは、本発明手法を用いない限り不可能であろう。さらに、反復計算を継続することにより、反復回数１２回で１．０Ｅ-10を超すところまで誤差を減じうることが分かる。大規模問題で、ここまで誤差を小さくできるポアソン方程式の処理手順は、皆無であると言える。
【符号の説明】
【００６９】
Ａ 構造格子モデル
【技術分野】
【０００１】
本発明は、例えばコンピュータを用いた流体解析等に好適に使用できる構造格子を用いたシミュレーション方法に関する。
【背景技術】
【０００２】
任意の空間を流れる流体を解析する場合、コンピュータを用いたシミュレーションが行われる。このシミュレーションでは、任意の三次元空間を、構造格子を用いてメッシュ化し、各格子点での流体の圧力、温度及び／又は流速等が計算される。そして、この中の圧力に関する方程式は、流体を非圧縮であると定義できる場合、通常、ポアソン方程式を解くことにより行われる。
【０００３】
前記ポアソン方程式の計算方法として、例えばヤコビ法、ガウス・ザイデル法、ＳＯＲ法などの反復法が知られている。反復法は、初期値を与え（手順１）、反復計算を行い（手順２）、収束判定条件を満たす解になるまで手順２を繰り返す（手順３）、ことで行われる。しかしながら、上記方法は、解くべきマトリクスサイズが大きくなると、極端に収束性が悪化し、多くの時間を必要とする。
【０００４】
このため、近年では、反復数がメッシュ規模に依存しないよう、ＡＭＧ（Algebraic Multigrid Methods）法が多用される傾向がある。周知のように、ＡＭＧ法は、直交格子だけでなく非構造格子も解くことができるように、幾何マルチグリッド法を拡張したものである。ＡＭＧ法は、マルチグリッド法と同様、反復数がメッシュサイズに依存しないという特徴があるため、現在、最も高速な連立一次方程式の解法と言え、市販されている汎用ソルバーにも、多く採用されている（例えば、ＡＮＳＹＳ社の「Fluent」、ＣＤ−adapco社の「ＳＴＡＲ−ＣＤ」等）。
【０００５】
しかしながら、ＡＭＧ法は、大規模マトリクスを解くときの収束性に優れる一方、コンピュータに大規模なメモリーサイズを要求する。このため、大きなメモリーを搭載したコンピュータでなければ、ＡＭＧ法で大規模な計算ができないという問題があった。
【０００６】
一方、必要なメモリーサイズが少なくても大規模マトリクスの計算が可能で収束性の良い方法として、直接法の一つであるブロックサイクリックリダクション法が知られている。しかしながら、ブロックサイクリックリダクション法は、ポアソン方程式（例えば請求項１の式１）で、kが不変量(constant value)である場合を前提としており、kが一定値ではなく、マトリクスの中で一カ所でも違う値をとってしまい、マトリクスの外に出せない様なケースでは、ブロックサイクリックリダクション法）を用いて解く事はできない。
【０００７】
従って、ブロックサイクリックリダクション法を正しく適用するためには、シミュレーションで用いられるkの値を一定にすること、つまり、陽解法を用いる（クーラン条件（CFL条件）を守る）必要があり、計算時の時間刻み幅を小さくとる必要があるので、計算コストが大きくなってしまうという欠点がある。
【０００８】
なお、CFL条件（Courant-Friedrichs-Lewy Condition）とは、コンピュータシミュレーションの計算（数値解析）において、「情報が伝播する速さ」を「実際の現象で波や物理量が伝播する速さ」よりも早くしなければならないという必要条件のことである。例えば、離散格子系において波動を扱う場合に、その運動方程式の数値解を求める際に用いる時間ステップ（Δt）の値は、実際の波動が隣り合う格子に伝達するまでの時間よりも小さくなければならない。もし（Δt）の値がCFL条件の上限を超えると、数値発散が生じてしまい、物理的に意味の無い解を得てしまう。また、格子の間隔が小さくなると、時間ステップの上限値も減少する。そして、実際の波(特性速度)の速さをC、計算上の時間ステップ幅をΔt、格子幅をΔxとすれば、 移流方程式等の場合、 Δx/Δt>C (情報が伝播する速さ>実際の波等の速さ)がCFL条件となる。このCFL条件は、陽解法の時間進展を行う際に用いられる条件であり、この条件を回避するために、しばしば陰解法が用いられる場合がある。 陰解法を用いることでCFL条件を回避しかつ緩和できる理由としては様々な説明が存在するが、最も簡潔に説明すると、陽解法は1ステップ前の自分の周りのごくわずかな格子点のみから情報を得て 次の時間の値を決めるのに対して、陰解法は1ステップ前の(ほぼ)全ての格子点の情報を処理して 次の時間の値を決めるためであり、CFL条件におけるΔxが実質的に巨大になるためである。
【先行技術文献】
【特許文献】
【０００９】
【特許文献１】特開平５−１３５０９１号公報
【特許文献２】特開２０００−２９３５０８号公報
【発明の概要】
【発明が解決しようとする課題】
【００１０】
本発明は、以上のような問題点に鑑み案出なされたもので、二次元又は三次元空間内に、均一又は不均一構造格子を設定し、この構造格子に関連付けられた物理量及び条件に基づいて特定形式のポアソン方程式を設定し、このポアソン方程式を、ブロックサイクリックリダクション法を近似的に計算することを基本として、不均一構造格子を前提とするポアソン方程式についても小さなメモリーサイズのコンピュータで迅速に計算することが可能な不均一構造格子を用いたシミュレーション方法を提供することを目的としている。
【課題を解決するための手段】
【００１１】
本発明のうち請求項１記載の発明は、二次元又は三次元空間内に、構造格子のモデルを設定する工程と、
この構造格子のモデルに関連付けられた物理量及び条件に基づいて式（１）のポアソン方程式を設定する工程と、
該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、
前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、
この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、
前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする。
【数２】
【００１２】
また請求項２記載の発明は、構造格子は、二次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法である。
【００１３】
また請求項３記載の発明は、前記構造格子は、三次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法である。
【００１４】
また請求項４記載の発明では、前記計算工程は、予め前記式（２）の方程式をフーリエ変換することによって、次元数を二次元に落とすことを特徴とする請求項３記載の構造格子を用いたシミュレーション方法である。
【発明の効果】
【００１５】
本発明によれば、二次元又は三次元空間内に、構造格子で空間を離散化したモデルを設定する工程と、この構造格子のモデルに関連付けられた物理量及び条件に基づいて特定のポアソン方程式を設定する工程と、該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする。
【００１６】
即ち、ブロックサイクリックリダクション法は直接法であり、本来では、１回の計算で誤差のほとんど無い解が得られるが、k≠constant型のポアソン方程式については、ブロックサイクリックリダクション法の適用条件を満たさず、そのような計算はできない。しかし、本発明では、誤差の概念を導入し、この誤差を、ブロックサイクリックリダクション法を用いて反復計算し、補正によって少しずつ小さくすることを特徴の一つとする。つまり、結果的に、k≠constant型のポアソン方程式を、前記ブロックサイクリックリダクション法を利用して少ない反復回数で解くことを可能にしている。
【００１７】
従って、k≠constantを前提とするポアソン方程式についても、ブロックサイクリックリダクション法を適用できるため、小さなメモリーサイズのコンピュータで計算することができ、かつ、陰解法を適用できる様になるため、大きな時間刻み幅で非定常計算を実施できる様になる。このため、本発明方法では、計算コストの削減にもつながる。これにより、例えば複雑な乱流が生じる空気流体解析等を、高い空間解像度、時間解像度で解くことができ、ゴルフボールのディンプル開発等において、効率化やコスト削減に役立つ。また、ブロックサイクリックリダクション法では、反復回数が少ないため、計算結果に桁落ち等による計算誤差が積み上がる可能性も低く、計算精度の向上にも繋がる。
【図面の簡単な説明】
【００１８】
【図１】本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。
【図２】（ａ）は均一構造格子、（ｂ）は不均一構造格子をそれぞれ説明する斜視図である。
【図３】（ａ）は三次元の構造格子の一つのセルの斜視図、（ｂ）はその正面図である。
【図４】（ａ）は本発明を説明するための二次元の構造格子の例であり、（ｂ）は均一構造格子、（ｃ）は不均一構造格子の例を示す。
【図５】構造格子の互いのマトリクス要素の関係を示す。
【図６】図５の構造格子の関係を示すパターンをマトリクスで示す。
【図７】図６のマトリクスをブロックごとに示す。
【図８】図７のマトリクスのブロックに番号を付したものを示す。
【図９】図８のマトリクスを３×３のブロックで区分し５×５の三重対角マトリクスで示す。
【図１０】（ａ）、（ｂ）は前記マトリクスを用いた方程式の一例を示す。
【図１１】有限差分法を説明する線図である。
【図１２】有限差分法の一般式を示す。
【図１３】図４（ａ）に示したｋ＝１の第１のケースについての有限差分法に基づいた方程式を示す。
【図１４】図１３の方程式を書き直したものである。
【図１５】ポアソン方程式を、サイクリックリダクション法を用いて計算する処理手順の一例を示すフローチャートである。
【図１６】図１４の方程式から誤差ｒを求める形に書き直したものである。
【図１７】誤差を小さくするために修正すべき量φ’corrを求める方程式である。
【図１８】φの計算結果を示す。
【図１９】２回目の反復計算での誤差ｒの計算結果を示す。
【図２０】２回目の反復計算での誤差を小さくするために修正すべき量φ’corrの計算結果を示す。
【図２１】２回目の反復計算でのφの計算結果を示す。
【図２２】３回目の反復計算での誤差ｒの計算結果を示す。
【図２３】実施例、比較例の計算に用いられた不均一構造格子の三次元空間モデルの斜視図である。
【図２４】格子密度５００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２５】格子密度１０００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２６】格子密度１５００万個の実施例及び比較例の収束性評価結果を示すグラフである。
【図２７】実施例及び比較例の計算時間を示すグラフである。
【図２８】格子密度１０００万個のケースで、反復回数を増やしたときの誤差の変化を示すグラフである。
【図２９】実施例及び比較例の計算時にアロケートするメモリー容量の比較を示すグラフである。
【図３０】（ａ）、（ｂ）は反復計算の回数を変えたときのポアソン方程式の解を比較した結果を示す。
【図３１】１億個の格子で分割されかつ構造格子の不等分割度も高いケースでのテスト結果を示す。
【図３２】図３１のモデルの収束性評価結果を示す。
【発明を実施するための形態】
【００１９】
以下、本発明の実施の一形態が図面に基づき説明される。
図１は、本実施形態のシミュレーション方法の概略を示すフローチャートである。本実施形態のシミュレーション方法は、言うまでもなくコンピュータ装置を用いて行われ、二次元又は三次元空間（図１では三次元空間）内に、均一又は不均一な格子で構成された不均一構造格子のモデルを設定する工程（Ｓ１）と、この不均一構造格子のモデルに関連付けられた物理量及び条件に基づいて特定式で表されるポアソン方程式を設定する工程（Ｓ２）と、該ポアソン方程式を解くことにより前記物理量を計算する計算工程とを含む（Ｓ３）。
【００２０】
前記二次元又は三次元空間は、例えば直交座標系、円筒座標系又は極座標系が使用される。前記座標系のいずれのかの二次元又は三次元空間内で、構造格子（オイラーメッシュ）のモデルが設定され、それらの情報がコンピュータ装置（図示省略）に記憶される。このモデルは、格子で分割された空間と、そこに与えられた境界条件や変数などを含み、例えば流体シミュレーションでは流体の流れ場を表す。そして、このモデルから、例えば各格子点での圧力や流速等の物理量が計算される。
【００２１】
本明細書において、「均一構造格子」とは、図２（ａ）に示されるように、前記三次元空間Ａに配置されたメッシュの最小単位である各構造格子ａ１…（この格子が区画する最小空間を「セル」と言うことがある。）が、各次元の軸（三次元の場合、ｘ、ｙ及びｚ軸）について、一定の長さで分割されたものとする。この直交座標系の例では、結果的に、各セルは同一形状かつ同一容積を持つ。一方、「不均一構造格子」は、図２（ｂ）に示されるように、三次元空間Ａに配置された構造格子が、各次元の軸（三次元の場合、ｘ、ｙ及びｚ軸）のいずれかについて、非一定の長さで分割されたものとする。この直交座標系の例では、各セルは、結果的に容積又は形状が異なる少なくとも２種類の構造格子ａ２、ａ３を含んで分割される。
【００２２】
不均一構造格子のモデルでは、例えば、図２（ｂ）のように、特定の空間Ａ１をより詳細に解析したい場合、その空間だけを相対的に小さな構造格子ａ２で分割する一方、詳細な解析が特に必要ない空間Ａ２については、相対的に大きな構造格子ａ３で分割することができる。従って、不均一構造格子のモデルは、必要な計算精度を維持しつつ計算コストを減じるのに役立つ。
本発明方法は、均一又は不均一構造格子のいずれについても適用可能である。即ち、k≠constantのモデルについても適用することが可能である。
【００２３】
次に、工程Ｓ２では、前記コンピュータ装置により、下記式（１）のポアソン方程式が設定される。コンピュータ装置は、作業者によって入力された構造格子のモデルの仕様、境界条件及び必要なパラメータに基づいて前記ポアソン方程式を設定する。
【数３】
【００２４】
式（１）において、係数ｋは各格子に定められた定数であり、CFL条件を満たさない、陰解法を用いた場合、非一定である。また、φは計算される変数であり、流体シミュレーションの場合、圧力Ｐや速度三成分（ｕ、ｖ、ｗ）等になるが、これらを代表して一つのφという変数で示される。さらに、"src"はソースタームの略号であり、流体シミュレーションの場合、質量流束（mass flux）の項となり、通常、上記変数φの関数ｆ（φ）となる。
【００２５】
また、上記ｋは、各構造格子の容積VolをＡｐで除して次のように表すこともできる。
ｋ＝Vol／Ａｐ
ここで、Ａｐは、運動量の式（momentum equation）の対角項成分である。また、Ａｐは質量流束の関数であり、この値は一次元の流体シミュレーションの場合、下式のように表すことができ（例えば、「Computational Methods for Fluid Dynamics 3rd Edition」 J.H.Ferziger、M.Peric著）、セル毎に異なる値となる。
Ａｐ＝−（ＡＥ＋ＡＷ）＋ρ／Δｔ
ここで、
ＡＥ＝（ρｕ／２Δｘ）−（Γ／Δｘ2）
Ａｗ＝−（ρｕ／２Δｘ）−（Γ／Δｘ2）
ρ：流体の密度（音速以下では一定）、
ｕ：流体の速度
Γ：流体の粘性
上記式において、Ａpが流体の速度(ｕ)の関数になっている。また、速度ｕは、一様流でない限り、空間で分布を持っている。従って、k≠constant となることがわかる。他方、陽解法の場合、 ｋ＝ρ／Ａｐ、かつ、Ａｐ＝ρ／Δｔ となるため、ｋ＝Δｔとなる。そして、陽解法の場合、通常、時間刻みΔｔは常に一定であるため、k=constant となり得る。
【００２６】
また、構造格子が三次元空間に定義された流体シミュレーションを陰解法で解く場合、Ａｐは次のように求められる。先ず、各セルは、図３（ａ）に示されるように、一つのセル隣り合う６つのセル、即ち、ノース（Ｎ）、サウス（Ｓ）、ウエスト（Ｗ）、イースト（Ｅ）、フロント（Ｆ）及びバック（Ｂ）を持つ。そして、今、各々の隣接するセルの中心（黒丸部分）間の距離を図３（ｂ）に示されるように、「Dist_」に上記括弧内の符号を付して表す。また、セルの上記６つの面Ｎ、Ｓ、Ｗ、Ｅ、Ｆ及びＢは、それぞれの表面積及びこの表面積中心での流体の速度を有し、これらは、「Area_」及び「Velocity_」に上記括弧内の符号を付して表す。そして、各セルの表面での質量流量（「Flux_」に上記括弧内の符号を付して表される。）は、次の式で与えられ、これらは流速が異なるため、質量流量は同一値にはならない。
Flux_Ｎ＝ρ×Area_Ｎ×Velocity_Ｎ
Flux_Ｓ＝ρ×Area_Ｓ×Velocity_Ｓ
Flux_Ｗ＝ρ×Area_Ｗ×Velocity_Ｗ
Flux_Ｅ＝ρ×Area_Ｅ×Velocity_Ｅ
Flux_Ｆ＝ρ×Area_Ｆ×Velocity_Ｆ
Flux_Ｂ＝ρ×Area_Ｂ×Velocity_Ｂ
【００２７】
また、流体の粘性をΓで表すと、非対角項成分（非対角要素：off diagonal element）Ａ＿＊については、以下の式が成り立つ。
A_Ｎ＝−Γ×（Area_Ｎ／Dist_Ｎ）＋Min（Flux_Ｎ，0）
A_Ｓ＝−Γ×（Area_Ｓ／Dist_Ｓ）＋Min（Flux_Ｓ，0）
A_Ｗ＝−Γ×（Area_Ｗ／Dist_Ｗ）＋Min（Flux_Ｗ，0）
A_Ｅ＝−Γ×（Area_Ｅ／Dist_Ｅ）＋Min（Flux_Ｅ，0）
A_Ｆ＝−Γ×（Area_Ｆ／Dist_Ｆ）＋Min（Flux_Ｆ，0）
A_Ｂ＝−Γ×（Area_Ｂ／Dist_Ｂ）＋Min（Flux_Ｂ，0）
なお、上記式において、「Min（Flux_＊,０）」は、引数Flux_＊がマイナス値ならFlux_＊を、プラスの値であれば０をそれぞれ戻り値として返す関数である。
【００２８】
上記の変数を用いれば、構造格子が三次元の場合には、Ａｐは、次の式で表すことができる。
Ａｐ＝−１×（A_Ｎ＋A_Ｓ＋A_Ｗ＋A_Ｅ＋A_Ｆ＋A_Ｂ）
ここで、各面の質量流束は同一ではないため、各々のセルのＡｐも異なったものになり、その結果、k≠constant になる。
【００２９】
陽解法で解く場合、場合、ｋ＝ρ／Ａｐ、 Ａｐ＝ρ／Δｔ となり、k＝Δｔ＝constant となるため、下記式（３）のように、ポアソン方程式の外に出せる。従って、この場合のポアソン方程式は、通常、下式のようにｋの無い形式になる。一方、陰解法を用いて解く場合、ｋ＝Vol／Ａｐ Ａｐ＝−（ＡＥ＋ＡＷ）＋ρ／Δｔとなり、Ａpは流速(ｕ)の関数となる。通常、場所によって速度は変化しているので、ｋは一定にならない。このため、ｋをポアソン方程式の外に出せなくなる。
【数４】
【００３０】
そこで、本発明では、式（１）のようなポアソン方程式を、新たな手法で短時間かつ近似的に計算する。具体的には、ブロックサイクリックリダクション法を反復法のように用いて前記ポアソン方程式を近似的に解くことを特徴とする（計算工程Ｓ３）。
【００３１】
上で述べたように、ブロックサイクリックリダクション法は、直接法であり、１回の計算で殆ど誤差のない解を得ることができる。しかしながら、このブロックサイクリックリダクション法は、本来、k＝constantのポアソン方程式を解く事は可能であるが、k≠constantのポアソン方程式を直接解くことはできない。そこで、本実施形態では、誤差の概念を導入し、この誤差を、ブロックサイクリックリダクション法を用いて反復計算し、かつ、補正によって少しずつ小さくする。つまり、結果的に、k≠constantのポアソン方程式を、前記サイクリックリダクション法を利用して少ない回数の反復計算で解くことを可能にしている。
【００３２】
以下、陰解法では、何故、ブロックサイクリックリダクション法がそのままでは適用できないのかの説明を含めて、具体的な例を用いて説明する。
【００３３】
図４（ａ）には、二次元空間の例を示す。この例では、二次元空間が横３×縦５の計１５の構造格子で分割され、各セルには図のように番号が割り当てられている。なお、この二次元空間の領域は、縦１ｍ×横１ｍの正方形とする。
【００３４】
先ず、２つのケースについてポアソン方程式が設定される。第一のケースは均一構造格子のモデル（図４（ｂ））、第二のケースは不均一構造格子のモデル（図４（ｃ））とする。
【００３５】
ここでの目的の一つは、k≠constantのモデル場合、この例題がブロックサイクリックリダクション法を使って解くことができないことを示すことである。なお、この例題において、構造格子の右端の境界Ｗは、固定された値（＝１０．０）の条件が定義される。また、残りの３つの境界には０の勾配が与えられる。さらに、ソースタームとして、セル１について−１００が、その他のセル２乃至１５については０がそれぞれ定義される。なお、この例題は、単なるマトリクス演算であり、実際の物理現象を再現しているものではない。
【００３６】
また、上記２つのケースについてｋの値がそれぞれ定義される。即ち、図４（ｂ）に示されるように、第１のケースの均一構造格子については全てｋ＝１で一定とする。また、図４（ｃ）に示されるように、第２の不均一構造格子のｋは、中央のセルは全て０．７５、その両側のセルは全て１．０として非一定とする。
【００３７】
また、図５に示される各々の矢印は、これから解くことになるポアソン方程式のマトリクス要素の関係を示している。さらに、図６には、図５の構造格子の関係を示すパターンをマトリクスで示している。縦軸及び横軸の数字は、各セルの番号を意味している。図６から明らかなように、このマトリクスの中において、ゼロでない要素（＃で表されている）は、前記の各々の矢印に対応している。さらに、図７から明らかなように、図６のマトリクスの中には、いくつかの３×３のサブマトリクスがあることが分かる。これらのサブマトリクスの一つを、「ブロック」と呼ぶ。また、これらのブロックには、図８に示されるように、名前を付けることとする。
【００３８】
さらに、図８のマトリクスを３×３のブロックで区分すると、図９に示されるような５×５の三重対角マトリクスと等価であることがわかる。つまり、図１０（ａ）に示される方程式は、図１０（ｂ）の方程式と等価と言える。
【００３９】
次に、図１０（ｂ）の方程式を解くために、本実施形態では、有限差分法（Finite Difference Method）が用いられる。有限差分法では、例えば図１１に示されるように、各セルについて、隣り合う全てのセルとの関係が考慮される。即ち、例えば、セル８について見ると、セル５、７、９及び１１との関係が考慮される。また、有限差分法では、「Computational methods for fluid dynamics」Ferziger J H peric 著の５０頁等で示されるように、図１２で示される一般式が用いられる。
【００４０】
ここで、本例題の場合、Δｘ＝１／３≒０．３３３３３３ｍ、Δｙ＝１／５＝０．２ｍとなる。よって、これらの値と先に述べた境界条件とを用いることにより、図４（ｂ）に示したｋ＝１の第１のケースについては、図１３の方程式に書き表すことができる。同様に、図４（ｃ）に示した第２のケースについては、図１４の方程式に書き表すことができる。
【００４１】
図１３及び図１４から明らかなように、第１のケースでは、図８のブロックＡ２乃至Ａ５に相当する対角のブロックは、Ａi＝ａiＩの形式で表すことができる（なお、ａiはスカラー，Ｉは単位マトリクスである。）。しかしながら、図１４の第２のケースでは、対角項のブロックＡ２乃至Ａ５は、Ａi＝ａiＩの形式で表すことができない。
【００４２】
一方、方程式に（ブロック）サイクリックリダクション法を適用するためには、その方程式が下記の３条件、即ち、Ａ乃至Ｃのブロックが次の形式であることが必要である（この点については、例えば「A Direct Method for the Discrete Solution of Separable Elliptic Equations」 Paul N Swarztrauber著を参照のこと）。
（条件１）Ａi＝ａi×Ｉ
（条件２）Ｂi＝Ｂ＋ｂi×Ｉ
（条件３）Ｃi＝ｃi×Ｉ
ここで、ａi、ｂi及びｃiはスカラー，Ｉは単位マトリクスである。
【００４３】
第１及び第２のケースが、上記３条件を満たすか否かについて検討する。第１のケースでは、条件１乃至３を全て満たしている。従って、均一構造格子（ｋ＝１で一定）のポアソン方程式については、サイクリックリダクション法を適用し、直接法として計算することができる。しかしながら、第２のケース（ｋ≠一定）では、ブロックＡ及びＣは、スカラー×単位マトリクスの形式にはならないため、上記条件１及び３を満たさない。従って、不均一構造格子（ｋ≠一定）については、そのままでは、（ブロック）サイクリックリダクション法を直接法として適用することはできない。
【００４４】
本発明では、上述のような不均一構造格子（ｋ≠一定）のポアソン方程式を、サイクリックリダクション法を反復法的に用いて高速に計算するものであり、その処理手順の一例が図１５に示されている。
【００４５】
図１５に示されるように、本実施形態では、先ず、予め定義された誤差ｒが計算される。該誤差ｒは、本来ゼロとなるべきものであるから、下式（４）のように定義されかつ計算される（ステップＳ３１）。なお、この式（４）は、請求項１の式（２）に対応している。
誤差ｒ＝src−Ａ・φ …式（４）
ここで、srcは前記ソースターム、符号Ａは図１４で示したｋ≠一定のときのポアソン方程式のマトリクス、φは計算される変数であるが初回計算時にはゼロとし、二回目以降は下式（５）で与えられるものとする（この式（５）は、請求項１のφi+1＝φi＋φ'i+1に相当している。）。
φ＝φold＋φ’corr …（５）
また、φoldはφの前回計算値、φ’corrは、前記誤差を小さくするために修正すべき量を示す修正パラメータ（請求項１ではφ'i+1で表されている。）であり、下記式（６）で定義される（なお、請求項１の形式では、式（７）で定義される）。
φ’corr＝Ａuni-1・［ｒ／ｋ_bar］
【数５】
【００４６】
上記式において、Ａuni-1は、図１３に示した第１のケースの均一構造格子のマトリクスＡ（すなわち、第２のケースとは境界条件、格子の分割数等は同一であるが、ｋ＝一定である点のみが異なる第１のケースの均一構造格子のマトリクス）の逆行列であり、ｋ_barは、各セルのｋの平均値である。このφ’corrの計算式において、前記マトリックスＡは、サイクリックリダクション法の適用条件１乃至３の全てを満たすので、その逆行列であるマトリックスＡuni-1も、適用条件１乃至３の全てを満たす。従って、前記修正パラメータφ’corrの計算には、サイクリックリダクション法を直接法として適用でき、１回の計算で精度良く修正パラメータの解を得ることができる。
【００４７】
図１４、図１６の方程式のように、１回目の計算において、誤差ｒはソースタームと等しくなる。この例題では、誤差ｒの総和Σ｜ｒi｜は１００＋９０×５＝５５０である。
【００４８】
次に、本実施形態では、誤差ｒの総和Σ｜ｒi｜が予め定めた許容範囲か否かが判断され（ステップＳ３２）。そして、許容範囲内であれば処理を終える（ステップＳ３２でＹ）。
【００４９】
他方、誤差ｒの総和Σ｜ｒi｜が許容範囲内でない場合（ステップＳ３２でＮ）、上記式に基づいてφ’corr、及び新たなφが計算され、ステップＳ３１に戻って誤差ｒが再度計算される。図１７には、このような方程式が示される。また、図１８には、変数φの計算結果が示されている。
【００５０】
以上のステップが、誤差ｒの総和が許容範囲内になるまで反復計算される。例えば、図１９乃至図２１には、２回目のイタレーションを実行した結果が示される。このときの誤差ｒの総和は６８．７５０９である。
【００５１】
さらに、図２２には、３回目の反復計算を実行したときの誤差ｒの計算式が示されている。このときの誤差ｒの総和は８．８６６９である。このように、わずか３回のイタレーションを実行することで、誤差ｒの総和を顕著に減らすことができる。そして、この誤差ｒの総和が予め定めた許容範囲内になったときの変数φ（請求項１の式（２）ではφi+1）が求めようとする解である。この許容範囲については、シミュレーションの内容や目的に応じて適宜決定することができる。
【００５２】
以上の説明から明らかなように、ポアソン方程式を、陰解法を用いて解くときには、k≠constantとなり、直接法であるサイクリックリダクション法でそのまま解くことはできない。しかしながら、前記した誤差の概念を導入し、この誤差の計算に、不均一構造格子モデル（k≠constant）の境界条件と同じ条件の均一構造格子モデル（即ち、ｋ＝constant）のマトリックスの逆行列を用いることで、サイクリックリダクション法を適用可能とし、かつ、この誤差を補正によって少しずつ小さくすることで、結果的に、k≠constantのポアソン方程式を、サイクリックリダクション法を利用して計算可能としている。従って、ＡＭＧソルバーなどを用いた場合に比して、少ないメモリーで計算時間及び計算コストを大幅に低減することができる。
【００５３】
上記実施形態では、前記構造格子は、二次元のものを示したが、三次元空間として設定されてもよい。この場合、前記計算工程では、予め前記式（１）の方程式をフーリエ変換する。これにより、式（１）のマトリックスの次元数を二次元に落とすことができ、以降は、上で説明した手順がそのまま適用できる。
【００５４】
具体的に説明すると、上記式（７）は、行列式で表すと、式（８）のようになる。なお、式（８）の各要素Ｔiは、例えば図１４のマトリクスＡを意味している。
【数６】
【００５５】
上記式（８）のマトリクスを、フーリエ変換すると、式（９）のマトリクスに変換できる。
【数７】
【００５６】
式（９）の一つ一つの要素を抜き出すと、式（１０）のように表せ、この要素は、ブロックサイクリックリダクション法を用いて上記手順にて計算することができる。なお、
上記式（８）のΦｉは、φが集まって、一回り大きくなったマトリクスを、ブロックとして捉えて表現されたものである。
【００５７】
また、
式（８）で使われている符号ｓｉは、ソースタームであり、［ｒｉ／ｋ_bar］に相当する。さらに、式（８）の４×４の行列が、前記Ａuni-1 に相当する。なお、式（９）及び（１０）のチルダは、フーリエ空間に行ったことを示す。
【実施例】
【００５８】
以下、本発明をより具現化した実施例について説明する。
＜実施例＞
本発明方法を適用して不均一構造格子モデルのポアソン方程式を計算する。
＜比較例１＞
Fluent(ver.6.3.26)に搭載されているＡＭＧ法ソルバー（Additive corrective Multi-Grid, Coarsen by 2 cells, Pre-smoothing iteration=0, post smoothing iterations =2, smoother = ガウス・ザイデル, V-cycle ）を用いて、実施例と同一の不均一構造格子モデルのポアソン方程式を計算する。
【００５９】
＜評価方法＞
実施例及び比較例のテストに使用した計算機の仕様は次の通りである。
ハードウエア：HP社のデスクトップマシン（HP xw9300 Workstation /CPU AMD Dual-core Opteron275 (2.2GHz))
ＯＳ：Linux(x86_64-redhat-linux)
メモリー搭載量：３２GB
実施例のプログラムソフトはC++で作成されている。また、上記ハードウエアのＣＰＵは２つのコアが搭載されているが、計算は一つのコアだけで実行させている。
【００６０】
＜問題設定＞
図２３には、実施例、比較例の計算に用いられた不均一構造格子の三次元空間モデルが示される。該モデルは円筒状をなし、ポアソン方程式の係数ｋは、空間における格子の座標の関数として、ｋ＝√（ｘ2＋ｙ2＋ｚ2）として定義されている。従って、実施例、比較例では、下記の式（１１）のポアソン方程式を計算することになる。
【数８】
【００６１】
また、空間の格子密度は、５００万個、１０００万個及び１５００万個の３つのケースとした。そして、誤差が最初の計算値の１／１０００になるまでに要する反復回数及び計算時間が評価された。
【００６２】
＜テスト結果＞
収束性評価結果（最初のステップでの誤差の値を１で表している。）は図２４乃至図２６に示される通りである。各セルの誤差は先に述べた通り、例えばsrc−Ａ・φで表され、 誤差の総和はΣ｜ｒi｜ で表される。
【００６３】
テストの結果、FluentのＡＭＧ法ソルバーを使用した比較例では、１．０Ｅ-3まで誤差を減らすのに、最低４回の反復計算が必要になった。一方、実施例では、全てのモデルにおいて３回で１．０Ｅ-3よりも小さな誤差になっていることが確認できる。また、計算時間の比較を図２７に示すが、実施例は、比較例の約半分以下の短い時間で収束解が得られていることが分かる。
【００６４】
次に、格子密度１０００万個のモデルについて、さらに反復計算の回数を増やすと、誤差の値がどうなるかを調べた。結果は、図２８に示される通りである。実施例では、リニアに誤差を小さくできるが、比較例では途中からほとんど誤差が減らなくなることが分かる。これも、本発明の利点の一つである。
【００６５】
さらに、図２９には、各実施例での計算にアロケートされるメモリー容量の比較が示される。格子数をＮとすると、比較例の必要なArrayは、７×Ｎ×２＝１４Ｎとなる。一方、 実施例の必要なArrayは、９×Ｎである。従って、実施例は比較例の約６５％のメモリー使用量で計算できることが分かる。
【００６６】
次に、それぞれのポアソン方程式の解を比較した結果が、図３０（ａ）及び（ｂ）に示される。理論的には全ての領域における値が１００になっているのが正解である。図３０のから明らかなように、実施例では３回反復計算させると、既に収束解の１００が得られている。しかし、比較例では２０００回反復させても１００未満であり、誤差が小さくなっていないことがわかる。
【００６７】
さらに、大規模（１億個）の格子で分割されかつ構造格子の不等分割度も高いケースでのテスト結果を示す。不均一構造格子のモデルは図３１に示される通りである。また、ポアソン方程式は、先の例と同じ形である。収束性評価結果（最初のステップでの誤差値を１としている）は図３２に示される通りである。
【００６８】
実施例では、この様に大規模なマトリクスにおいても、たったの５回の反復計算で、誤差を１×Ｅ-8まで減少させることができる。この大規模問題で、この様な速さで収束解を得るのは、本発明手法を用いない限り不可能であろう。さらに、反復計算を継続することにより、反復回数１２回で１．０Ｅ-10を超すところまで誤差を減じうることが分かる。大規模問題で、ここまで誤差を小さくできるポアソン方程式の処理手順は、皆無であると言える。
【符号の説明】
【００６９】
Ａ 構造格子モデル
【特許請求の範囲】
【請求項１】
二次元又は三次元空間内に、格子で構成された構造格子のモデルを設定する工程と、
この構造格子のモデルに関連付けられた物理量及び条件に基づいて式（１）のポアソン方程式を設定する工程と、
該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、
前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、
この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、
前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする構造格子を用いたシミュレーション方法。
【数１】
【請求項２】
前記構造格子は、二次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法。
【請求項３】
前記構造格子は、三次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法。
【請求項４】
前記計算工程は、予め前記式（２）の方程式をフーリエ変換することによって、次元数を二次元に落とすことを特徴とする請求項３記載の構造格子を用いたシミュレーション方法。
【請求項１】
二次元又は三次元空間内に、格子で構成された構造格子のモデルを設定する工程と、
この構造格子のモデルに関連付けられた物理量及び条件に基づいて式（１）のポアソン方程式を設定する工程と、
該ポアソン方程式に基づいて前記物理量を計算する計算工程とを含み、
前記計算工程は、式（２）で定義される誤差を、ブロックサイクリックリダクション法を用いて計算する誤差計算工程と、
この誤差が予め定めた許容範囲内か否かを判断する工程とを含み、
前記誤差が予め定められた許容範囲内になるまで、変数φを補正パラメータを用いて補正しながら前記誤差計算工程を繰り返すことにより、前記式（１）のポアソン方程式を近似的に解くことを特徴とする構造格子を用いたシミュレーション方法。
【数１】
【請求項２】
前記構造格子は、二次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法。
【請求項３】
前記構造格子は、三次元空間として設定される請求項１記載の構造格子を用いたシミュレーション方法。
【請求項４】
前記計算工程は、予め前記式（２）の方程式をフーリエ変換することによって、次元数を二次元に落とすことを特徴とする請求項３記載の構造格子を用いたシミュレーション方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３２】
【図３０】
【図３１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３２】
【図３０】
【図３１】
【公開番号】特開２０１２−８３９５８（Ｐ２０１２−８３９５８Ａ）
【公開日】平成２４年４月２６日（２０１２．４．２６）
【国際特許分類】
【出願番号】特願２０１０−２２９８２４（Ｐ２０１０−２２９８２４）
【出願日】平成２２年１０月１２日（２０１０．１０．１２）
【出願人】（０００１８３２３３）住友ゴム工業株式会社 (3,458)
【Ｆターム（参考）】
【公開日】平成２４年４月２６日（２０１２．４．２６）
【国際特許分類】
【出願日】平成２２年１０月１２日（２０１０．１０．１２）
【出願人】（０００１８３２３３）住友ゴム工業株式会社 (3,458)
【Ｆターム（参考）】
[ Back to top ]