(57)【要約】【課題】新しい演算操作と手順を有し且つ簡略で判りやすい物理量のシミュレーション方法を提供する。【解決手段】計算空間に格子点を配置し変数を格子点上に置く。３格子点上の変数によって2次関数近似を行うことにより格子点上の１階微分係数を演算する。近似解を求めるには非線形計画法の罰金関数を使用する。全格子点について、物理量を記述する微分方程式を罰金関数の制約条件式に組み込み、制約条件式に先で演算した１階微分係数を代入する。罰金関数において、制約条件式の残差を取る目的関数をゼロ近傍に収束させることにより近似解を求める。【選択図】図２
2006039850.tif 000002
【特許請求の範囲】
【請求項１】
所定の微分方程式の変数として記述される物理量を演算するシミュレーション方法において、
計算空間に配置され且つ変数が与えられた格子点のうちの少なくとも３つの格子点上の変数により近似関数を作成し、各格子点における少なくとも一次の微分係数を求める第一のステップと、
前記微分方程式を組み込んだ制約条件式と前記制約条件式の残差を求める目的関数から構成される罰金関数の前記制約条件式に前記微分係数を代入し、前記目的関数がゼロ近傍に収束するように前記変数を変化させながら、前記変数の近似解を求める第二のステップとを備えることを特徴とするシミュレーション方法。
【請求項２】
前記第二のステップにおいて、各格子点の変数を変化させる方向を決めるための方向係数を演算し、その際、複数の格子点それぞれの変数の方向係数を並列演算で求めることを特徴とする請求項１に記載のシミュレーション方法。
【請求項３】
各格子点に複数の変数が与えられ、各変数に対応する微分方程式がそれぞれ組み込まれた複数の制約条件式を含む罰金関数が作成される場合、各制約条件式に重みが付されることを特徴とする請求項１又は２に記載のシミュレーション方法。
【請求項４】
所定の微分方程式の変数として記述される物理量を演算するシミュレーションプログラムにおいて、
計算空間に配置され且つ変数が与えられた格子点のうちの少なくとも３つの格子点上の変数により近似関数を作成し、各格子点における少なくとも一次の微分係数を求める第一の処理と、
前記微分方程式を組み込んだ制約条件式と前記制約条件式の残差を求める目的関数から構成される罰金関数の前記制約条件式に前記微分係数を代入し、前記目的関数がゼロ近傍に収束するように前記変数を変化させながら、前記変数の近似解を求める第二の処理とを備えることを特徴とするシミュレーションプログラム。
【請求項５】
前記第二の処理において、各格子点の変数を変化させる方向を決めるための方向係数を演算し、その際、複数の格子点それぞれの変数の方向係数を並列演算で求めることを特徴とする請求項４に記載のシミュレーションプログラム。
【請求項６】
各格子点に複数の変数が与えられ、各変数に対応する微分方程式がそれぞれ組み込まれた複数の制約条件式を含む罰金関数が作成される場合、各制約条件式に重みが付されることを特徴とする請求項４又は５に記載のシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、所定の物理量をシミュレーション演算するシミュレーション方法及びそれをコンピュータ上で実行するためのコンピュータプログラムに関し、特に、従来の計算が離散化の手順を経るのに替わって、罰金関数を導入した新しい演算操作と手順を提供するものである。
【背景技術】
【０００２】
例えば、流体の物理量である流速や圧力などをシミュレーション演算する場合、従来の演算方法は、流速や圧力など各物理量の状態を表す数学的モデル（微分方程式）を選択し、それを適切な離散化の方法により離散化することで微分方程式を代数方程式の系で近似して解くものであった。流体の物理量の場合、数学的モデルは、偏微分方程式や積分形微分方程式と境界条件を検討し、圧縮性か非圧縮性か、粘性か非粘性か、層流か乱流かなどの適切なモデルを解析対象の事象より選択される。
【０００３】
また、離散化の方法としては、差分法、有限体積法、有限要素法が主要なものとしてあげられる。代数方程式の解法としては、ガウス消去法、LU分解やＴＤＭＡ法などの直接解法や反復解法などがある。また、代数方程式の解法においては、直接解法を取り入れた反復解法を用いて逐次求めることが一般的である。
【０００４】
図１は、２次元直交座標系のナビエ・ストークスの式にｋ−εの乱流モデルを取り入れた流体の数学的モデルに対する従来のシミュレーション演算の手順を示す。変数（物理量）は、ｘ方向流速を u 、ｙ方向流速を v 、圧力を p 、乱流エネルギーを ｋ 、散逸率をεとしている。差分法を使い変数を支配する代数方程式を、逐次演算、反復し収束させ近似解を求める。
【０００５】
具体的には、ステップＳ１では、変数配列の大きさ及び方程式の定数を設定する。ステップＳ２では、格子点に添字番号と座標値を与える。ステップＳ３では、変数ｕ、ｖ、ｐ、ｋ、εの適切な初期値と境界条件を与える。ステップＳ４では、ｘ方向運動方程式を与えられた出発値より例えばＴＤＭＡ法（Tri-Diagonal Matrix Algorithm）にて数回緩和し、変数ｕの新しい近似値を得る。ステップＳ５では、ｙ方向運動方程式を与えられた出発値よりＴＤＭＡ法にて数回緩和し、変数ｖの新しい近似値を得る。ステップＳ６では、連続の式とｘ、ｙ方向運動方程式より圧力補正式を、与えられた圧力出発値よりＴＤＭＡ法にて数回緩和して得る。そして、これを基に変数ｐ、ｕ、ｖの新しい近似値を演算する。ステップＳ７では、乱流エネルギー輸送方程式を、与えられた出発値よりＴＤＭＡ法にて数回緩和し、乱流エネルギーｋの新しい近似値を得る。ステップＳ８では、散逸率εの輸送方程式を、与えられた出発値よりＴＤＭＡ法にて数回緩和し、散逸率εの新しい近似値を得る。ステップＳ９では、流体物性の粘度に乱流粘性係数を加えて実行粘性を演算する。ステップＳ１０では、各変数の収束性を判断し、収束が不完全な場合は、ステップＳ４に戻る。また、収束性が満足された場合は、ステップＳ１１において、各変数の最終値を出力する。
【発明の開示】
【発明が解決しようとする課題】
【０００６】
このように、従来の演算方法では、物理量を記述する微分方程式（数学的モデル）を離散化し、微分方程式を代数方程式の系で近似して解くことが必要である。離散化については、解くべき課題の幾何形状に関係し複雑となり、演算量が膨大となる。また、解くべき課題に対して、その都度代数方程式を個々に作成する必要があり、汎用性がない。
【０００７】
そこで、本発明の目的は、汎用性が有り且つより簡便に物理量を演算するシミュレーション方法及びそれをコンピュータ上で実行するシミュレーションプログラムを提供することにある。
【課題を解決するための手段】
【０００８】
上記目的を達成するための本発明における第一のシミュレーション方法は、所定の微分方程式の変数として記述される物理量を演算するシミュレーション方法において、計算空間に配置され且つ変数が与えられた格子点のうちの少なくとも３つの格子点上の変数により近似関数を作成し、各格子点における少なくとも一次の微分係数を求める第一のステップと、前記微分方程式を組み込んだ制約条件式と前記制約条件式の残差を求める目的関数から構成される罰金関数の前記制約条件式に前記微分係数を代入し、前記目的関数がゼロ近傍に収束するように前記変数を変化させながら、前記変数の近似解を求める第二のステップとを備えることを特徴とする。
【０００９】
また、本発明における第二のシミュレーション方法は、上記第一のシミュレーション方法の第二のステップにおいて、各格子点の変数を変化させる方向を決めるための方向係数を演算し、その際、複数の格子点それぞれの変数の方向係数を並列演算で求めることを特徴とする。
【００１０】
さらに、本発明における第三のシミュレーション方法は、上記第一又は第二のシミュレーション方法において、各格子点に複数の変数が与えられ、各変数に対応する微分方程式がそれぞれ組み込まれた複数の制約条件式を含む罰金関数が作成される場合、各制約条件式に重みが付されることを特徴とする。
【発明の効果】
【００１１】
本発明では、非線形計画法における罰金関数の制約条件式に、物理量の状態を記述する微分方程式を直接組み込むことにより、離散化の過程を経ないで近似解を求めることができる。
【００１２】
また、微分方程式がそのまま罰金関数の一部になるため、解くべき課題が、幾何形状と境界条件が千変万化したものであっても、幾何形状と境界条件ごとに罰金関数を作成する必要なく、一つ作成し、それに格子点の配置と境界条件を与えるだけよく、処理のパターン化が可能である。離散化を行う場合は、解くべき課題に対して、その都度代数方程式を個々に作成する必要があるが、その必要がなくなる。
【００１３】
変数の方向性を決める方向係数の演算について、並列演算を行えば演算時間を短縮できる。
【発明を実施するための最良の形態】
【００１４】
以下、図面を参照して本発明の実施の形態について説明する。しかしながら、かかる実施の形態例が、本発明の技術的範囲を限定するものではない。
【００１５】
本発明の実施の形態における物理量のシミュレーション方法は、ある物理量の状態を表す微分方程式を、非線形計画法の罰金関数に等式として制約条件式に組み込み、罰金関数を最大化する演算操作により、目的関数である前記等式の残差を０近傍に収束させることにより近似解を得る演算手順を持つものである。数学的モデルである微分方程式を罰金関数に直接等式として組み込むため、従来の計算で必要であった離散化の過程が不要となる。
【００１６】
さらに具体的には、まず、計算空間に格子点が配置される。格子点はコロケート格子を利用する。コロケート格子は、例えば流体の数値計算において、物理量である速度、圧力などの全ての変数を格子点上に配置する方法である。そして、計算空間に配置された格子点上に変数が配置され、３格子点以上を使って１階以上微分可能な変数の近似関数を作成する。通常は、３格子点による２次関数近似を行い、１階の微分係数を作成する。もちろん、３次関数以上の近似を行い、２階以上の微分係数を用いてもよい。各格子点について、この微分係数を罰金関数の制約条件式に使うので３格子点が最小幾何形状単位となる。最小幾何形状単位が小さいため、格子点の疎密の配置が比較的自由に行える。二次関数近似の場合、この近似関数は３格子点の範囲内で有効な１階微分係数を有するため、３格子点において罰金関数の制約条件式の１階微分式の部分に直接代入できる。このため離散化が不要となる。格子点が端点の場合は端点の１次微分を使用するが、中間の場合は中間点の１階微分を使用することにより精度が向上する。
【００１７】
非線形計画法である罰金関数は、下記数式１に示される目的関数と制約条件式にて構成する。ただし、数式１は１変数の場合を表す。右辺第１項が目的関数であり、第２項が制約条件式である。制約条件式Ai(x)が等式の場合はh[Ai(x)]＝１とし、Ai(x)に非負条件を課すならばAi(x)＞０およびAi(x)＝０にてh[Ai(x)]＝０、
Ai(x)＜０にてh[Ai(x)]＝１とする。Rは計算ステップの進行とともに増大させる。
Mは制約条件式の個数である。
【００１８】
【数１】
2006039850.tif 000003
罰金関数は、制約条件のもとで目的関数を最大化するために使われる。最大化する手法については、SUMT(Sequential Unconstrained Maximization Techniqe)法が知られている。SUMT法については、例えば、「オペレーションズ・リサーチ入門  ４＝整数計画法と非線形計画法」（H.M.ワグナー著、森村英典ほか監訳、培風館、昭和６１年、１６５−１７０頁）に詳細に説明されている。例えば、流体の物理量演算に罰金関数を適用することは、計算空間に配置された各格子点について制約条件式に連続の式やナビエ・ストークスの運動方程式などを代入し、目的関数としてはこれら等式の残差の２乗和の負をとる。罰金関数を最大化する演算操作により目的関数を最大化する。即ち、残差を０の近傍に収束させることによって近似解を得る。
【００１９】
計算の収束状況を把握するために、各制約条件の収束率を導入する。２次元直交座標系における連続の式では、下記数式２に示す各格子点の残差の和を、各格子点の成分項∂u/∂x と∂v/∂ｙの各絶対値の和で除した値を収束率Rnと定義する。Nnは計算空間において連続の式を考慮すべき格子点の数である。他の運動方程式なども同様とする。この収束率は計算の収束状況の指標となる。
【００２０】
【数２】
2006039850.tif 000004
収束率の導入により目的関数を修正する。下記数式３に示すように各格子点における制約条件式の残差Onに各格子点の収束率Rnを乗じた数値の２乗On2・Rn2を全格子点についての負の和をとることをもって修正目的関数とする。目的関数の修正についてはこの他にも色々考えられるであろう。Nは格子点の個数である。
【００２１】
【数３】
2006039850.tif 000005
多変数になる数値流体の計算における罰金関数c(x|R)mは、制約条件式に重みWiをつけることにより対応できる。数式１を変形した下記数式４に多変数の罰金関数を示す。数式４では制約条件の個数iと目的関数の個数が同数となっているが、厳密には制約条件には非負条件式があるのでこれらを除く等式部分が目的関数の個数となる。もちろん、制約条件、目的関数とも各格子点についての総和とする。
【００２２】
制約条件式の重みWiの操作について、２つの制約条件式A1(x)、A2(x)がある場合を示す。基準となる制約条件式Ａ1(X)とする。W1の初期値W1sを設定する。以下の操作を罰金関数最大化のステップの最初（後述のステップサイズｔk＝０）に行う。制約条件式の残差の２乗和ΣO12、ΣO22を求める。W2は、ΣO12／ΣO22とし残差の２乗和の大きさを制約条件式A1(x)に合わせる。収束率R1＞R2ならば、W1s＋δを新しい重みとし、R1＜R2ならば、W1s−δとする。ただし、最小重みW1minは有限の値とする。以上により、収束率R1の大小によって重みW1が自動調節される。
【００２３】
【数４】
2006039850.tif 000006
本発明の演算において演算負荷の大きな部分を占めるのが、変数の方向性を決める部分である。新しい変数の点xjk+1は現在の変数の点xjkを使って下記数式５に示される。
【００２４】
【数５】
2006039850.tif 000007
ｔkはステップサイズであり,ｄjkは最急上昇法を使えば、罰金関数のxjについての偏微分係数であり方向係数とも呼ばれる。方向係数を解析的に導出することは困難であり下記の数式６を使って数値計算する。
【００２５】
【数６】
2006039850.tif 000008
罰金関数c(x|R) xj+δは、変数xjをxj+δに変化させるのみでその他の全ての変数は同じである。このため、方向係数ｄjkの演算においては逐次計算は不要であり、並列演算を行えば計算時間を短縮できる。これは罰金関数の最大化を行うのに通常数回の罰金関数の演算を行うのに較べ、方向係数ｄjkの演算には格子点上の全変数の個数に等しい回数の罰金関数の演算を行うため演算負荷が大きいからである。
【００２６】
以下、上述した本発明の実施の形態における演算手法を用いた流体物理量のシミュレーション方法の手順について説明する。
【００２７】
図２は、本発明の実施の形態における物理量のシミュレーション方法の手順の例を示す図であって、図１と同等の条件により、２次元直交座標系のナビエ・ストークスの式にｋ−εの乱流モデルを取り入れた流体の数学的モデルに対するシミュレーション演算の手順を例に説明する。以下に説明するように、罰金関数の演算は、５種類の罰金関数があり変数uについての罰金関数が最大化されれば、変数vについての罰金関数の最大化を行うというように、逐次最大化させることにより収束させる手順をとるものとする。
【００２８】
ステップＳ１０１では、変数配列の大きさ及び方程式の定数を決める。ステップＳ１０２では、格子点に添字番号と座標値を与え、３格子点のつながりを定義する。ステップＳ１０３では、変数ｕ、ｖ、ｐ、ｋ、εの適切な初期値と境界条件を与える。ステップＳ１０４では、３格子点上の変数より二次近似関数を作成し、格子点の一次微分を演算する。ステップＳ１０５では、ステップＳ１０４で演算した一次微分係数を罰金関数の制約条件式に代入し、罰金関数の値を得る。罰金関数には次の５種類がある。（１）連続の式とｘ方向運動方程式を制約条件に組み込み、変数ｕのみを変化させる罰金関数。（２）連続の式とｙ方向運動方程式を制約条件に組み込み、変数ｖのみを変化させる罰金関数。（３）ｘ、ｙ運動方程式を制約条件に組み込み、変数ｐのみ変化させる罰金関数。（４）乱流エネルギーｋの輸送方程式を制約条件に組み込み、変数ｋのみ変化させる罰金関数。（５）散逸率εの輸送方程式を制約条件に組み込み、変数εのみ変化させる罰金関数。
【００２９】
ステップＳ１０６では、罰金関数の最大化に使用する変数の方向係数を演算する。演算は、罰金関数の最大化の過程において最初の１回のみ行い、同一最大化の範囲内では、同じ方向係数を使用する。ステップＳ１０７では、罰金関数を最大化する。その方法は、出発値の変数から方向係数にステップサイズを乗じて新しい変数値とし罰金関数を演算する。ステップサイズを操作することによって、最大の罰金関数値を探す。ステップＳ１０８では、罰金関数は凸の負の関数と考えられるため、あらかじめ３点のステップサイズで罰金関数を求めておき、２次関数近似して頂点のステップサイズを新しい値とする。ステップＳ１０９では、罰金関数が最大近傍に収束したなら、そのときのステップサイズより新しい変数値を演算する。ステップＳ１１０では、罰金関数は、上述の通り５つあり、通常、ステップＳ１０５における（１）〜（５）の順に演算する。定義した収束率に基づいて演算の過程を監視し、収束率が一定の値になれば、収束したものとして演算を終了する。ステップＳ１１１では、収束率の改善ができなくなったことにより演算操作を終了し、各変数の最終値を出力する。
【００３０】
図３は、本発明の実施の形態におけるシミュレーション装置のブロック構成図である。本実施の形態例シミュレーション装置は、例えば、パーソナルコンピュータのような汎用コンピュータ装置であってもよいし、シミュレーション専用のコンピュータ装置であってもよい。
【００３１】
シミュレーション装置１０のシミュレーションプログラム格納部１０２は、本実施の形態例における上記演算手順を実行するためのシミュレーションプログラムを記憶する記憶手段であり、ＣＰＵ１００が、本実施の形態に特徴的なシミュレーションプログラムを実行する。ＲＡＭ１０４は一時的なデータを格納する記憶手段である。シミュレータプログラムの実行による処理結果は、データとしてハードディスク装置に記憶されるとともに、ディスプレイ１０６に画像として出力されたり、又は、プリンタ１０８から紙媒体で出力される。
【００３２】
ところで、流体力学における２次元の流れ函数ψは、ｘ方向流速uとｙ方向流速ｖより数式７の関係がある。この関係式を制約条件式とすれば、下記数式８の罰金関数c(ψ|R)が作成できる。これを最大化すれば流れ函数の近似解を求めることができる。Nfは格子点の個数である。後にこの罰金関数を使用して実際に流れ函数を求める。
【００３３】
【数７】
2006039850.tif 000009
【００３４】
【数８】
2006039850.tif 000010
罰金関数を使用した場合の解の振動安定性をみる。１次元の対流・拡散問題において速度u、拡散係数Γが一定と仮定すると下記数式９が得られる。
【００３５】
【数９】
2006039850.tif 000011
この数式は解析解を有し、ペクレ数PeをρuＬ／Γとし、Pe＝５，１０，２０，４０における罰金関数を使用した近似解と解析解を図４に示す。格子点は等間隔の２１点とした。近似解に振動の発生はないが、Peが更に大きくなれば振動が発生する傾向がある。一般的には本実施の形態における演算方法にコロケート格子が使用できるといえる。
【００３６】
上記本実施の形態では、流体の物理量についてのシミュレーション演算について詳細を記述したが、他の分野（例えば電磁気）の物理量においても本実施の形態における演算操作と手順を適用できる。
【実施例】
【００３７】
実施例として、ｋ−ε乱流モデルと壁法則を取り入れた２次元のナビエ・ストークスの運動方程式を数学的モデルとする、円弧翼の流れを計算する。円弧翼の曲率と長さの比は２．５：１とした。境界条件は入り口の流速と圧力のみを一定としその他は自由流れとした。レイノルズ数は１．２５×１０5である。収束率は、連続の式で７．３％であり、ｘ、ｙ方向の運動方程式で７．１％と７．３％である。乱れエネルギーｋの輸送方程式と散逸率εの輸送方程式では１３．４％と１０．５％である。
【００３８】
図５に円弧翼形状と格子点の配置を示す。図６に速度ベクトルｕ、ｖを示す。円弧翼の上下面に壁法則を適用した。図７に罰金関数を使用して求めた流れ函数による流線を示す。図８には、円弧翼の壁法則で配置した格子点の圧力ｐの分布を示す。
【００３９】
図９は乱流エネルギーｋの分布を表すものであり、乱流エネルギーｋの大きさを円の径の大きさで格子点上に示している。円弧翼廻りに最大部があるのがわかる。図１０には、散逸率εの分布を乱流エネルギーｋと同様の方法で表す。
【図面の簡単な説明】
【００４０】
【図１】従来における物理量演算のシミュレーション方法の手順を示す図である。
【図２】本発明の実施の形態における物理量のシミュレーション方法の手順の例を示す図である。
【図３】本発明の実施の形態におけるシミュレーション装置のブロック構成図である。
【図４】１次元対流・拡散方程式の解を示す図である。
【図５】円弧翼の流れの翼形状と流入条件と格子点の配置を表す図である。
【図６】円弧翼の流れの演算結果であり、速度ベクトルｕ、ｖを示す図である。
【図７】図５の速度ベクトルｕ、ｖから、罰金関数によって計算した流れ函数φの流線図を示す。
【図８】円弧翼の流れの演算結果であり、円弧翼上下面の境界層に配置した格子点の圧力ｐを示す。
【図９】円弧翼の流れの演算結果であり、ｋ−εモデルの乱流エネルギーｋの分布図である。
【図１０】円弧翼の流れの演算結果であり、ｋ−εモデルの散逸率εの分布図である。
【符号の説明】
【００４１】
１：円弧翼、２：格子点、１０：シミュレーション装置、１００：ＣＰＵ、１０２：シミュレーションプログラム格納部、１０４：ＲＡＭ、１０６：ディスプレイ、１０８：プリンタ
【図１】
2006039850.tif 000012
【図２】
2006039850.tif 000013
【図３】
2006039850.tif 000014
【図４】
2006039850.tif 000015
【図５】
2006039850.tif 000016
【図６】
2006039850.tif 000017
【図７】
2006039850.tif 000018
【図８】
2006039850.tif 000019
【図９】
2006039850.tif 000020
【図１０】
2006039850.tif 000021