図面 (/)

技術 シミュレーション方法、シミュレーション装置、及びシミュレーションプログラム

出願人 住友重機械工業株式会社
発明者 大西良孝市嶋大路宮崎修司
出願日 2016年2月2日 (4年8ヶ月経過) 出願番号 2016-017588
公開日 2017年3月9日 (3年7ヶ月経過) 公開番号 2017-049971
状態 特許登録済
技術分野 特定用途計算機
主要キーワード 線形体 アルゴリズム適用 中央領 物質科学 差分近似 シミュレーション対象物 古典力学 スカラー場
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(2017年3月9日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (15)

課題

機構弾性現象の動解析連成させて熱伝導現象のシミュレーションを行う際の計算時間の長大化を抑制することが可能なシミュレーション方法を提供する。

解決手段

複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う。この際に、空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式同形の式に変形可能な運動方程式数値計算を行うことにより、シミュレーション対象物の熱伝導現象のシミュレーションを行う。

概要

背景

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

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

概要

機構弾性現象の動解析連成させて熱伝導現象のシミュレーションを行う際の計算時間の長大化を抑制することが可能なシミュレーション方法を提供する。複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う。この際に、空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式数値計算を行うことにより、シミュレーション対象物の熱伝導現象のシミュレーションを行う。

目的

本発明の目的は、熱伝導現象を取扱う際の計算時間の長大化を抑制することが可能なシミュレーション方法を提供する

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

(分野番号表示ON)※整理標準化データをもとに当社作成

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

複数の粒子を含むシミュレーション対象物機構弾性現象と、熱伝導現象との連成シミュレーションを行う方法であって、空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式同形の式に変形可能な運動方程式数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行うシミュレーション方法

請求項2

前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または0を設定して前記運動方程式の数値計算を行う請求項1に記載のシミュレーション方法。

請求項3

前記運動方程式の数値計算のタイムステップごとに、温度の速度を0に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を0に設定する請求項2に記載のシミュレーション方法。

請求項4

前記運動方程式の数値計算のタイムステップごとに、前記熱伝導現象のシミュレーションで求められた温度の計算値分布定常状態に到達したか否かを判定する工程と、前記熱伝導現象のシミュレーションで求められた温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う工程とを有する請求項2または3に記載のシミュレーション方法。

請求項5

前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻し、シミュレーション開始からのタイムステップ数が基準値以上になったら、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に固定する請求項1に記載のシミュレーション方法。

請求項6

前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが同一である場合、前記運動方程式の減衰係数を現時点の減衰係数より小さくし、前記シミュレーション対象物の温度の速度の符号と、温度の加速度の符号とが異なる場合、前記運動方程式の減衰係数を、前記熱伝導方程式から求まる本来の値に戻すとともに、前記運動方程式の質量に相当するパラメータである仮想質量を、現時点の仮想質量より増加させる請求項1に記載のシミュレーション方法。

請求項7

前記仮想質量が上限値を超えると、前記仮想質量を前記上限値に設定する請求項6に記載のシミュレーション方法。

請求項8

複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行うシミュレーション装置であって、空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式を記憶する記憶装置と、前記記憶装置に記憶されている前記運動方程式の数値演算を行うことにより、前記シミュレーション対象物の温度分布を求める中央処理ユニットと、前記中央処理ユニットで行われた数値演算の結果を出力する出力装置とを有するシミュレーション装置。

請求項9

前記運動方程式は減衰係数を含み、前記中央処理ユニットは、前記運動方程式の減衰係数に、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数よりも小さい値、または0を設定して前記運動方程式の数値計算を行う請求項8に記載のシミュレーション装置。

請求項10

前記中央処理ユニットは、前記運動方程式の数値計算のタイムステップごとに、温度の速度を0に設定する条件が満たされるか否かを判定し、前記条件が満たされる場合には、前記運動方程式の数値計算の次のタイムステップにおいて、温度の速度を0に設定する請求項9に記載のシミュレーション装置。

請求項11

前記中央処理ユニットは、前記運動方程式の数値計算のタイムステップごとに、前記運動方程式の数値演算で算出された温度の計算値の分布が定常状態に到達したか否かを判定する工程と、前記運動方程式の数値演算で算出された温度の計算値の分布が前記定常状態に到達したと判定された場合には、前記運動方程式の減衰係数を、前記シミュレーション対象物に適用される前記熱伝導方程式から求まる本来の減衰係数に設定して、前記運動方程式の数値計算を行う請求項9または10に記載のシミュレーション装置。

請求項12

複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションをコンピュータに実行させるシミュレーションプログラムであって、空間的な温度分布の項、及び温度の時間微分の項に関して、熱伝導方程式と同形の式に変形可能な運動方程式の数値計算を行うことにより、前記シミュレーション対象物の熱伝導現象のシミュレーションを行う機能を実現するシミュレーションプログラム。

技術分野

背景技術

0002

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

0003

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

先行技術

0004

特許第5441422号公報

発明が解決しようとする課題

0005

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

0006

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

0007

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

課題を解決するための手段

0008

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

0009

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

0010

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

発明の効果

0011

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

図面の簡単な説明

0012

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

実施例

0013

図1に、実施例によるシミュレーション方法のフローチャートを示す。実施例においては、複数の粒子を含むシミュレーション対象物の機構弾性現象と、熱伝導現象との連成シミュレーションを行う。機構弾性現象のシミュレーションでは、シミュレーション対象物を粒子の集まり表現し、ある時間刻み幅Δtで各粒子の動解析を行う。この動解析には、例えば、MD(molecular dynamics)法、SPH(Smoothed Particle Hydrodynamics)法、MPS(Moving Particle Simulation)法、DEM(Distinct Element Method)等の粒子法を適用可能である。各粒子をメッシュ節点とみなせば、有限要素法(FEM)を適用することも可能である。

0014

熱伝導現象のシミュレーションでは、各粒子に温度のパラメータを持たせ、動解析と同一の時間刻み幅Δtで温度の数値計算を行う。温度の数値計算を行う際に、熱伝導方程式を解く代わりに、熱伝導方程式と等価な運動方程式を解くことにより、温度場の解析を行う。ここで、「等価」とは、熱伝導方程式の空間的な温度分布の項と、温度の時間微分の項とに関して、熱伝導方程式と同形の式に変形可能であることを意味する。運動方程式における質点の位置を、「温度」に対応させると、温度の時間変化の傾き(温度の一次微分)が運動方程式における速度に対応する。本明細書において、温度の時間変化の傾きを「温度の速度」ということとする。

0015

温度場の解析において、温度変化の時定数可変とすることにより、温度変化の傾きを実際の温度変化の傾きより急峻にする。これにより、温度の計算結果が定常状態に到達するまでの時間を短くすることができる。ここで、「定常状態」とは、各粒子の位置における温度の速度が0であることを意味する。

0016

シミュレーションにおいて温度変化の傾きを急峻にすると、シミュレーションで求められた温度がオーバシュートして、振動波形が現れる場合がある。実施例では、この振動短期緩和させる手法を適用することにより、温度場が定常状態に達するまでの時間の短縮化を図っている。振動を短期に緩和させる手法として、例えばFIRE(fast inertial relaxation engine)と呼ばれる方法を用いることができる。FIREに関しては、Erik Bitzek, Pekka Koskinen, Franz Gahler, Michael Moseler, and Peter Gumbsch, “Structural Relaxation Made Simple”, Physical Review Letters, 97, 170201 (2006)に説明されている。

0017

図1を参照しながら、実施例によるシミュレーション方法について説明する。ステップ10において、シミュレーション対象物を構成する各粒子の位置、速度、温度(粒子の位置における物体の温度)等、シミュレーションを行うに当たっての初期条件、及び拘束条件境界条件)を設定する。

0018

ステップ11において、各粒子の動解析を1タイムステップだけ実行する。例えば、各粒子の運動方程式を数値的に解くことにより、時間刻み幅Δt後における粒子の位置及び速度を求める。ステップ12において、動解析の結果に基づいて、発熱量を算出する。発熱量Pは、異なる部材の接触部において、以下の式を用いて算出される。
P=μ|F・vr|
ここで、μは摩擦係数、Fは接触部に加わる力、vrは相対速度である。

0019

ステップ13において、各粒子の発熱量を用いて、熱伝導方程式と等価な運動方程式を解くことにより、時間刻み幅Δt後における各粒子の温度を算出する。

0020

図2に、ステップ13のより詳細なフローチャートを示す。まず、ステップ131において、温度場の計算結果が定常状態に達したか否かを判定する。定常状態に達したか否かの判定は、例えば各粒子の温度の計算値の変化量と、判定閾値とを比較することにより行う。全ての粒子の温度の計算値の変化量が、判定閾値以下である場合、温度場が定常状態に達したと判定される。

0021

温度場が定常状態に達していない場合、ステップ132において、熱伝導方程式と等価な運動方程式の減衰項を規定するパラメータ(減衰係数)を、本来の値よりも小さくするか、0にする。温度場が定常状態に達するまでは、運動方程式の減衰係数を、本来の値よりも小さくするか、0にして運動方程式を解く。これにより、温度の変化を、本来の変化より急峻にすることができる。

0022

温度場が定常状態に達した場合、ステップ133において、熱伝導方程式と等価な運動方程式の減衰係数を、本来の値と等しくする。減衰係数の本来の値は、後に詳細に説明するように、シミュレーション対象物の密度と比熱との積に基づいて決定される。

0023

ステップ132またはステップ133の後、ステップ134において、熱伝導方程式と等価な運動方程式の数値計算を1タイムステップ実行する。このとき、運動方程式の減衰係数として、ステップ132またはステップ133で設定された値を用いる。ステップ132の後、図1のフローチャートに戻る。

0024

図1のステップ14において、シミュレーションの終了条件が満たされているか否かを判定する。シミュレーションの終了条件として、例えば、動解析の実行ステップ数等を採用することができる。シミュレーションの終了条件が満たされている場合には、シミュレーションを終了する。シミュレーションの終了条件が満たされていない場合には、ステップ15において、温度場が定常状態に達したか否かを判定する。

0025

ステップ15では、例えば各粒子の温度の速度が基準値以下になったら、温度場が定常状態に達したと判定される。温度場が定常状態に達したと判定された場合、ステップ11に戻って、動解析の次の1タイムステップを実行する。温度場が定常状態に達していないと判定された場合、ステップ16において、温度の速度を0に設定する条件が満たされているか否かを判定する。この判定は、熱伝導方程式と等価な運動方程式の数値計算のタイムステップごとに行われる。

0026

ステップ16において、判定条件が満たされていないと判定された場合、ステップ11に戻って、動解析を1タイムステップ実行する。ステップ16において、判定基準が満たされていると判定された場合、ステップ17において、各粒子の温度の速度を0に設定した後、ステップ11に戻って、動解析を1タイムステップ実行する。ステップ16及びステップ17の手法は、FIREと呼ばれる。

0027

次に、図3を参照して、FIREについて説明する。図3は、シミュレーション対象物のある注目点における温度の速度のシミュレーション結果の一例を示す。細い実線aは、通常の熱伝導方程式を数値的に解いて得られる温度の速度の一例を示す。図3には、時間の経過とともに、温度が徐々に低下する例が示されている。

0028

図3破線bは、ステップ13(図1)において熱伝導方程式と等価な運動方程式を数値的に解いて得られる温度の速度の一例を示す。後述するように、運動方程式の減衰係数は、定常状態における温度場にほとんど影響を与えない。従って、定常状態における温度分布を求めたい場合には、減衰係数を任意に設定することが可能である。

0029

温度場が定常状態になるまでの数値計算のステップ数を少なくするために、減衰係数を小さくすることが好ましい。減衰係数が小さいと、同一の外力が加わっても温度の計算値の変化が大きくなる。このとき、弾性項の影響が相対的に大きくなるため、温度の速度の計算値に振動が現れる。

0030

減衰係数を小さくして、熱伝導方程式と等価な運動方程式を解くと、破線bで示すように、数値演算によって求められた温度の計算結果の変化が急峻になる。さらに、温度の計算値は、振動しながら定常状態の温度に近づく。

0031

FIREを適用した場合の温度の計算値を太い実線cで示す。時刻t1の時点で、ステップ16(図1)の条件が満たされたと仮定する。時刻t1において、温度の速度が0に再設定(ステップ17)される。このため、時刻t1以降における温度の計算値の速度が、一旦緩やかになる。その結果、温度の計算値の振動が抑制され、定常状態に達するまでの時間が短くなる。

0032

ステップ16(図1)の条件の一例について説明する。各粒子において、温度の1次微分と2次微分と質量との積を求める。すべての粒子について求められた積の和を求める。ステップ16の条件として、この和が負であるこという条件が採用される。粒子ごとに求められた積が負であることは、各粒子の温度の速度の向きと加速度の向きとが異なること、すなわち温度変化が抑制されることを意味する。従って、この条件が満たされた場合、温度場が定常状態に近づきつつあると判断することができる。

0033

図4に、実施例によるシミュレーション方法を実行するシミュレーション装置(コンピュータ)のブロック図を示す。中央処理ユニット20、記憶装置21、入力装置22、及び出力装置23がデータバスライン24を介して相互に接続されている。

0034

入力装置22から、シミュレーションを実行するための前提条件となるパラメータが入力される。このパラメータには、例えば、粒子間ポテンシャル外場、密度、比熱、熱伝導係数等が含まれる。ステップ10(図1)で決定される粒子の初期状態、及び拘束条件も、入力装置22から入力される。入力装置22には、例えばキーボードを用いることができる。

0035

記憶装置21に、制御プログラム、シミュレーションプログラム、及びシミュレーションで数値計算される運動方程式を定義する種々のパラメータ等が格納されている。シミュレーションプログラムは、図1に示したシミュレーションを行う機能を実現する。中央処理ユニット20は、制御プログラム及びシミュレーションプログラムを実行する。シミュレーションが終了すると、中央処理ユニット20は、シミュレーション結果を出力装置23に出力する。例えば、三次元座標系内に各粒子の位置を表示し、各粒子の温度を色の濃淡で表すことができる。

0036

次に、ステップ13(図1)において、熱伝導方程式と等価な運動方程式を解くことにより温度場の解析を行うことができる理由について説明する。

0037

古典的スカラー場φのラグランジアンLは、例えばAdvanced Quantum Mechanics, (Addison Wesley, 1967)に説明されているように、下記の式で表される。



ここでσs[φ(x)]は拘束条件を表している。λsは、ラグランジェの未定定数である。δは、ディラックのδ関数である。δに付された指数3は、及びdxに付された指数3は、三次元を意味する。Σは、xに関して足し合わせることを意味する。xがxs以外のときは、δ関数の値は0である。式(1)において、ρ(x)、κ、Q(x)は、特定の物理量を意味していない。

0038

拘束条件σs[φ(x)]として、ディレクリ境界条件、ノイマン境界条件を採用することができる。ディレクリ境界条件は、下記の式で表すことができる。




式(2)は、x=xsにおいてφ(x)が定数であることを意味する。

0039

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




式(3)は、x=xsにおいてφ(x)の傾きが定数であることを意味する。φ(xs)をφsと表している。∇sは、x=xsの位置における傾きを表す。φs0は、拘束条件に基づく定数である。

0040

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




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




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

0041

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




ここで、γiは減衰係数である。式(4)、式(6)から、i番目の粒子に対する以下の運動方程式が得られる。




ここで、式(7)の右辺第1項は弾性項を記述し、第2項は減衰項を記述し、第3項は外力を記述し、第4項は拘束条件を記述している。右辺の第1項のΣのj=±1は、i番目の粒子の両側の粒子について和を取ることを意味する。右辺の第4項のΣのsは、拘束条件の数に相当する。

0042

(d2/dt2)φiが0になる定常状態において、式(7)は、以下のように変形できる。

0043

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

0044

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




ここで、ρは密度、Cvは比熱、Tは温度、tは時間、Kは熱伝導率、Qは単位体積当たりの発熱量を表す。

0045

式(8)のφi、γi、κi、及びQiを、それぞれ温度T、密度ρと比熱Cvとの積、熱伝導率K、及び発熱量Qとみなせば、式(8)が式(9)の熱伝導方程式と同形であることがわかる。具体的には、式(8)の空間的なφiの分布の項、及びφiの時間微分の項を、それぞれ式(9)の空間的な温度分布の項、及び温度の時間微分の項に対応付けることができる。このため、ρiが十分小さいという条件の下で、式(7)の運動方程式が式(9)の熱伝導方程式と等価であるといえる。

0046

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




ここで、式(10)の右辺第1項のΣは、i番目の粒子に接触している全ての粒子について足し合わせることを意味する。ΔSijはボロノイ多面体iとjとの接触面の面積を表す。miは下記の式で表される。




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

0047

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

0048

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

0049

式(10)のγiを、密度ρiと比熱Cviとの積に置き換え、φiを温度Tiに置き換えることにより、以下の式が得られる。Cviは、i番目の粒子の位置の比熱を表す。

0050

運動方程式と同形の式(14)を解くことにより、温度場のシミュレーションを行うことができる。式(10)のφiを式(14)においてTiに置き換えたため、式(14)のmiの次元は、式(10)、(11)のmiの次元とは異なる。仮想質量miの次元は、式(14)の右辺のΔViを含む項から、密度ρi、比熱Cv、体積素ΔVi、及び時間刻み幅dtの積の次元と同一であることがわかる。すなわち、仮想質量miの次元は、[J/K][s]となる。

0051

式(12)から、仮想質量miの満たすべき条件は、以下の式で与えられる。

0052

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

0053

式(15)及び式(16)から、一例として、miを下記の式のように設定することが推奨される。

0054

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

0055

ステップ13(図1)において、運動方程式(7)の減衰係数γi、式(14)においてはρiCviを0にすれば、温度Tの速度の算出結果が大きくなる。すなわち、温度Tの時間変化が急峻になる。このため、温度の計算値が定常状態における温度に達するまでの時間を短くすることができる。

0056

次に、ステップ132(図2)において、減衰係数γiまたはρiCviを、本来の値とは異なる値、または0に設定しても、定常状態における正確な温度を求めることが可能なことを示す。

0057

図5に示した棒状のモデルを考える。体積dViの中央領域30に、単位体積あたりの発熱量Qiの熱源が存在する。中央領域30から、面積Siの界面を通って、両側の中間領域31、32に向って熱伝導が生じる。中間領域31、32の端に、それぞれ端部領域33、34が連続している。端部領域33、34が、それぞれ面積Sjの端面において温度T0の熱浴に接している。

0058

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




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

0059

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

0060

式(18)〜式(20)から、以下の式が導出される。

0061

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

0062

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

0063

減衰係数ρiCviを0にすると振動が減衰せず、いつまでも振動が継続してしまう。実施例においては、FIREを適用することにより、図3の実線cで示したように、温度の計算値の振動を短時間で収束させている。

0064

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

0065

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

0066

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

0067

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

0068

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

0069

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

0070

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

0071

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

0072

図7に示した例では、例えば温度の計算値が250K以上の領域では、ステップ16(図1)において条件が満たされていないと判定されている。このため、温度の計算値が急峻に低下している。温度の計算値が250K以下になった最初のタイムステップで、ステップ16(図1)において条件が満たされていると判定されている。このため、温度の計算値の低下測度が急激に小さくなっている。

0073

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

0074

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

0075

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

0076

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

0077

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

0078

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

0079

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

0080

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




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

0081

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

0082

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

0083

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

0084

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

0085

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

0086

具体的には、熱伝導方程式と等価な運動方程式(式(14))の減衰係数(Tiドット係数)であるρiCviΔViを現時点の値より小さくする。一例として、減衰係数を、現時点の値の0.5倍とする。減衰係数を小さくすることにより、温度変化をより急峻にすることができる。

0087

ステップ1351において加速アルゴリズム適用条件が満たされていないと判定された場合、すなわち、温度Tiの速度の符号と温度Tiの加速度の符号とが異なる場合、ステップ1353において、熱伝導方程式と等価な運動方程式(式(14))の減衰係数(Tiドットの係数)であるρiCviΔViを、本来の値に戻す。

0088

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

0089

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

0090

図9に、シミュレーション対象である簡易的なベアリングの断面図を示す。シミュレーション対象のベアリングは、内輪50、外輪51、及び両者の間に配置された12個のころ52で構成される。外輪51の外周面は、温度300Kの外気に触れており、内輪50を強制的に3600rpmの回転速度で回転させた。外輪51の外周面と外気との界面における熱伝達率は5,000W/m2/Kとした。内輪50、外輪51、及びころ52の熱伝導率として、一般的なステンレス鋼の値を用いた。

0091

粒子iと粒子jとの接触力をfij、粒子iと粒子jとの相対速度をvij、粒子iと粒子jとの間の摩擦係数をμijで表したとき、粒子iの単位体積当たりの発熱量Qiは、以下の式で表される。

0092

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

0093

図10Aと図10Bとを比較すると、図10Aのシミュレーション結果の方が、定常状態に到達するまでの時間が著しく短いことが分かる。これは、図8のステップ1352において、減衰係数を小さくしたことの効果である。図10Bのシミュレーション結果では、全体温度が319Kに収束しているが、図10Aのシミュレーション結果では、定常状態に到達した後も、幅1.0K程度の範囲内で全体温度が振動していることがわかる。

0094

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

0095

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

0096

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

0097

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

0098

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

0099

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

0100

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

0101

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

0102

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

0103

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

0104

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

0105

図14に、図13に示した実施例によるシミュレーション方法を用いて図9に示したベアリングの温度変化を求めた結果を示す。横軸は経過時間を単位「s」で表し、縦軸は温度を単位「K」で表す。図14の太い実線dが、図13に示した実施例によるシミュレーション方法を用いた場合の全体温度の計算結果を示し、細い実線eが、図8に示した参考例によるシミュレーション方法を用いた場合の全体温度の計算結果を示す。

0106

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

0107

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

0108

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

0109

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

0110

10〜17 ステップ
20中央処理ユニット
21記憶装置
22入力装置
23出力装置
24データバスライン
30中央領域
31、32 中間領域
33、34 端部領域
41、42シミュレーション対象の物体
50内輪
51外輪
52コロ
131〜134、136 ステップ
1351〜1355 ステップ

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

(分野番号表示ON)※整理標準化データをもとに当社作成

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ