(57)【特許請求の範囲】
【請求項１】
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
前記運動方程式の数値計算において、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定するシミュレーション方法。
【請求項２】
前記運動方程式の数値計算のタイムステップごとに、温度の速度を０に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を０に設定する請求項１に記載のシミュレーション方法。
【請求項３】
前記運動方程式の数値計算のタイムステップごとに、前記熱伝導現象のシミュレーションで求められた温度の計算値の分布が定常状態に到達したか否かを判定する工程と、
前記熱伝導現象のシミュレーションで求められた温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う工程と
を有する請求項１または２に記載のシミュレーション方法。
【請求項４】
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
前記運動方程式の数値計算を行う際に、
前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、
前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻し、
シミュレーション開始からのタイムステップ数が基準値以上になったら、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に固定するシミュレーション方法。
【請求項５】
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
前記運動方程式の数値計算を行う際に、
前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、
前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻すとともに、前記運動方程式の質量に相当するパラメータである仮想質量を、現時点の仮想質量より増加させるシミュレーション方法。
【請求項６】
前記仮想質量が上限値を超えると、前記仮想質量を前記上限値に設定する請求項５に記載のシミュレーション方法。
【請求項７】
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行うシミュレーション装置であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式を記憶する記憶装置と、
前記記憶装置に記憶されている前記運動方程式の数値演算を行うことにより、前記シミュレーション対象物の温度分布を求める中央処理ユニットと、
前記中央処理ユニットで行われた数値演算の結果を出力する出力装置と
を有し、
前記運動方程式は減衰係数を含み、
前記中央処理ユニットは、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行うシミュレーション装置。
【請求項８】
前記中央処理ユニットは、前記運動方程式の数値計算のタイムステップごとに、温度の速度を０に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を０に設定する請求項７に記載のシミュレーション装置。
【請求項９】
前記中央処理ユニットは、
前記運動方程式の数値計算のタイムステップごとに、前記運動方程式の数値演算で算出された温度の計算値の分布が定常状態に到達したか否かを判定する工程と、
前記運動方程式の数値演算で算出された温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う請求項７または８に記載のシミュレーション装置。
【請求項１０】
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションをコンピュータに実行させるシミュレーションプログラムであって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行う機能、及び
前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行う機能を実現するシミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション対象の熱伝導現象をシミュレーションするシミュレーション方法、シミュレーション装置、及びシミュレーションプログラムに関する。
【背景技術】
【０００２】
古典力学や量子力学等を基に、コンピュータを用いて物質科学全般の現象を探るための方法として、分子動力学法（ＭＤ法）に基づくシミュレーション方法が知られている。マクロスケール系を扱うためにＭＤ法を発展させた繰り込み群分子動力学法（ＲＭＤ法）が提案されている。シミュレーション対象の熱伝導現象について、ＭＤ法やＲＭＤ法では、通常、格子振動（フォノン）による熱伝導しか取り扱うことができない。このため、熱伝導現象に自由電子が大きく寄与する金属の熱伝導現象をＭＤ法やＲＭＤ法でシミュレーションすると、実際の現象からかけ離れた結果が得られる場合が多い。
【０００３】
下記の特許文献１に、ＭＤ法やＲＭＤ法において熱伝導現象をより正確に扱うことができるシミュレーション方法が開示されている。特許文献１に開示されたシミュレーション方法においては、粒子に温度のパラメータを持たせる。ある時刻における各粒子の温度に基づいて熱伝導方程式を解くことにより、次の時刻における各粒子の温度を算出する。熱伝導方程式を解くことにより、実際のシミュレーション対象の比熱、熱伝導率を反映したシミュレーションが可能になる。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特許第５４４１４２２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
一般的に、熱伝導現象による温度変化は、機構的な現象や弾性現象（以下、「機構弾性現象」という。）による粒子位置の変化に比べて緩やかである。熱伝導現象をシミュレーションするときの時間刻み幅は、機構弾性現象をシミュレーションするときの時間刻み幅に比べて長く設定される。ところが、熱伝導現象と機構弾性現象との両方を取扱う連成シミュレーションにおいては、機構弾性現象をシミュレーションするための短い時間刻み幅で、熱伝導現象もシミュレーションされる。
【０００６】
温度場の解析においては、一般的に、定常状態における温度分布に興味が示される。温度分布が定常状態に達するまで、短い時間刻み幅で熱伝導現象をシミュレーションすると、計算時間が長大化する。
【０００７】
本発明の目的は、熱伝導現象を取扱う際の計算時間の長大化を抑制することが可能なシミュレーション方法を提供することである。本発明の他の目的は、計算時間の長大化を抑制することが可能なシミュレーション装置を提供することである。本発明のさらに他の目的は、計算時間の長大化を抑制することが可能なシミュレーションプログラムを提供することである。
【課題を解決するための手段】
【０００８】
本発明の一観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
前記運動方程式の数値計算において、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定するシミュレーション方法が提供される。
【０００９】
本発明の他の観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行うシミュレーション装置であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式を記憶する記憶装置と、
前記記憶装置に記憶されている前記運動方程式の数値演算を行うことにより、前記シミュレーション対象物の温度分布を求める中央処理ユニットと、
前記中央処理ユニットで行われた数値演算の結果を出力する出力装置と
を有し、
前記運動方程式は減衰係数を含み、
前記中央処理ユニットは、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行うシミュレーション装置が提供される。
【００１０】
本発明のさらに他の観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションをコンピュータに実行させるシミュレーションプログラムであって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行う機能、及び
前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行う機能を実現するシミュレーションプログラムが提供される。
【発明の効果】
【００１１】
熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、熱伝導現象のシミュレーションを行うことができる。運動方程式の数値計算において、種々の効率化、高速化手法を適用することができる。
【図面の簡単な説明】
【００１２】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２は、図１に示したフローチャートのステップ１３のより詳細なフローチャートである。
【図３】図３は、シミュレーション対象物のある注目点における温度の速度のシミュレーション結果の一例を示すグラフである。
【図４】図４は、実施例によるシミュレーション方法を実行するシミュレーション装置のブロック図である。
【図５】図５は、熱伝導方程式と等価な運動方程式の減衰係数を、本来の値とは異なる値に設定しても、定常状態における正確な温度を求めることが可能なことを説明するためのモデルの斜視図である。
【図６】図６Ａは、シミュレーション対象物の斜視図であり、図６Ｂは、定常状態に達した時のシミュレーション対象物の長さ方向に関する温度分布のシミュレーション結果を示すグラフである。
【図７】図７は、シミュレーション対象物（図６Ａ）の左端に位置する１つの粒子の温度の速度のシミュレーション結果を示すグラフである。
【図８】図８は、参考例によるシミュレーション方法のフローチャートである。
【図９】図９は、シミュレーション対象である簡易的なベアリングの断面図である。
【図１０】図１０Ａは、図９の参考例によるシミュレーション方法を用いて求めた全体温度の変化を示すグラフであり、図１０Ｂは、減衰係数を本来の値に固定してシミュレーションを行って求めた全体温度の変化を示すグラフである。
【図１１】図１１は、他の実施例によるシミュレーション方法のフローチャートである。
【図１２】図１２は、図１１に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示すグラフである。
【図１３】図１３は、さらに他の実施例によるシミュレーション方法のフローチャートである。
【図１４】図１４は、図１３に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示すグラフである。
【発明を実施するための形態】
【００１３】
図１に、実施例によるシミュレーション方法のフローチャートを示す。実施例においては、複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う。機構弾性現象のシミュレーションでは、シミュレーション対象物を粒子の集まりで表現し、ある時間刻み幅Δｔで各粒子の動解析を行う。この動解析には、例えば、ＭＤ(molecular dynamics)法、ＳＰＨ(Smoothed Particle Hydrodynamics)法、ＭＰＳ(Moving Particle Simulation)法、ＤＥＭ(Distinct Element Method)等の粒子法を適用可能である。各粒子をメッシュの節点とみなせば、有限要素法（ＦＥＭ）を適用することも可能である。
【００１４】
熱伝導現象のシミュレーションでは、各粒子に温度のパラメータを持たせ、動解析と同一の時間刻み幅Δｔで温度の数値計算を行う。温度の数値計算を行う際に、熱伝導方程式を解く代わりに、熱伝導方程式と等価な運動方程式を解くことにより、温度場の解析を行う。ここで、「等価」とは、熱伝導方程式の空間的な温度分布の項と、温度の時間微分の項とに関して、熱伝導方程式と同形の式に変形可能であることを意味する。運動方程式における質点の位置を、「温度」に対応させると、温度の時間変化の傾き（温度の一次微分）が運動方程式における速度に対応する。本明細書において、温度の時間変化の傾きを「温度の速度」ということとする。
【００１５】
温度場の解析において、温度変化の時定数を可変とすることにより、温度変化の傾きを実際の温度変化の傾きより急峻にする。これにより、温度の計算結果が定常状態に到達するまでの時間を短くすることができる。ここで、「定常状態」とは、各粒子の位置における温度の速度が０であることを意味する。
【００１６】
シミュレーションにおいて温度変化の傾きを急峻にすると、シミュレーションで求められた温度がオーバシュートして、振動波形が現れる場合がある。実施例では、この振動を短期に緩和させる手法を適用することにより、温度場が定常状態に達するまでの時間の短縮化を図っている。振動を短期に緩和させる手法として、例えばＦＩＲＥ(fast inertial relaxation engine)と呼ばれる方法を用いることができる。ＦＩＲＥに関しては、Erik Bitzek, Pekka Koskinen, Franz Gahler, Michael Moseler, and Peter Gumbsch, “Structural Relaxation Made Simple”, Physical Review Letters, 97, 170201 (2006)に説明されている。
【００１７】
図１を参照しながら、実施例によるシミュレーション方法について説明する。ステップ１０において、シミュレーション対象物を構成する各粒子の位置、速度、温度（粒子の位置における物体の温度）等、シミュレーションを行うに当たっての初期条件、及び拘束条件（境界条件）を設定する。
【００１８】
ステップ１１において、各粒子の動解析を１タイムステップだけ実行する。例えば、各粒子の運動方程式を数値的に解くことにより、時間刻み幅Δｔ後における粒子の位置及び速度を求める。ステップ１２において、動解析の結果に基づいて、発熱量を算出する。発熱量Ｐは、異なる部材の接触部において、以下の式を用いて算出される。
Ｐ＝μ｜Ｆ・ｖｒ｜
ここで、μは摩擦係数、Ｆは接触部に加わる力、ｖｒは相対速度である。
【００１９】
ステップ１３において、各粒子の発熱量を用いて、熱伝導方程式と等価な運動方程式を解くことにより、時間刻み幅Δｔ後における各粒子の温度を算出する。
【００２０】
図２に、ステップ１３のより詳細なフローチャートを示す。まず、ステップ１３１において、温度場の計算結果が定常状態に達したか否かを判定する。定常状態に達したか否かの判定は、例えば各粒子の温度の計算値の変化量と、判定閾値とを比較することにより行う。全ての粒子の温度の計算値の変化量が、判定閾値以下である場合、温度場が定常状態に達したと判定される。
【００２１】
温度場が定常状態に達していない場合、ステップ１３２において、熱伝導方程式と等価な運動方程式の減衰項を規定するパラメータ（減衰係数）を、本来の値よりも小さくするか、０にする。温度場が定常状態に達するまでは、運動方程式の減衰係数を、本来の値よりも小さくするか、０にして運動方程式を解く。これにより、温度の変化を、本来の変化より急峻にすることができる。
【００２２】
温度場が定常状態に達した場合、ステップ１３３において、熱伝導方程式と等価な運動方程式の減衰係数を、本来の値と等しくする。減衰係数の本来の値は、後に詳細に説明するように、シミュレーション対象物の密度と比熱との積に基づいて決定される。
【００２３】
ステップ１３２またはステップ１３３の後、ステップ１３４において、熱伝導方程式と等価な運動方程式の数値計算を１タイムステップ実行する。このとき、運動方程式の減衰係数として、ステップ１３２またはステップ１３３で設定された値を用いる。ステップ１３２の後、図１のフローチャートに戻る。
【００２４】
図１のステップ１４において、シミュレーションの終了条件が満たされているか否かを判定する。シミュレーションの終了条件として、例えば、動解析の実行ステップ数等を採用することができる。シミュレーションの終了条件が満たされている場合には、シミュレーションを終了する。シミュレーションの終了条件が満たされていない場合には、ステップ１５において、温度場が定常状態に達したか否かを判定する。
【００２５】
ステップ１５では、例えば各粒子の温度の速度が基準値以下になったら、温度場が定常状態に達したと判定される。温度場が定常状態に達したと判定された場合、ステップ１１に戻って、動解析の次の１タイムステップを実行する。温度場が定常状態に達していないと判定された場合、ステップ１６において、温度の速度を０に設定する条件が満たされているか否かを判定する。この判定は、熱伝導方程式と等価な運動方程式の数値計算のタイムステップごとに行われる。
【００２６】
ステップ１６において、判定条件が満たされていないと判定された場合、ステップ１１に戻って、動解析を１タイムステップ実行する。ステップ１６において、判定基準が満たされていると判定された場合、ステップ１７において、各粒子の温度の速度を０に設定した後、ステップ１１に戻って、動解析を１タイムステップ実行する。ステップ１６及びステップ１７の手法は、ＦＩＲＥと呼ばれる。
【００２７】
次に、図３を参照して、ＦＩＲＥについて説明する。図３は、シミュレーション対象物のある注目点における温度の速度のシミュレーション結果の一例を示す。細い実線ａは、通常の熱伝導方程式を数値的に解いて得られる温度の速度の一例を示す。図３には、時間の経過とともに、温度が徐々に低下する例が示されている。
【００２８】
図３の破線ｂは、ステップ１３（図１）において熱伝導方程式と等価な運動方程式を数値的に解いて得られる温度の速度の一例を示す。後述するように、運動方程式の減衰係数は、定常状態における温度場にほとんど影響を与えない。従って、定常状態における温度分布を求めたい場合には、減衰係数を任意に設定することが可能である。
【００２９】
温度場が定常状態になるまでの数値計算のステップ数を少なくするために、減衰係数を小さくすることが好ましい。減衰係数が小さいと、同一の外力が加わっても温度の計算値の変化が大きくなる。このとき、弾性項の影響が相対的に大きくなるため、温度の速度の計算値に振動が現れる。
【００３０】
減衰係数を小さくして、熱伝導方程式と等価な運動方程式を解くと、破線ｂで示すように、数値演算によって求められた温度の計算結果の変化が急峻になる。さらに、温度の計算値は、振動しながら定常状態の温度に近づく。
【００３１】
ＦＩＲＥを適用した場合の温度の計算値を太い実線ｃで示す。時刻ｔ１の時点で、ステップ１６（図１）の条件が満たされたと仮定する。時刻ｔ１において、温度の速度が０に再設定（ステップ１７）される。このため、時刻ｔ１以降における温度の計算値の速度が、一旦緩やかになる。その結果、温度の計算値の振動が抑制され、定常状態に達するまでの時間が短くなる。
【００３２】
ステップ１６（図１）の条件の一例について説明する。各粒子において、温度の１次微分と２次微分と質量との積を求める。すべての粒子について求められた積の和を求める。ステップ１６の条件として、この和が負であるこという条件が採用される。粒子ごとに求められた積が負であることは、各粒子の温度の速度の向きと加速度の向きとが異なること、すなわち温度変化が抑制されることを意味する。従って、この条件が満たされた場合、温度場が定常状態に近づきつつあると判断することができる。
【００３３】
図４に、実施例によるシミュレーション方法を実行するシミュレーション装置（コンピュータ）のブロック図を示す。中央処理ユニット２０、記憶装置２１、入力装置２２、及び出力装置２３がデータバスライン２４を介して相互に接続されている。
【００３４】
入力装置２２から、シミュレーションを実行するための前提条件となるパラメータが入力される。このパラメータには、例えば、粒子間ポテンシャル、外場、密度、比熱、熱伝導係数等が含まれる。ステップ１０（図１）で決定される粒子の初期状態、及び拘束条件も、入力装置２２から入力される。入力装置２２には、例えばキーボードを用いることができる。
【００３５】
記憶装置２１に、制御プログラム、シミュレーションプログラム、及びシミュレーションで数値計算される運動方程式を定義する種々のパラメータ等が格納されている。シミュレーションプログラムは、図１に示したシミュレーションを行う機能を実現する。中央処理ユニット２０は、制御プログラム及びシミュレーションプログラムを実行する。シミュレーションが終了すると、中央処理ユニット２０は、シミュレーション結果を出力装置２３に出力する。例えば、三次元座標系内に各粒子の位置を表示し、各粒子の温度を色の濃淡で表すことができる。
【００３６】
次に、ステップ１３（図１）において、熱伝導方程式と等価な運動方程式を解くことにより温度場の解析を行うことができる理由について説明する。
【００３７】
古典的スカラー場φのラグランジアンＬは、例えばAdvanced Quantum Mechanics, (Addison Wesley, 1967)に説明されているように、下記の式で表される。
【数１】
0006679161.tif 000002
ここでσｓ［φ（ｘ）］は拘束条件を表している。λｓは、ラグランジェの未定定数である。δは、ディラックのδ関数である。δに付された指数３は、及びｄｘに付された指数３は、三次元を意味する。Σは、ｘに関して足し合わせることを意味する。ｘがｘｓ以外のときは、δ関数の値は０である。式（１）において、ρ（ｘ）、κ、Ｑ（ｘ）は、特定の物理量を意味していない。
【００３８】
拘束条件σｓ［φ（ｘ）］として、ディレクリ境界条件、ノイマン境界条件を採用することができる。ディレクリ境界条件は、下記の式で表すことができる。
【数２】
0006679161.tif 000003
式（２）は、ｘ＝ｘｓにおいてφ（ｘ）が定数であることを意味する。
【００３９】
ノイマン境界条件は、下記の式で表すことができる。
【数３】
0006679161.tif 000004
式（３）は、ｘ＝ｘｓにおいてφ（ｘ）の傾きが定数であることを意味する。φ（ｘｓ）をφｓと表している。∇ｓは、ｘ＝ｘｓの位置における傾きを表す。φｓ０は、拘束条件に基づく定数である。
【００４０】
連続体のラグランジアンを定義する式（１）を、一次元の規則格子上で離散化し、離散化したスカラー場をφｉと表記する。離散化したスカラー場φｉのラグランジアンＬは、以下の式で表される。
【数４】
0006679161.tif 000005
ここで、ΔＶｉは線形体積素を表す。∇ｉ±１は差分演算子を表す。差分演算子は、以下の式で定義される。
【数５】
0006679161.tif 000006
ここで、Δは、粒子間の距離を表す。
【００４１】
スカラー場φｉを、粒子の属性、例えば位置、速度、温度等とみなせば、φｉは以下のラグランジェの方程式に従う。
【数６】
0006679161.tif 000007
ここで、γｉは減衰係数である。式（４）、式（６）から、ｉ番目の粒子に対する以下の運動方程式が得られる。
【数７】
0006679161.tif 000008
ここで、式（７）の右辺第１項は弾性項を記述し、第２項は減衰項を記述し、第３項は外力を記述し、第４項は拘束条件を記述している。右辺の第１項のΣのｊ＝±１は、ｉ番目の粒子の両側の粒子について和を取ることを意味する。右辺の第４項のΣのｓは、拘束条件の数に相当する。
【００４２】
（ｄ２／ｄｔ２）φｉが０になる定常状態において、式（７）は、以下のように変形できる。
【数８】
0006679161.tif 000009
【００４３】
式（８）の右辺第１項は、κ∇２φ（ｘ）の差分近似を表しているから、式（８）は拡散方程式に帰着する。式（８）は、式（７）の左辺を０と仮定して得られたものである。従って、式（７）の左辺のρｉを十分小さすれば、式（６）の解が、拡散方程式（８）の解に漸近することがわかる。すなわち、定微分方程式の数値解法を用いて式（６）を解くことにより、近似的に、式（８）の拡散方程式を解くことができる。定微分方程式の数値解法として、ベレル(Verlet)法、リープフロッグ法、ギア(Gear)法等が挙げられる。
【００４４】
熱伝導方程式は、以下の式で表される。
【数９】
0006679161.tif 000010
ここで、ρは密度、Ｃｖは比熱、Ｔは温度、ｔは時間、Ｋは熱伝導率、Ｑは単位体積当たりの発熱量を表す。
【００４５】
式（８）のφｉ、γｉ、κｉ、及びＱｉを、それぞれ温度Ｔ、密度ρと比熱Ｃｖとの積、熱伝導率Ｋ、及び発熱量Ｑとみなせば、式（８）が式（９）の熱伝導方程式と同形であることがわかる。具体的には、式（８）の空間的なφｉの分布の項、及びφｉの時間微分の項を、それぞれ式（９）の空間的な温度分布の項、及び温度の時間微分の項に対応付けることができる。このため、ρｉが十分小さいという条件の下で、式（７）の運動方程式が式（９）の熱伝導方程式と等価であるといえる。
【００４６】
ここまでは、古典的スカラー場φが一次元の規則格子上で離散化された場合について説明した。以上の議論を、三次元の複雑形状に拡張するためには、式（７）に有限体積法を適用すればよい。式（７）に有限体積法を適用すると、以下の式が得られる。
【数１０】
0006679161.tif 000011
ここで、式（１０）の右辺第１項のΣは、ｉ番目の粒子に接触している全ての粒子について足し合わせることを意味する。ΔＳｉｊはボロノイ多面体ｉとｊとの接触面の面積を表す。ｍｉは下記の式で表される。
【数１１】
0006679161.tif 000012
ｍｉは、運動方程式の質量に相当するパラメータであるため、ｍｉを仮想質量と呼ぶこととする。
【００４７】
仮想質量ｍｉは、数値積分の安定条件から、以下の下限条件を満足する必要がある。
【数１２】
0006679161.tif 000013
【００４８】
時間刻み幅ｄｔについては、機構弾性現象の解析を行う際の粒子間距離（メッシュの寸法）及び音速に基づいて、最適値を決定することできる。機構弾性現象に基づいて決定された時間刻み幅ｄｔ、及び式（１２）から、仮想質量ｍｉの下限値が決定される。さらに、上述のように、ρｉが十分小さいという条件の下で、式（６）の解が拡散方程式（８）の解に漸近する。このため、ρｉは、シミュレーション対象の真の密度ρよりも十分小さいことが好ましい。すなわち、以下の式が満たされることが好ましい。
【数１３】
0006679161.tif 000014
【００４９】
式（１０）のγｉを、密度ρｉと比熱Ｃｖｉとの積に置き換え、φｉを温度Ｔｉに置き換えることにより、以下の式が得られる。Ｃｖｉは、ｉ番目の粒子の位置の比熱を表す。
【数１４】
0006679161.tif 000015
【００５０】
運動方程式と同形の式（１４）を解くことにより、温度場のシミュレーションを行うことができる。式（１０）のφｉを式（１４）においてＴｉに置き換えたため、式（１４）のｍｉの次元は、式（１０）、（１１）のｍｉの次元とは異なる。仮想質量ｍｉの次元は、式（１４）の右辺のΔＶｉを含む項から、密度ρｉ、比熱Ｃｖ、体積素ΔＶｉ、及び時間刻み幅ｄｔの積の次元と同一であることがわかる。すなわち、仮想質量ｍｉの次元は、［Ｊ／Ｋ］［ｓ］となる。
【００５１】
式（１２）から、仮想質量ｍｉの満たすべき条件は、以下の式で与えられる。
【数１５】
0006679161.tif 000016
【００５２】
さらに、仮想質量ｍｉの満たすべき条件は、種々の数値計算を行った結果、経験的に以下の式で与えられる。
【数１６】
0006679161.tif 000017
【００５３】
式（１５）及び式（１６）から、一例として、ｍｉを下記の式のように設定することが推奨される。
【数１７】
0006679161.tif 000018
【００５４】
以上、説明したように、熱伝導方程式（９）と等価な運動方程式（１４）を解くことにより温度場の解析を行うことができる。
【００５５】
ステップ１３（図１）において、運動方程式（７）の減衰係数γｉ、式（１４）においてはρｉＣｖｉを０にすれば、温度Ｔの速度の算出結果が大きくなる。すなわち、温度Ｔの時間変化が急峻になる。このため、温度の計算値が定常状態における温度に達するまでの時間を短くすることができる。
【００５６】
次に、ステップ１３２（図２）において、減衰係数γｉまたはρｉＣｖｉを、本来の値とは異なる値、または０に設定しても、定常状態における正確な温度を求めることが可能なことを示す。
【００５７】
図５に示した棒状のモデルを考える。体積ｄＶｉの中央領域３０に、単位体積あたりの発熱量Ｑｉの熱源が存在する。中央領域３０から、面積Ｓｉの界面を通って、両側の中間領域３１、３２に向って熱伝導が生じる。中間領域３１、３２の端に、それぞれ端部領域３３、３４が連続している。端部領域３３、３４が、それぞれ面積Ｓｊの端面において温度Ｔ０の熱浴に接している。
【００５８】
中央領域３０から一方の中間領域３１に流れる単位面積、単位時間あたりの熱エネルギｊＡは、以下の式で表される。
【数１８】
0006679161.tif 000019
端部領域３３の温度をＴで表し、端部領域３３と熱浴との熱伝達係数をｈで表すと、端部領域３３から熱浴に流れる単位面積、単位時間あたりの熱エネルギｊＢは、以下の式で表される。
【数１９】
0006679161.tif 000020
【００５９】
中央領域３０から中間領域３１に流れる単位時間あたりの熱エネルギと、端部領域３３から熱浴に流れる単位時間あたりの熱エネルギとは等しい。このため、以下の式が成立する。
【数２０】
0006679161.tif 000021
【００６０】
式（１８）〜式（２０）から、以下の式が導出される。
【数２１】
0006679161.tif 000022
【００６１】
式（２１）から、定常状態における表面温度は、熱源からの発熱量と、シミュレーション対象物の表面積のみで決まることがわかる。定常状態における表面温度は、密度及び比熱に依存しない。従って、式（１４）の減衰係数ρｉＣｖｉを任意の値、または０として運動方程式（１４）を解いても、定常状態における表面温度を求めることが可能である。
【００６２】
減衰係数ρｉＣｖｉを小さくすると、図３の破線ｂで示したように、温度の計算値に振動が現れる。減衰係数ρｉＣｖｉの大きさに応じて振動が減衰し、温度の計算値が定常状態の温度に収束する。
【００６３】
減衰係数ρｉＣｖｉを０にすると振動が減衰せず、いつまでも振動が継続してしまう。実施例においては、ＦＩＲＥを適用することにより、図３の実線ｃで示したように、温度の計算値の振動を短時間で収束させている。
【００６４】
温度場が定常状態に達するまでの期間、すなわち減衰係数を本来の値と異なる値に設定してステップ１３を実行している期間は、熱伝導方程式と等価な運動方程式を解くときの時間刻み幅の実質的な意味が、ステップ１１（図１）における動解析で適用される時間刻み幅と異なる。このため、シミュレーション対象物の機構弾性現象の時間発展と、熱伝導現象の時間発展とが、非同期で進行することになる。温度場が定常状態に達した後、減衰係数として本来の値を設定すると、シミュレーション対象の機構弾性現象の時間発展と、熱伝導現象の時間発展とが同期することになる。
【００６５】
上記実施例によるシミュレーション方法の効果を確認するために、物体の温度場のシミュレーションを行った。次に、図６Ａ、図６Ｂ、及び図７を参照して、このシミュレーション結果について説明する。
【００６６】
図６Ａに、シミュレーション対象物の斜視図を示す。四角柱状の物体４１の左側の端面と、四角柱状の物体４２の右側の端面とが熱接触している。物体４１及び物体４２の長さは共に２０ｍｍであり、長さ方向に直交する断面は正方形であり、面積は同一である。
【００６７】
物体４１の右端から物体４１に熱エネルギが流入する。物体４１の左端と物体４２の右端との接触面を介して、物体４１から物体４２に熱が伝達される。物体４２の左端から、熱エネルギが流出する。物体４１に流入する単位時間あたりの熱エネルギと、物体４２の左端から流出する単位時間あたりの熱エネルギとは等しい。このため、物体４１と物体４２との熱収支が０に保たれる。物体４１及び４２の温度場は、最終的に定常状態に達する。すなわち、温度の速度が０になる。
【００６８】
図６Ｂに、定常状態に達した時の物体４１及び４２の長さ方向に関する温度分布のシミュレーション結果を示す。物体４２の左端の位置座標が５ｍｍ、物体４１と４２との接触面の位置座標が２５ｍｍ、物体４１の右端の位置座標が４５ｍｍである。座標４５ｍｍの位置から熱エネルギが物体４１に流入し、座標５ｍｍの位置から熱エネルギが流出するため、温度勾配が負になっている。座標２５ｍｍの位置に、物体４１と４２との接触面の熱伝達係数に基づく温度の不連続が現れている。
【００６９】
図７に、物体４２（図６Ａ）の左端に位置する１つの粒子の温度の速度のシミュレーション結果を示す。横軸は、シミュレーションのタイムステップを表し、縦軸は、温度を単位「Ｋ」で表す。図７の破線は、一般的な有限要素法を用いてシミュレーションを行った結果を示し、実線は、実施例による方法を用いてシミュレーションを行った結果を示す。実線上に付された黒丸記号は、１回のタイムステップ終了時点における温度のシミュレーション結果を示す。
【００７０】
物体４２の左端の温度が定常状態の温度に到達するまでの期間、実施例によるシミュレーションで算出された温度の速度が、有限要素法を用いて算出された温度の速度よりも大きいことがわかる。このため、少ないタイムステップ数で熱伝導現象を定常状態に到達させることができる。これは、ステップ１３２（図２）において熱伝導方程式と等価な運動方程式の減衰係数を０に設定したことの効果である。
【００７１】
運動方程式の減衰係数を０に設定すると、通常はシミュレーション結果に振動波形が現れる。ところが、図７に示したシミュレーション結果には、振動波形が現れていない。これは、熱伝導現象が定常状態に到達しておらず、ステップ１６（図１）において条件が満たされると判定された場合に、ステップ１７（図１）において、各粒子の温度の速度を０に設定する処理（ＦＩＲＥ）を実行したことの効果である。
【００７２】
図７に示した例では、例えば温度の計算値が２５０Ｋ以上の領域では、ステップ１６（図１）において条件が満たされていないと判定されている。このため、温度の計算値が急峻に低下している。温度の計算値が２５０Ｋ以下になった最初のタイムステップで、ステップ１６（図１）において条件が満たされていると判定されている。このため、温度の計算値の低下測度が急激に小さくなっている。
【００７３】
温度の計算値が２４０Ｋ近傍に達すると、ステップ１３１（図２）において熱伝導現象が定常状態になったと判定されている。このため、熱伝導方程式と等価な運動方程式の減衰係数に本来の値が設定される（図２のステップ１３３）。これにより、物体４２の左端の粒子の本来の温度変化が、温度の計算値の変化に反映される。
【００７４】
図７に示したように、実施例によるシミュレーション方法を用いることにより、温度場が定常状態に到達するまでのタイムステップ数を少なくすることができる。実際に実施例による方法でシミュレーションを行ったところ、従来の有限要素法を用いる場合に比べて、２０倍以上高速に温度場の定常状態が得られることが分かった。
【００７５】
温度場をシミュレーションするときの好ましい時間刻み幅は、機構弾性現象をシミュレーションするときの好ましい時間刻み幅の１００００倍程度である。従って、温度場の解析と、機構弾性現象の解析とを連成させて、同一の時間刻み幅で実行する方法では、温度場が定常状態に到達するまでの数値計算のタイムステップ数が膨大になる。このため、温度場が定常状態に到達するまでシミュレーションを行うことは現実的には不可能であった。
【００７６】
上記実施例では、熱伝導方程式に等価な運動方程式の減衰係数を任意に設定することができるため、温度変化の時定数を機構弾性現象の時定数に近づけることが可能である。これにより、熱伝導現象の解析と、機構弾性現象の解析とを連成させて、熱伝導現象が定常状態に到達するまでシミュレーションを行うことが可能になる。
【００７７】
次に、他の実施例について説明する前に、図８〜図１０Ｂを参照して、参考例によるシミュレーション方法について説明する。以下、図１〜図７を参照して説明した実施例との相違点について説明し、共通の構成については説明を省略する。
【００７８】
図８に、参考例によるシミュレーション方法のフローチャートを示す。本実施例によるシミュレーション方法のステップ１０〜１２は、図１に示した実施例のステップ１０〜１２と同一である。本実施例によるシミュレーション方法のステップ１３が、図１に示した実施例のステップ１３と異なる。以下、ステップ１３の処理について説明する。
【００７９】
まず、ステップ１３５１において、温度場の各粒子が加速アルゴリズム適用条件を満たしているか否かを、粒子ごとに判定する。
【００８０】
以下、加速アルゴリズム適用条件について説明する。判定パラメータＰｉを、下記のように定義する。
【数２２】
0006679161.tif 000023
Ｆｉは、粒子ｉの温度Ｔｉの２次微分と仮想質量ｍｉとの積である。これは、粒子ｉの温度を変化させるための力に相当する。Ｆｉが正のとき、粒子ｉに、温度を上昇させる向きの力が作用していると考えられ、Ｆｉが負のとき、粒子ｉに、温度を低下させる向きの力が作用していると考えられる。Ｆｉの符号は、温度Ｔｉの加速度の符号と同一である。ｖｉは、粒子ｉの温度の速度に相当する。
【００８１】
判定パラメータＰｉが正のとき、加速アルゴリズム適用条件が満たされ、判定パラメータＰｉが負または０のとき、加速アルゴリズム適用条件が満たされていないと判定される。
【００８２】
次に、判定パラメータＰｉの物理的意味について説明する。判定パラメータＰｉが正であるということは、温度Ｔｉの加速度の符号と温度Ｔｉの速度の符号とが同一であることを意味する。すなわち、判定パラメータＰｉが正であるとき、温度Ｔｉが上昇しているとき、粒子ｉに、温度Ｔｉの上昇速度をより速める力が作用している。逆に、温度Ｔｉが低下しているとき、粒子ｉに、温度Ｔｉの下降速度をより速める力が作用している。この状態は、粒子ｉの温度Ｔｉが定常状態には至っておらず、定常状態に向かって変化している途上であると考えられる。
【００８３】
判定パラメータＰｉが負であるということは、温度Ｔｉの上昇速度を遅くする力、または温度Ｔｉの下降速度を遅くする力が粒子ｉに作用していることを意味する。従って、判定パラメータＰｉが負であるとき、粒子ｉの温度Ｔｉが定常状態に近づいた状態であると考えられる。
【００８４】
判定パラメータＰｉが正であるとき、加速アルゴリズム適用条件が満たされていると判定され、判定パラメータＰｉが負であるとき、加速アルゴリズム適用条件が満たされていないと判定される。判定パラメータＰｉが０の場合には、加速アルゴリズム適用条件が満たされていると判定してもよいし、満たされていないと判定してもよい。本実施例においては、判定パラメータＰｉが０のとき、加速アルゴリズム適用条件が満たされていないと判定される。
【００８５】
ステップ１３５１において加速アルゴリズム適用条件が満たされていると判定された場合、すなわち、温度Ｔｉの速度の符号と温度Ｔｉの加速度の符号とが同一である場合、ステップ１３５２において、熱伝導方程式と等価な運動方程式の減衰係数を現時点の値より小さくする。
【００８６】
具体的には、熱伝導方程式と等価な運動方程式（式（１４））の減衰係数（Ｔｉドットの係数）であるρｉＣｖｉΔＶｉを現時点の値より小さくする。一例として、減衰係数を、現時点の値の０．５倍とする。減衰係数を小さくすることにより、温度変化をより急峻にすることができる。
【００８７】
ステップ１３５１において加速アルゴリズム適用条件が満たされていないと判定された場合、すなわち、温度Ｔｉの速度の符号と温度Ｔｉの加速度の符号とが異なる場合、ステップ１３５３において、熱伝導方程式と等価な運動方程式（式（１４））の減衰係数（Ｔｉドットの係数）であるρｉＣｖｉΔＶｉを、本来の値に戻す。
【００８８】
ステップ１３５２及びステップ１３５３の後、ステップ１３５４において、修正後の減衰係数を用いて運動方程式（式（１４））を１タイムステップ実行する。その後、ステップ１４において、シミュレーションの終了条件が満たされているか否かを判定する。判定条件は、図１に示したステップ１４と同一とすることができる。シミュレーションの終了条件が満たされている場合には、シミュレーションを終了する。シミュレーションの終了条件が満たされていない場合には、ステップ１１に戻ってシミュレーションを継続する。
【００８９】
次に、図８に示したシミュレーション方法を用いて実際にシミュレーションを行った結果について説明する。簡易的なベアリングを動作させたときの温度変化のシミュレーションを行った。
【００９０】
図９に、シミュレーション対象である簡易的なベアリングの断面図を示す。シミュレーション対象のベアリングは、内輪５０、外輪５１、及び両者の間に配置された１２個のころ５２で構成される。外輪５１の外周面は、温度３００Ｋの外気に触れており、内輪５０を強制的に３６００ｒｐｍの回転速度で回転させた。外輪５１の外周面と外気との界面における熱伝達率は５，０００Ｗ／ｍ２／Ｋとした。内輪５０、外輪５１、及びころ５２の熱伝導率として、一般的なステンレス鋼の値を用いた。
【００９１】
粒子ｉと粒子ｊとの接触力をｆｉｊ、粒子ｉと粒子ｊとの相対速度をｖｉｊ、粒子ｉと粒子ｊとの間の摩擦係数をμｉｊで表したとき、粒子ｉの単位体積当たりの発熱量Ｑｉは、以下の式で表される。
【数２３】
0006679161.tif 000024
【００９２】
図１０Ａに、図９の参考例によるシミュレーション方法を用いて求めた全体温度の変化のグラフを示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。比較のために、減衰係数を本来の値に固定してシミュレーションを行った結果を、図１０Ｂに示す。ここで、「全体温度」は、温度場を構成する全粒子の温度の平均値を意味する。
【００９３】
図１０Ａと図１０Ｂとを比較すると、図１０Ａのシミュレーション結果の方が、定常状態に到達するまでの時間が著しく短いことが分かる。これは、図８のステップ１３５２において、減衰係数を小さくしたことの効果である。図１０Ｂのシミュレーション結果では、全体温度が３１９Ｋに収束しているが、図１０Ａのシミュレーション結果では、定常状態に到達した後も、幅１．０Ｋ程度の範囲内で全体温度が振動していることがわかる。
【００９４】
以下、定常状態に達した後も全体温度が振動する原因について考察する。図６Ａに示したシミュレーション対象においては、熱エネルギの流入及び流出が時間的に一定に保たれていた。これに対し、図９に示したシミュレーション対象においては、発熱場所及び発熱量が時間的に変化する。このため、定常状態に到達した後も、シミュレーション対象物の内部では、熱流が時間的空間的に変化する。その結果、個々の粒子に着目すると、ステップ１３５１（図８）の加速アルゴリズム適用条件が満たされなくなった後（定常状態に到達したと考えられた後）も、加速アルゴリズム適用条件が再度満たされてしまう場合が生じる。このように、加速アルゴリズム適用条件が満たされる状態と満たされない状態とが繰り返し発生するため、全体温度が振動すると考えられる。
【００９５】
次に、図１１及び図１２を参照して、全体温度の振動を無くすことが可能な実施例について説明する。
【００９６】
図１１に、本実施例によるシミュレーション方法のフローチャートを示す。以下、図８に示したフローチャートとの相違点について説明し、共通の手順については説明を省略する。
【００９７】
図１１に示した実施例では、ステップ１３５１の前にステップ１３６において、減衰係数を固定するか否かが判定される。減衰係数を固定すると判定された場合には、ステップ１３５３において、減衰係数を本来の値に戻す。すなわち、減衰係数を現時点の値より小さくするステップ１３５２は実行されない。一度、減衰係数が固定されると、その後は減衰係数が固定された状態が維持される。
【００９８】
減衰係数を固定するか否かの判定条件として、シミュレーション開始からの経過時間やタイムステップ数を用いることができる。この場合、経過時間やタイムステップ数が判定上限値を超えると、減衰係数を本来の値に固定すればよい。全体温度が定常状態に達していると判定された時点で、減衰係数を本来の値に固定してもよい。
【００９９】
図１２に、図１１に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。時刻ｔ２において減衰係数が固定される。減衰係数が固定されるまでは、全体温度が振動しているが、時刻ｔ２で減衰係数が固定された後は、この振動が発生していないことが分かる。
【０１００】
減衰係数が固定されるまでは、図８に示した参考例と同様にシミュレーションが行われるため、温度場が定常状態に達するまでのシミュレーション時間を短縮することができる。さらに、本実施例では、減衰係数を固定した後、全体温度の振動を防止することができる。
【０１０１】
次に、図１３及び図１４を参照して、さらに他の実施例について説明する。以下、図８に示した参考例との相違点について説明し、共通の手順については説明を省略する。
【０１０２】
図１３に、本実施例によるシミュレーション方法のフローチャートを示す。本実施例においては、ステップ１３５３が実行された後、またはその前に、さらに、ステップ１３５５が実行される。ステップ１３５５においては、熱運動方程式と等価な運動方程式の質量に相当するパラメータである仮想質量ｍｉを、現時点の仮想質量ｍｉより増加させる。
【０１０３】
ステップ１３５１で加速アルゴリズム適用条件が満たされていると判定されると、減衰係数は、ステップ１３５３で本来の値に戻された後でも、再度本来の値より小さくされる。これに対し、仮想質量ｍｉは、一旦増加されると、ステップ１３５１で加速アルゴリズム適用条件が満たされていると判定された場合でも元の値に戻されない。すなわち、仮想質量ｍｉは増加したままである。
【０１０４】
以下、ステップ１３５５の手順について、より具体的に説明する。ステップ１３５５では、仮想質量ｍｉを例えば現時点の値の１．０５倍にする。ただし、増加後の仮想質量ｍｉが、予め決められている上限値ｍｉｍａｘを超える場合には、仮想質量ｍｉを、上限値ｍｉｍａｘと等しくする。従って、仮想質量ｍｉが上限値ｍｉｍａｘを超えることはない。
【０１０５】
図１４に、図１３に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。図１４の太い実線ｄが、図１３に示した実施例によるシミュレーション方法を用いた場合の全体温度の計算結果を示し、細い実線ｅが、図８に示した参考例によるシミュレーション方法を用いた場合の全体温度の計算結果を示す。
【０１０６】
定常状態に至るまでの時間は、図８に示した参考例及び図１３に示した実施例の場合で、ほぼ同一である。図８に示した参考例の場合には、定常状態に達した後に、全体温度に振動が残っているのに対し、図１３に示した実施例の場合には、この振動が抑制されていることがわかる。
【０１０７】
図１３に示した実施例において、定常状態に達した後、仮想質量ｍｉが増加し、一旦増加した仮想質量ｍｉは減少しない。仮想質量ｍｉが増加すると、各粒子ｉの温度変動が緩やかになる。その結果、定常状態に達した後の全体温度の振動が抑制される。また、温度場が定常状態に達するまでは、ステップ１３５１（図１３）において加速アルゴリズム適用条件が満たされる頻度は低い。従って、温度場が定常状態に達する前に仮想質量ｍｉが大幅に増加することは回避される。このため、温度場が定常状態に達するまでは、急峻な温度変化が得られる。
【０１０８】
全体温度が振動することなく、一定の値に収束した後、仮想質量ｍｉと減衰係数とを本来の値に戻すことにより、図１０Ｂの場合と同様の終状態が得られる。
【０１０９】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【０１１０】
１０〜１７  ステップ
２０  中央処理ユニット
２１  記憶装置
２２  入力装置
２３  出力装置
２４  データバスライン
３０  中央領域
３１、３２  中間領域
３３、３４  端部領域
４１、４２  シミュレーション対象の物体
５０  内輪
５１  外輪
５２  コロ
１３１〜１３４、１３６  ステップ
１３５１〜１３５５  ステップ
【図１】
0006679161.tif 000025
【図２】
0006679161.tif 000026
【図３】
0006679161.tif 000027
【図４】
0006679161.tif 000028
【図５】
0006679161.tif 000029
【図６】
0006679161.tif 000030
【図７】
0006679161.tif 000031
【図８】
0006679161.tif 000032
【図９】
0006679161.tif 000033
【図１０】
0006679161.tif 000034
【図１１】
0006679161.tif 000035
【図１２】
0006679161.tif 000036
【図１３】
0006679161.tif 000037
【図１４】
0006679161.tif 000038