JP6679161B2 - シミュレーション方法、シミュレーション装置、及びシミュレーションプログラム - Google Patents

Description

本発明は、シミュレーション対象の熱伝導現象をシミュレーションするシミュレーション方法、シミュレーション装置、及びシミュレーションプログラムに関する。

古典力学や量子力学等を基に、コンピュータを用いて物質科学全般の現象を探るための方法として、分子動力学法（ＭＤ法）に基づくシミュレーション方法が知られている。マクロスケール系を扱うためにＭＤ法を発展させた繰り込み群分子動力学法（ＲＭＤ法）が提案されている。シミュレーション対象の熱伝導現象について、ＭＤ法やＲＭＤ法では、通常、格子振動（フォノン）による熱伝導しか取り扱うことができない。このため、熱伝導現象に自由電子が大きく寄与する金属の熱伝導現象をＭＤ法やＲＭＤ法でシミュレーションすると、実際の現象からかけ離れた結果が得られる場合が多い。

下記の特許文献１に、ＭＤ法やＲＭＤ法において熱伝導現象をより正確に扱うことができるシミュレーション方法が開示されている。特許文献１に開示されたシミュレーション方法においては、粒子に温度のパラメータを持たせる。ある時刻における各粒子の温度に基づいて熱伝導方程式を解くことにより、次の時刻における各粒子の温度を算出する。熱伝導方程式を解くことにより、実際のシミュレーション対象の比熱、熱伝導率を反映したシミュレーションが可能になる。

特許第５４４１４２２号公報

一般的に、熱伝導現象による温度変化は、機構的な現象や弾性現象（以下、「機構弾性現象」という。）による粒子位置の変化に比べて緩やかである。熱伝導現象をシミュレーションするときの時間刻み幅は、機構弾性現象をシミュレーションするときの時間刻み幅に比べて長く設定される。ところが、熱伝導現象と機構弾性現象との両方を取扱う連成シミュレーションにおいては、機構弾性現象をシミュレーションするための短い時間刻み幅で、熱伝導現象もシミュレーションされる。

温度場の解析においては、一般的に、定常状態における温度分布に興味が示される。温度分布が定常状態に達するまで、短い時間刻み幅で熱伝導現象をシミュレーションすると、計算時間が長大化する。

本発明の目的は、熱伝導現象を取扱う際の計算時間の長大化を抑制することが可能なシミュレーション方法を提供することである。本発明の他の目的は、計算時間の長大化を抑制することが可能なシミュレーション装置を提供することである。本発明のさらに他の目的は、計算時間の長大化を抑制することが可能なシミュレーションプログラムを提供することである。

本発明の一観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
前記運動方程式の数値計算において、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定するシミュレーション方法が提供される。

本発明の他の観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行うシミュレーション装置であって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式を記憶する記憶装置と、
前記記憶装置に記憶されている前記運動方程式の数値演算を行うことにより、前記シミュレーション対象物の温度分布を求める中央処理ユニットと、
前記中央処理ユニットで行われた数値演算の結果を出力する出力装置と
を有し、
前記運動方程式は減衰係数を含み、
前記中央処理ユニットは、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行うシミュレーション装置が提供される。

本発明のさらに他の観点によると、
複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションをコンピュータに実行させるシミュレーションプログラムであって、
空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行う機能、及び
前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行う機能を実現するシミュレーションプログラムが提供される。

熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、熱伝導現象のシミュレーションを行うことができる。運動方程式の数値計算において、種々の効率化、高速化手法を適用することができる。

図１は、実施例によるシミュレーション方法のフローチャートである。 図２は、図１に示したフローチャートのステップ１３のより詳細なフローチャートである。 図３は、シミュレーション対象物のある注目点における温度の速度のシミュレーション結果の一例を示すグラフである。 図４は、実施例によるシミュレーション方法を実行するシミュレーション装置のブロック図である。 図５は、熱伝導方程式と等価な運動方程式の減衰係数を、本来の値とは異なる値に設定しても、定常状態における正確な温度を求めることが可能なことを説明するためのモデルの斜視図である。 図６Ａは、シミュレーション対象物の斜視図であり、図６Ｂは、定常状態に達した時のシミュレーション対象物の長さ方向に関する温度分布のシミュレーション結果を示すグラフである。 図７は、シミュレーション対象物（図６Ａ）の左端に位置する１つの粒子の温度の速度のシミュレーション結果を示すグラフである。 図８は、参考例によるシミュレーション方法のフローチャートである。 図９は、シミュレーション対象である簡易的なベアリングの断面図である。 図１０Ａは、図９の参考例によるシミュレーション方法を用いて求めた全体温度の変化を示すグラフであり、図１０Ｂは、減衰係数を本来の値に固定してシミュレーションを行って求めた全体温度の変化を示すグラフである。 図１１は、他の実施例によるシミュレーション方法のフローチャートである。 図１２は、図１１に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示すグラフである。 図１３は、さらに他の実施例によるシミュレーション方法のフローチャートである。 図１４は、図１３に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示すグラフである。

図１に、実施例によるシミュレーション方法のフローチャートを示す。実施例においては、複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う。機構弾性現象のシミュレーションでは、シミュレーション対象物を粒子の集まりで表現し、ある時間刻み幅Δｔで各粒子の動解析を行う。この動解析には、例えば、ＭＤ(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)に説明されているように、下記の式で表される。

ここでσ_ｓ［φ（ｘ）］は拘束条件を表している。λ_ｓは、ラグランジェの未定定数である。δは、ディラックのδ関数である。δに付された指数３は、及びｄｘに付された指数３は、三次元を意味する。Σは、ｘに関して足し合わせることを意味する。ｘがｘ_ｓ以外のときは、δ関数の値は０である。式（１）において、ρ（ｘ）、κ、Ｑ（ｘ）は、特定の物理量を意味していない。

拘束条件σ_ｓ［φ（ｘ）］として、ディレクリ境界条件、ノイマン境界条件を採用することができる。ディレクリ境界条件は、下記の式で表すことができる。

式（２）は、ｘ＝ｘ_ｓにおいてφ（ｘ）が定数であることを意味する。

ノイマン境界条件は、下記の式で表すことができる。

式（３）は、ｘ＝ｘ_ｓにおいてφ（ｘ）の傾きが定数であることを意味する。φ（ｘ_ｓ）をφ_ｓと表している。∇_ｓは、ｘ＝ｘ_ｓの位置における傾きを表す。φ_ｓ０は、拘束条件に基づく定数である。

連続体のラグランジアンを定義する式（１）を、一次元の規則格子上で離散化し、離散化したスカラー場をφ_ｉと表記する。離散化したスカラー場φｉのラグランジアンＬは、以下の式で表される。

ここで、ΔＶ_ｉは線形体積素を表す。∇_ｉ±１は差分演算子を表す。差分演算子は、以下の式で定義される。

ここで、Δは、粒子間の距離を表す。

スカラー場φ_ｉを、粒子の属性、例えば位置、速度、温度等とみなせば、φ_ｉは以下のラグランジェの方程式に従う。

ここで、γ_ｉは減衰係数である。式（４）、式（６）から、ｉ番目の粒子に対する以下の運動方程式が得られる。

ここで、式（７）の右辺第１項は弾性項を記述し、第２項は減衰項を記述し、第３項は外力を記述し、第４項は拘束条件を記述している。右辺の第１項のΣのｊ＝±１は、ｉ番目の粒子の両側の粒子について和を取ることを意味する。右辺の第４項のΣのｓは、拘束条件の数に相当する。

（ｄ^２／ｄｔ^２）φ_ｉが０になる定常状態において、式（７）は、以下のように変形できる。

式（８）の右辺第１項は、κ∇^２φ（ｘ）の差分近似を表しているから、式（８）は拡散方程式に帰着する。式（８）は、式（７）の左辺を０と仮定して得られたものである。従って、式（７）の左辺のρ_ｉを十分小さすれば、式（６）の解が、拡散方程式（８）の解に漸近することがわかる。すなわち、定微分方程式の数値解法を用いて式（６）を解くことにより、近似的に、式（８）の拡散方程式を解くことができる。定微分方程式の数値解法として、ベレル(Verlet)法、リープフロッグ法、ギア(Gear)法等が挙げられる。

熱伝導方程式は、以下の式で表される。

ここで、ρは密度、Ｃ_ｖは比熱、Ｔは温度、ｔは時間、Ｋは熱伝導率、Ｑは単位体積当たりの発熱量を表す。

式（８）のφ_ｉ、γ_ｉ、κ_ｉ、及びＱ_ｉを、それぞれ温度Ｔ、密度ρと比熱Ｃ_ｖとの積、熱伝導率Ｋ、及び発熱量Ｑとみなせば、式（８）が式（９）の熱伝導方程式と同形であることがわかる。具体的には、式（８）の空間的なφ_ｉの分布の項、及びφ_ｉの時間微分の項を、それぞれ式（９）の空間的な温度分布の項、及び温度の時間微分の項に対応付けることができる。このため、ρ_ｉが十分小さいという条件の下で、式（７）の運動方程式が式（９）の熱伝導方程式と等価であるといえる。

ここまでは、古典的スカラー場φが一次元の規則格子上で離散化された場合について説明した。以上の議論を、三次元の複雑形状に拡張するためには、式（７）に有限体積法を適用すればよい。式（７）に有限体積法を適用すると、以下の式が得られる。

ここで、式（１０）の右辺第１項のΣは、ｉ番目の粒子に接触している全ての粒子について足し合わせることを意味する。ΔＳ_ｉｊはボロノイ多面体ｉとｊとの接触面の面積を表す。ｍ_ｉは下記の式で表される。

ｍ_ｉは、運動方程式の質量に相当するパラメータであるため、ｍ_ｉを仮想質量と呼ぶこととする。

仮想質量ｍ_ｉは、数値積分の安定条件から、以下の下限条件を満足する必要がある。

時間刻み幅ｄｔについては、機構弾性現象の解析を行う際の粒子間距離（メッシュの寸法）及び音速に基づいて、最適値を決定することできる。機構弾性現象に基づいて決定された時間刻み幅ｄｔ、及び式（１２）から、仮想質量ｍ_ｉの下限値が決定される。さらに、上述のように、ρ_ｉが十分小さいという条件の下で、式（６）の解が拡散方程式（８）の解に漸近する。このため、ρ_ｉは、シミュレーション対象の真の密度ρよりも十分小さいことが好ましい。すなわち、以下の式が満たされることが好ましい。

式（１０）のγ_ｉを、密度ρ_ｉと比熱Ｃ_ｖｉとの積に置き換え、φ_ｉを温度Ｔ_ｉに置き換えることにより、以下の式が得られる。Ｃ_ｖｉは、ｉ番目の粒子の位置の比熱を表す。

運動方程式と同形の式（１４）を解くことにより、温度場のシミュレーションを行うことができる。式（１０）のφ_ｉを式（１４）においてＴ_ｉに置き換えたため、式（１４）のｍ_ｉの次元は、式（１０）、（１１）のｍ_ｉの次元とは異なる。仮想質量ｍ_ｉの次元は、式（１４）の右辺のΔＶ_ｉを含む項から、密度ρ_ｉ、比熱Ｃ_ｖ、体積素ΔＶ_ｉ、及び時間刻み幅ｄｔの積の次元と同一であることがわかる。すなわち、仮想質量ｍｉの次元は、［Ｊ／Ｋ］［ｓ］となる。

式（１２）から、仮想質量ｍ_ｉの満たすべき条件は、以下の式で与えられる。

さらに、仮想質量ｍ_ｉの満たすべき条件は、種々の数値計算を行った結果、経験的に以下の式で与えられる。

式（１５）及び式（１６）から、一例として、ｍ_ｉを下記の式のように設定することが推奨される。

以上、説明したように、熱伝導方程式（９）と等価な運動方程式（１４）を解くことにより温度場の解析を行うことができる。

ステップ１３（図１）において、運動方程式（７）の減衰係数γ_ｉ、式（１４）においてはρ_ｉＣ_ｖｉを０にすれば、温度Ｔの速度の算出結果が大きくなる。すなわち、温度Ｔの時間変化が急峻になる。このため、温度の計算値が定常状態における温度に達するまでの時間を短くすることができる。

次に、ステップ１３２（図２）において、減衰係数γ_ｉまたはρ_ｉＣ_ｖｉを、本来の値とは異なる値、または０に設定しても、定常状態における正確な温度を求めることが可能なことを示す。

図５に示した棒状のモデルを考える。体積ｄＶ_ｉの中央領域３０に、単位体積あたりの発熱量Ｑ_ｉの熱源が存在する。中央領域３０から、面積Ｓ_ｉの界面を通って、両側の中間領域３１、３２に向って熱伝導が生じる。中間領域３１、３２の端に、それぞれ端部領域３３、３４が連続している。端部領域３３、３４が、それぞれ面積Ｓ_ｊの端面において温度Ｔ_０の熱浴に接している。

中央領域３０から一方の中間領域３１に流れる単位面積、単位時間あたりの熱エネルギｊ_Ａは、以下の式で表される。

端部領域３３の温度をＴで表し、端部領域３３と熱浴との熱伝達係数をｈで表すと、端部領域３３から熱浴に流れる単位面積、単位時間あたりの熱エネルギｊ_Ｂは、以下の式で表される。

中央領域３０から中間領域３１に流れる単位時間あたりの熱エネルギと、端部領域３３から熱浴に流れる単位時間あたりの熱エネルギとは等しい。このため、以下の式が成立する。

式（１８）〜式（２０）から、以下の式が導出される。

式（２１）から、定常状態における表面温度は、熱源からの発熱量と、シミュレーション対象物の表面積のみで決まることがわかる。定常状態における表面温度は、密度及び比熱に依存しない。従って、式（１４）の減衰係数ρ_ｉＣ_ｖｉを任意の値、または０として運動方程式（１４）を解いても、定常状態における表面温度を求めることが可能である。

減衰係数ρ_ｉＣ_ｖｉを小さくすると、図３の破線ｂで示したように、温度の計算値に振動が現れる。減衰係数ρ_ｉＣ_ｖｉの大きさに応じて振動が減衰し、温度の計算値が定常状態の温度に収束する。

減衰係数ρ_ｉＣ_ｖｉを０にすると振動が減衰せず、いつまでも振動が継続してしまう。実施例においては、ＦＩＲＥを適用することにより、図３の実線ｃで示したように、温度の計算値の振動を短時間で収束させている。

温度場が定常状態に達するまでの期間、すなわち減衰係数を本来の値と異なる値に設定してステップ１３を実行している期間は、熱伝導方程式と等価な運動方程式を解くときの時間刻み幅の実質的な意味が、ステップ１１（図１）における動解析で適用される時間刻み幅と異なる。このため、シミュレーション対象物の機構弾性現象の時間発展と、熱伝導現象の時間発展とが、非同期で進行することになる。温度場が定常状態に達した後、減衰係数として本来の値を設定すると、シミュレーション対象の機構弾性現象の時間発展と、熱伝導現象の時間発展とが同期することになる。

上記実施例によるシミュレーション方法の効果を確認するために、物体の温度場のシミュレーションを行った。次に、図６Ａ、図６Ｂ、及び図７を参照して、このシミュレーション結果について説明する。

図６Ａに、シミュレーション対象物の斜視図を示す。四角柱状の物体４１の左側の端面と、四角柱状の物体４２の右側の端面とが熱接触している。物体４１及び物体４２の長さは共に２０ｍｍであり、長さ方向に直交する断面は正方形であり、面積は同一である。

物体４１の右端から物体４１に熱エネルギが流入する。物体４１の左端と物体４２の右端との接触面を介して、物体４１から物体４２に熱が伝達される。物体４２の左端から、熱エネルギが流出する。物体４１に流入する単位時間あたりの熱エネルギと、物体４２の左端から流出する単位時間あたりの熱エネルギとは等しい。このため、物体４１と物体４２との熱収支が０に保たれる。物体４１及び４２の温度場は、最終的に定常状態に達する。すなわち、温度の速度が０になる。

図６Ｂに、定常状態に達した時の物体４１及び４２の長さ方向に関する温度分布のシミュレーション結果を示す。物体４２の左端の位置座標が５ｍｍ、物体４１と４２との接触面の位置座標が２５ｍｍ、物体４１の右端の位置座標が４５ｍｍである。座標４５ｍｍの位置から熱エネルギが物体４１に流入し、座標５ｍｍの位置から熱エネルギが流出するため、温度勾配が負になっている。座標２５ｍｍの位置に、物体４１と４２との接触面の熱伝達係数に基づく温度の不連続が現れている。

図７に、物体４２（図６Ａ）の左端に位置する１つの粒子の温度の速度のシミュレーション結果を示す。横軸は、シミュレーションのタイムステップを表し、縦軸は、温度を単位「Ｋ」で表す。図７の破線は、一般的な有限要素法を用いてシミュレーションを行った結果を示し、実線は、実施例による方法を用いてシミュレーションを行った結果を示す。実線上に付された黒丸記号は、１回のタイムステップ終了時点における温度のシミュレーション結果を示す。

物体４２の左端の温度が定常状態の温度に到達するまでの期間、実施例によるシミュレーションで算出された温度の速度が、有限要素法を用いて算出された温度の速度よりも大きいことがわかる。このため、少ないタイムステップ数で熱伝導現象を定常状態に到達させることができる。これは、ステップ１３２（図２）において熱伝導方程式と等価な運動方程式の減衰係数を０に設定したことの効果である。

運動方程式の減衰係数を０に設定すると、通常はシミュレーション結果に振動波形が現れる。ところが、図７に示したシミュレーション結果には、振動波形が現れていない。これは、熱伝導現象が定常状態に到達しておらず、ステップ１６（図１）において条件が満たされると判定された場合に、ステップ１７（図１）において、各粒子の温度の速度を０に設定する処理（ＦＩＲＥ）を実行したことの効果である。

図７に示した例では、例えば温度の計算値が２５０Ｋ以上の領域では、ステップ１６（図１）において条件が満たされていないと判定されている。このため、温度の計算値が急峻に低下している。温度の計算値が２５０Ｋ以下になった最初のタイムステップで、ステップ１６（図１）において条件が満たされていると判定されている。このため、温度の計算値の低下測度が急激に小さくなっている。

温度の計算値が２４０Ｋ近傍に達すると、ステップ１３１（図２）において熱伝導現象が定常状態になったと判定されている。このため、熱伝導方程式と等価な運動方程式の減衰係数に本来の値が設定される（図２のステップ１３３）。これにより、物体４２の左端の粒子の本来の温度変化が、温度の計算値の変化に反映される。

図７に示したように、実施例によるシミュレーション方法を用いることにより、温度場が定常状態に到達するまでのタイムステップ数を少なくすることができる。実際に実施例による方法でシミュレーションを行ったところ、従来の有限要素法を用いる場合に比べて、２０倍以上高速に温度場の定常状態が得られることが分かった。

温度場をシミュレーションするときの好ましい時間刻み幅は、機構弾性現象をシミュレーションするときの好ましい時間刻み幅の１００００倍程度である。従って、温度場の解析と、機構弾性現象の解析とを連成させて、同一の時間刻み幅で実行する方法では、温度場が定常状態に到達するまでの数値計算のタイムステップ数が膨大になる。このため、温度場が定常状態に到達するまでシミュレーションを行うことは現実的には不可能であった。

上記実施例では、熱伝導方程式に等価な運動方程式の減衰係数を任意に設定することができるため、温度変化の時定数を機構弾性現象の時定数に近づけることが可能である。これにより、熱伝導現象の解析と、機構弾性現象の解析とを連成させて、熱伝導現象が定常状態に到達するまでシミュレーションを行うことが可能になる。

次に、他の実施例について説明する前に、図８〜図１０Ｂを参照して、参考例によるシミュレーション方法について説明する。以下、図１〜図７を参照して説明した実施例との相違点について説明し、共通の構成については説明を省略する。

図８に、参考例によるシミュレーション方法のフローチャートを示す。本実施例によるシミュレーション方法のステップ１０〜１２は、図１に示した実施例のステップ１０〜１２と同一である。本実施例によるシミュレーション方法のステップ１３が、図１に示した実施例のステップ１３と異なる。以下、ステップ１３の処理について説明する。

まず、ステップ１３５１において、温度場の各粒子が加速アルゴリズム適用条件を満たしているか否かを、粒子ごとに判定する。

以下、加速アルゴリズム適用条件について説明する。判定パラメータＰ_ｉを、下記のように定義する。

Ｆ_ｉは、粒子ｉの温度Ｔ_ｉの２次微分と仮想質量ｍ_ｉとの積である。これは、粒子ｉの温度を変化させるための力に相当する。Ｆ_ｉが正のとき、粒子ｉに、温度を上昇させる向きの力が作用していると考えられ、Ｆ_ｉが負のとき、粒子ｉに、温度を低下させる向きの力が作用していると考えられる。Ｆ_ｉの符号は、温度Ｔ_ｉの加速度の符号と同一である。ｖ_ｉは、粒子ｉの温度の速度に相当する。

判定パラメータＰ_ｉが正のとき、加速アルゴリズム適用条件が満たされ、判定パラメータＰ_ｉが負または０のとき、加速アルゴリズム適用条件が満たされていないと判定される。

次に、判定パラメータＰ_ｉの物理的意味について説明する。判定パラメータＰ_ｉが正であるということは、温度Ｔ_ｉの加速度の符号と温度Ｔ_ｉの速度の符号とが同一であることを意味する。すなわち、判定パラメータＰ_ｉが正であるとき、温度Ｔ_ｉが上昇しているとき、粒子ｉに、温度Ｔ_ｉの上昇速度をより速める力が作用している。逆に、温度Ｔ_ｉが低下しているとき、粒子ｉに、温度Ｔ_ｉの下降速度をより速める力が作用している。この状態は、粒子ｉの温度Ｔ_ｉが定常状態には至っておらず、定常状態に向かって変化している途上であると考えられる。

判定パラメータＰ_ｉが負であるということは、温度Ｔ_ｉの上昇速度を遅くする力、または温度Ｔ_ｉの下降速度を遅くする力が粒子ｉに作用していることを意味する。従って、判定パラメータＰ_ｉが負であるとき、粒子ｉの温度Ｔ_ｉが定常状態に近づいた状態であると考えられる。

判定パラメータＰ_ｉが正であるとき、加速アルゴリズム適用条件が満たされていると判定され、判定パラメータＰ_ｉが負であるとき、加速アルゴリズム適用条件が満たされていないと判定される。判定パラメータＰ_ｉが０の場合には、加速アルゴリズム適用条件が満たされていると判定してもよいし、満たされていないと判定してもよい。本実施例においては、判定パラメータＰ_ｉが０のとき、加速アルゴリズム適用条件が満たされていないと判定される。

ステップ１３５１において加速アルゴリズム適用条件が満たされていると判定された場合、すなわち、温度Ｔ_ｉの速度の符号と温度Ｔ_ｉの加速度の符号とが同一である場合、ステップ１３５２において、熱伝導方程式と等価な運動方程式の減衰係数を現時点の値より小さくする。

具体的には、熱伝導方程式と等価な運動方程式（式（１４））の減衰係数（Ｔｉドットの係数）であるρ_ｉＣ_ｖｉΔＶ_ｉを現時点の値より小さくする。一例として、減衰係数を、現時点の値の０．５倍とする。減衰係数を小さくすることにより、温度変化をより急峻にすることができる。

ステップ１３５１において加速アルゴリズム適用条件が満たされていないと判定された場合、すなわち、温度Ｔ_ｉの速度の符号と温度Ｔ_ｉの加速度の符号とが異なる場合、ステップ１３５３において、熱伝導方程式と等価な運動方程式（式（１４））の減衰係数（Ｔｉドットの係数）であるρ_ｉＣ_ｖｉΔＶ_ｉを、本来の値に戻す。

ステップ１３５２及びステップ１３５３の後、ステップ１３５４において、修正後の減衰係数を用いて運動方程式（式（１４））を１タイムステップ実行する。その後、ステップ１４において、シミュレーションの終了条件が満たされているか否かを判定する。判定条件は、図１に示したステップ１４と同一とすることができる。シミュレーションの終了条件が満たされている場合には、シミュレーションを終了する。シミュレーションの終了条件が満たされていない場合には、ステップ１１に戻ってシミュレーションを継続する。

次に、図８に示したシミュレーション方法を用いて実際にシミュレーションを行った結果について説明する。簡易的なベアリングを動作させたときの温度変化のシミュレーションを行った。

図９に、シミュレーション対象である簡易的なベアリングの断面図を示す。シミュレーション対象のベアリングは、内輪５０、外輪５１、及び両者の間に配置された１２個のころ５２で構成される。外輪５１の外周面は、温度３００Ｋの外気に触れており、内輪５０を強制的に３６００ｒｐｍの回転速度で回転させた。外輪５１の外周面と外気との界面における熱伝達率は５，０００Ｗ／ｍ^２／Ｋとした。内輪５０、外輪５１、及びころ５２の熱伝導率として、一般的なステンレス鋼の値を用いた。

粒子ｉと粒子ｊとの接触力をｆ_ｉｊ、粒子ｉと粒子ｊとの相対速度をｖ_ｉｊ、粒子ｉと粒子ｊとの間の摩擦係数をμ_ｉｊで表したとき、粒子ｉの単位体積当たりの発熱量Ｑｉは、以下の式で表される。

図１０Ａに、図９の参考例によるシミュレーション方法を用いて求めた全体温度の変化のグラフを示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。比較のために、減衰係数を本来の値に固定してシミュレーションを行った結果を、図１０Ｂに示す。ここで、「全体温度」は、温度場を構成する全粒子の温度の平均値を意味する。

図１０Ａと図１０Ｂとを比較すると、図１０Ａのシミュレーション結果の方が、定常状態に到達するまでの時間が著しく短いことが分かる。これは、図８のステップ１３５２において、減衰係数を小さくしたことの効果である。図１０Ｂのシミュレーション結果では、全体温度が３１９Ｋに収束しているが、図１０Ａのシミュレーション結果では、定常状態に到達した後も、幅１．０Ｋ程度の範囲内で全体温度が振動していることがわかる。

以下、定常状態に達した後も全体温度が振動する原因について考察する。図６Ａに示したシミュレーション対象においては、熱エネルギの流入及び流出が時間的に一定に保たれていた。これに対し、図９に示したシミュレーション対象においては、発熱場所及び発熱量が時間的に変化する。このため、定常状態に到達した後も、シミュレーション対象物の内部では、熱流が時間的空間的に変化する。その結果、個々の粒子に着目すると、ステップ１３５１（図８）の加速アルゴリズム適用条件が満たされなくなった後（定常状態に到達したと考えられた後）も、加速アルゴリズム適用条件が再度満たされてしまう場合が生じる。このように、加速アルゴリズム適用条件が満たされる状態と満たされない状態とが繰り返し発生するため、全体温度が振動すると考えられる。

次に、図１１及び図１２を参照して、全体温度の振動を無くすことが可能な実施例について説明する。

図１１に、本実施例によるシミュレーション方法のフローチャートを示す。以下、図８に示したフローチャートとの相違点について説明し、共通の手順については説明を省略する。

図１１に示した実施例では、ステップ１３５１の前にステップ１３６において、減衰係数を固定するか否かが判定される。減衰係数を固定すると判定された場合には、ステップ１３５３において、減衰係数を本来の値に戻す。すなわち、減衰係数を現時点の値より小さくするステップ１３５２は実行されない。一度、減衰係数が固定されると、その後は減衰係数が固定された状態が維持される。

減衰係数を固定するか否かの判定条件として、シミュレーション開始からの経過時間やタイムステップ数を用いることができる。この場合、経過時間やタイムステップ数が判定上限値を超えると、減衰係数を本来の値に固定すればよい。全体温度が定常状態に達していると判定された時点で、減衰係数を本来の値に固定してもよい。

図１２に、図１１に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。時刻ｔ_２において減衰係数が固定される。減衰係数が固定されるまでは、全体温度が振動しているが、時刻ｔ_２で減衰係数が固定された後は、この振動が発生していないことが分かる。

減衰係数が固定されるまでは、図８に示した参考例と同様にシミュレーションが行われるため、温度場が定常状態に達するまでのシミュレーション時間を短縮することができる。さらに、本実施例では、減衰係数を固定した後、全体温度の振動を防止することができる。

次に、図１３及び図１４を参照して、さらに他の実施例について説明する。以下、図８に示した参考例との相違点について説明し、共通の手順については説明を省略する。

図１３に、本実施例によるシミュレーション方法のフローチャートを示す。本実施例においては、ステップ１３５３が実行された後、またはその前に、さらに、ステップ１３５５が実行される。ステップ１３５５においては、熱運動方程式と等価な運動方程式の質量に相当するパラメータである仮想質量ｍ_ｉを、現時点の仮想質量ｍ_ｉより増加させる。

ステップ１３５１で加速アルゴリズム適用条件が満たされていると判定されると、減衰係数は、ステップ１３５３で本来の値に戻された後でも、再度本来の値より小さくされる。これに対し、仮想質量ｍｉは、一旦増加されると、ステップ１３５１で加速アルゴリズム適用条件が満たされていると判定された場合でも元の値に戻されない。すなわち、仮想質量ｍｉは増加したままである。

以下、ステップ１３５５の手順について、より具体的に説明する。ステップ１３５５では、仮想質量ｍｉを例えば現時点の値の１．０５倍にする。ただし、増加後の仮想質量ｍｉが、予め決められている上限値ｍ_ｉｍａｘを超える場合には、仮想質量ｍ_ｉを、上限値ｍ_ｉｍａｘと等しくする。従って、仮想質量ｍ_ｉが上限値ｍ_ｉｍａｘを超えることはない。

図１４に、図１３に示した実施例によるシミュレーション方法を用いて図９に示したベアリングの温度変化を求めた結果を示す。横軸は経過時間を単位「ｓ」で表し、縦軸は温度を単位「Ｋ」で表す。図１４の太い実線ｄが、図１３に示した実施例によるシミュレーション方法を用いた場合の全体温度の計算結果を示し、細い実線ｅが、図８に示した参考例によるシミュレーション方法を用いた場合の全体温度の計算結果を示す。

定常状態に至るまでの時間は、図８に示した参考例及び図１３に示した実施例の場合で、ほぼ同一である。図８に示した参考例の場合には、定常状態に達した後に、全体温度に振動が残っているのに対し、図１３に示した実施例の場合には、この振動が抑制されていることがわかる。

図１３に示した実施例において、定常状態に達した後、仮想質量ｍ_ｉが増加し、一旦増加した仮想質量ｍ_ｉは減少しない。仮想質量ｍ_ｉが増加すると、各粒子ｉの温度変動が緩やかになる。その結果、定常状態に達した後の全体温度の振動が抑制される。また、温度場が定常状態に達するまでは、ステップ１３５１（図１３）において加速アルゴリズム適用条件が満たされる頻度は低い。従って、温度場が定常状態に達する前に仮想質量ｍ_ｉが大幅に増加することは回避される。このため、温度場が定常状態に達するまでは、急峻な温度変化が得られる。

全体温度が振動することなく、一定の値に収束した後、仮想質量ｍ_ｉと減衰係数とを本来の値に戻すことにより、図１０Ｂの場合と同様の終状態が得られる。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

１０〜１７ ステップ
２０ 中央処理ユニット
２１ 記憶装置
２２ 入力装置
２３ 出力装置
２４ データバスライン
３０ 中央領域
３１、３２ 中間領域
３３、３４ 端部領域
４１、４２ シミュレーション対象の物体
５０ 内輪
５１ 外輪
５２ コロ
１３１〜１３４、１３６ ステップ
１３５１〜１３５５ ステップ

Claims (10)

1. 複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
   空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
   前記運動方程式の数値計算において、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定するシミュレーション方法。
2. 前記運動方程式の数値計算のタイムステップごとに、温度の速度を０に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を０に設定する請求項１に記載のシミュレーション方法。
3. 前記運動方程式の数値計算のタイムステップごとに、前記熱伝導現象のシミュレーションで求められた温度の計算値の分布が定常状態に到達したか否かを判定する工程と、
   前記熱伝導現象のシミュレーションで求められた温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う工程と
   を有する請求項１または２に記載のシミュレーション方法。
4. 複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
   空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
   前記運動方程式の数値計算を行う際に、
   前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、
   前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻し、
   シミュレーション開始からのタイムステップ数が基準値以上になったら、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に固定するシミュレーション方法。
5. 複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、
   空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行い、
   前記運動方程式の数値計算を行う際に、
   前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、
   前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻すとともに、前記運動方程式の質量に相当するパラメータである仮想質量を、現時点の仮想質量より増加させるシミュレーション方法。
6. 前記仮想質量が上限値を超えると、前記仮想質量を前記上限値に設定する請求項５に記載のシミュレーション方法。
7. 複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行うシミュレーション装置であって、
   空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式を記憶する記憶装置と、
   前記記憶装置に記憶されている前記運動方程式の数値演算を行うことにより、前記シミュレーション対象物の温度分布を求める中央処理ユニットと、
   前記中央処理ユニットで行われた数値演算の結果を出力する出力装置と
   を有し、
   前記運動方程式は減衰係数を含み、
   前記中央処理ユニットは、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行うシミュレーション装置。
8. 前記中央処理ユニットは、前記運動方程式の数値計算のタイムステップごとに、温度の速度を０に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を０に設定する請求項７に記載のシミュレーション装置。
9. 前記中央処理ユニットは、
   前記運動方程式の数値計算のタイムステップごとに、前記運動方程式の数値演算で算出された温度の計算値の分布が定常状態に到達したか否かを判定する工程と、
   前記運動方程式の数値演算で算出された温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う請求項７または８に記載のシミュレーション装置。
10. 複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションをコンピュータに実行させるシミュレーションプログラムであって、
    空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行う機能、及び
    前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または０を設定して前記運動方程式の数値計算を行う機能を実現するシミュレーションプログラム。
    

Images