JP7244409B2

JP7244409B2 - シミュレーション装置及びプログラム

Info

Publication number: JP7244409B2
Application number: JP2019227032A
Authority: JP
Inventors: 公則坂井
Original assignee: Sumitomo Heavy Industries Ltd
Current assignee: Sumitomo Heavy Industries Ltd
Priority date: 2019-12-17
Filing date: 2019-12-17
Publication date: 2023-03-22
Anticipated expiration: 2039-12-17
Also published as: JP2021096615A; US20210182457A1

Description

本発明は、シミュレーション装置及びプログラムに関する。

接触する２つの部材の構造解析において、クーロン摩擦力の大きさは、接触面に対して垂直方向に作用する垂直力と摩擦係数との積で表し、クーロン摩擦力の向きは、接触面の相対速度の方向に応じて決定される。

２つの部材の接触面の各部にすべりが生じているか否かのすべり状態を解析する状態解析方法が下記の特許文献１に開示されている。

特開２０１６－９０３２４号公報

動的陽解法を用いた粒子法、繰込み群分子動力学法では力の釣り合いを解かないため、接触面に位置する粒子（節点）は振動を繰り返しながら摺動する。接触する２つの面は、平均的には摺動方向に移動しているが、個々の粒子には速度の揺らぎが発生している。接触面に位置する個々の粒子について、揺らぎを持つ相対速度に基づいて摩擦力を決定すると、タイムステップごとの摩擦力の向きが様々な方向を向き、クーロン摩擦力を適切に表現できない。例えば、本来、クーロン摩擦力はマクロ的な摺動方向に作用するが、速度の揺らぎを考慮して求めた摩擦力は、摺動方向に対して直交する方向の成分も持つことになる。その結果、摺動方向に作用する摩擦力は、本来発生する摩擦力より小さくなってしまう。

本発明の目的は、適切なクーロン摩擦力を導入して、動的陽解法により２つの部材の構造解析を行うことが可能なシミュレーション装置及びプログラムを提供することである。

本発明の一観点によると、
シミュレーション条件が入力される入力装置と、
前記入力装置に入力されたシミュレーション条件に基づいて、部材を複数の粒子の集まりで表して粒子法を用いて構造解析を行う処理装置と、
出力装置と
を有し、
前記処理装置は、
解析対象の２つの部材が接触している面に位置する複数の粒子にそれぞれ作用する摩擦力の向きを、一方の部材の複数の粒子で定義される基準点と、他方の部材の複数の粒子で定義される基準点との、１タイムステップごとの相対変位ベクトルを積算して得られる積算変位ベクトルに基づいて決定し、
決定された摩擦力に基づいて、前記複数の粒子について運動方程式を解き、
解析結果を前記出力装置に出力するシミュレーション装置が提供される。

本発明の他の観点によると、
解析対象の２つの部材を複数の粒子の集まりで表した解析モデルにおいて、解析対象の２つの部材が接触している面に位置する複数の粒子にそれぞれ作用する摩擦力の向きを、一方の部材の複数の粒子で定義される基準点と、他方の部材の複数の粒子で定義される基準点との、１タイムステップごとの相対変位ベクトルを積算して得られる積算変位ベクトルに基づいて決定する機能と、
決定された摩擦力に基づいて、前記複数の粒子について運動方程式を解く機能と
をコンピュータに実現させるプログラムが提供される。

積算変位ベクトルに基づいて摩擦力の向きを決定することにより、実際の摩擦力を適切に再現したクーロン摩擦力を用いてシミュレーションを行うことができる。

図１Ａは、解析対象の２つの部材の一例を示す模式図であり、図１Ｂは、第１部材及び第２部材を、それぞれ複数の粒子の集まりで表した解析モデルの一部を示す図である。図２は、実施例によるシミュレーション装置のブロック図である。図３は、実施例によるシミュレーション方法のフローチャートである。図４Ａは、接触面に位置する複数の粒子によって定義される三角形要素を示す模式図であり、図４Ｂは、第１部材の複数の三角形要素を示す模式図である。図５Ａは、第１部材の粒子（図１Ｂ）に対する第２部材の粒子の相対的な位置の変化を模式的に示す図であり、図５Ｂは、静止摩擦状態のときに接触面に位置する粒子に作用する摩擦力の物理的意味を説明するための模式図である。図６は、シミュレーションモデルを示す模式図である。図７Ａは、第１部材と第２部材との接触面に発生する摩擦力の向きを、積算変位ベクトルと平行にしたときの摩擦力の大きさの時間変化を示すグラフであり、図７Ｂは、摩擦力の向きを、相対速度ベクトルの向きと平行にしたときの摩擦力の大きさの時間変化を示すグラフである。図８Ａは、積算変位ベクトルの大きさを判定上限値に制限することの効果を確認するために行ったシミュレーションにおける滑り部材の移動の軌跡を示すグラフであり、図８Ｂは、回転角度θ_ｖ、θ_Ｆの時間変化を示すグラフである。

図面を参照しながら、本発明の実施例によるシミュレーション方法、シミュレーション装置、及びプログラムについて説明する。

図１Ａは、解析対象の２つの部材の一例を示す模式図である。第１部材１１と第２部材１２とが、接触面１５で相互に接触している。接触面１５に対して垂直な方向に垂直力Ｆ_Ｎが印加されている。第１部材１１に対して第２部材１２を接触面１５に対して平行な方向に速度ｖで移動させる。このとき、接触面１５に摩擦力Ｆ_Ｓが発生する。第１部材１１に作用する摩擦力Ｆ_Ｓは速度ｖと同一の向きであり、第２部材１２に作用する摩擦力Ｆ_Ｓは速度ｖとは反対向きである。

図１Ｂは、第１部材１１及び第２部材１２を、それぞれ複数の粒子２１、２２の集まりで表した解析モデルの一部を示す図である。第１部材１１の複数の粒子２１がバネ２４で相互に接続されており、第２部材１２の複数の粒子２２がバネ２５で相互に接続されている。接触面１５に位置する粒子２１、２２に摩擦力Ｆが作用する。接触面１５に位置する粒子２１、２２の各々に作用する摩擦力Ｆを合計した値は、それぞれ第１部材１１、第２部材１２に作用する摩擦力Ｆ_Ｓに等しい。

第１部材１１の底面に位置する複数の粒子２１の位置を固定し、第２部材１２の最上面に位置する複数の粒子２２に、垂直力Ｆ_Ｎを分散して印加することにより、第２部材１２に印加する垂直力Ｆ_Ｎが再現される。第２部材１２の１つの側面に位置する複数の粒子２２を、強制的に速度ｖで移動させることにより、第１部材１１に対する第２部材１２の相対的な速度ｖが再現される。摩擦力Ｆについては、後に図４Ａ～図５Ｂを参照して説明する。

図２は、実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力装置５０、処理装置５１、出力装置５２、及び外部記憶装置５３を含む。入力装置５０から処理装置５１にシミュレーション条件等が入力される。さらに、オペレータから入力装置５０に各種指令（コマンド）等が入力される。入力装置５０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。

処理装置５１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法またはくりこみ群分子動力学法を用いたシミュレーションを行う。処理装置５１は、中央処理ユニット（ＣＰＵ）、主記憶装置（メインメモリ）等を含むコンピュータで実現される。コンピュータが実行するシミュレーションプログラムが、外部記憶装置５３に記憶されている。外部記憶装置５３には、例えばハードディスクドライブ（ＨＤＤ）、ソリッドステートドライブ（ＳＳＤ）等が用いられる。処理装置５１は、外部記憶装置５３に記憶されているプログラムを主記憶装置に読み出して実行する。

処理装置５１は、シミュレーション結果を出力装置５２に出力する。シミュレーション結果には、解析対象の部材を表す複数の粒子の状態、複数の粒子からなる粒子系の物理量の時間的変化等を表す情報が含まれる。出力装置５２は、例えば通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

図３は、実施例によるシミュレーション方法のフローチャートである。
まず、処理装置５１が、入力装置５０に入力されたシミュレーション条件を取得する（ステップＳ１）。シミュレーション条件には、第１部材１１及び第２部材１２の幾何学的形状や相対位置関係を定義する情報、第１部材１１及び第２部材１２の物性情報、摩擦係数、第１部材１１及び第２部材１２に加わる外力、速度等が含まれる。

処理装置５１は、シミュレーション条件を取得すると、シミュレーション対象の第１部材１１及び第２部材１２を複数の粒子の集まりで表した解析モデルを定義する（ステップＳ２）。その後、解析モデルの接触面１５（図１Ａ、図１Ｂ）に位置する複数の粒子２１、２２にそれぞれ作用する摩擦力Ｆを算出する（ステップＳ３）。摩擦力Ｆの求め方については、後に図４Ａ～図５Ｂを参照して説明する。

次に、各粒子２１、２２について運動方程式を解く（ステップＳ４）。このとき、接触面１５に位置する粒子については、摩擦力Ｆを作用させる。これにより、１タイムステップが経過した時点の各粒子２１、２２の状態が得られる。解析を継続する場合には、摩擦力Ｆを再度算出（ステップＳ３）し、運動方程式を解いて（ステップＳ４）タイムステップを進める。解析を終了する場合は、解析結果を出力装置５２に出力する（ステップＳ５）。

次に、図４Ａ～図５Ｂを参照して、接触面１５（図１Ｂ）に位置する粒子２１、２２に作用する摩擦力Ｆの算出方法について説明する。

まず、図４Ａを参照して、摩擦力を及ぼし合う第１部材１１の複数の粒子２１と、第２部材１２の複数の粒子２２との組み合わせについて説明する。

図４Ａは、接触面１５に位置する複数の粒子２１、２２によって定義される三角形要素を示す模式図である。三角形要素の生成は、周知の種々のアルゴリズムを用いて行うことができる。第１部材１１の３個の粒子２１によって１つの三角形要素３１が定義される。第２部材１２の複数の粒子によって定義される複数の三角形要素３２（図４Ａでは５個の三角形要素３２）が、平面視において、１つの三角形要素３１と部分的に重なっている。部分的に重なり合っている三角形要素３１と三角形要素３２との重心の相対変位によって三角形要素の間に摩擦力が発生すると考える。

第１部材１１の三角形要素３１は、この三角形要素３１に部分的に重なる第２部材１２の複数の三角形要素３２の各々から摩擦力を受ける。この摩擦力を合成することにより、三角形要素３１に作用する摩擦力Ｆｊを求めることができる。

図４Ｂは、第１部材１１の複数の三角形要素３１を示す模式図である。第１部材１１の１つの粒子２１を１つの頂点とする複数の三角形要素３１（図４Ｂでは６個の三角形要素３１）が定義される。着目する粒子２１Ａを含む６個の三角形要素３１の各々について、三角形要素３１に作用する摩擦力Ｆｊ（ｊ＝１、２、３、４、５、６）を、その頂点に位置する３個の粒子２１に、粒子２１の質量で重みづけして分配する。３個の粒子２１の質量が等しい場合には、粒子２１の各々に分配される摩擦力は三角形要素３１に作用する摩擦力Ｆｊの１／３になる。着目する粒子２１Ａに作用する摩擦力Ｆは、その粒子２１Ａに分配された摩擦力を合成することにより得られる。

次に、図５Ａ及び図５Ｂを参照して、三角形要素に作用する摩擦力の求め方について説明する。

図５Ａは、第２部材１２の１つの三角形要素３２（図４Ａ）の重心に対する第１部材１１の１つの三角形要素３１の重心３３の相対的な位置の変化を模式的に示す図である。三角形要素３１の３つの頂点に位置する粒子２１のそれぞれの速度ベクトルを、それぞれｖ１、ｖ２、ｖ３と表記し、それぞれの質量をｍ１、ｍ２、ｍ３と表記すると、重心３３の速度ベクトルｖｃは以下の式で表される。

第１部材１１の三角形要素３１の重心３３の速度ベクトルｖｃと、第２部材１２の三角形要素３２の重心の速度ベクトルとから、第２部材１２の１つの三角形要素３２（図４Ａ）の重心に対する第１部材１１の１つの三角形要素３１の重心３３の相対速度ベクトルを求めることができる。

重心３３の初期状態からの、接触面１５に平行な方向の１タイムステップごとの相対変位ベクトルを積算した積算変位ベクトルをｕと表記し、接触面１５に平行な方向の相対速度ベクトルをｖと表記する。初期状態における積算変位ベクトル及び相対速度ベクトルを、それぞれｕ（０）、ｖ（０）と表記し、ｉ番目のタイムステップを実行した後の積算変位ベクトル及び相対速度ベクトルを、それぞれｕ（ｉ）、ｖ（ｉ）と表記する。タイムステップの時間幅をｄｔと表記する。ｉ番目のタイムステップにおける相対変位ベクトルは、ｖ（ｉ－１）ｄｔと表すことができる。

ｉ番目のタイムステップ実行後の積算変位ベクトルｕ（ｉ）は、以下の式で表される。

積算変位ベクトルｕ（ｉ）の大きさが、あらかじめ決められている判定上限値ｕ_ｍａｘを超えたら、積算変位ベクトルｕ（ｉ）の大きさが判定上限値ｕ_ｍａｘに等しくなるように、積算変位ベクトルｕ（ｉ）を補正する。すなわち、下記の補正を行う。

積算変位ベクトルｕ（ｉ＋１）は、補正後の積算変位ベクトルｕ（ｉ）に基づいて算出される。図５Ａに示した例では、積算変位ベクトルｕ（２）の大きさが判定上限値ｕ_ｍａｘを超えているため、積算変位ベクトルｕ（２）の大きさをｕ_ｍａｘに補正している。積算変位ベクトルｕ（３）は、補正後の積算変位ベクトルｕ（２）に基づいて算出される。

次に、第１部材１１に対して第２部材１２が移動している動摩擦状態（滑り状態）のときに重心３３に作用する摩擦力Ｆｊについて説明する。ｉ番目のタイムステップ実行後の状態で、重心３３に作用する動摩擦状態の摩擦力Ｆｊ（ｉ）を、以下の式を用いて算出する。

ここで、μは動摩擦係数であり、Ｆ_Ｎ（ｉ）は三角形要素３１に作用する垂直力である。三角形要素３１に作用する垂直力Ｆ_Ｎ（ｉ）は、三角形要素３１と３２との重なり部分の面積に応じて算出される。なお、第２部材１２の三角形要素３２（図４Ａ）の重心に作用する摩擦力Ｆｊの向きは、第１部材１１の三角形要素３１の重心３３に作用する摩擦力Ｆｊの向きと反対向きである。式（４）は、摩擦力Ｆｊ（ｉ）の大きさが、動摩擦係数μと垂直力Ｆ_Ｎの大きさとの積に等しいことを意味する。さらに、摩擦力Ｆｊ（ｉ）の向きは、積算変位ベクトルｕ（ｉ）の向きと逆向きであることを意味する。

次に、第２部材１２に対して第１部材１１が固着している静止摩擦状態（固着状態）のときの摩擦力Ｆｊについて説明する。ｉ番目のタイムステップ実行後の状態で、第１部材１１の三角形要素３１の重心３３に作用する静止摩擦状態の摩擦力Ｆｊ（ｉ）として、以下の式を用いて算出する。

ここでＫは比例定数である。

次に、式（５）の物理的な意味について説明する。
図５Ｂは、静止摩擦状態のときに三角形要素３１、３２の重心３３、３４に作用する摩擦力Ｆｊの物理的意味を説明するための模式図である。静止摩擦状態のときには、重心３３と重心３４とが、バネ２６で接続されていると考える。式（５）の比例定数Ｋは、このバネ２６のバネ定数に相当する。すなわち、重心３４に対して重心３３が接触面１５に平行な方向に変位すると、バネ２６によって重心３３を元の位置に戻そうとする復元力が作用する。

次に、式（３）に含まれる判定上限値ｕ_ｍａｘの大きさの一例について説明する。
第１部材１１と第２部材１２とが相互に固着された静止摩擦状態のときに、第１部材１１に対して第２部材１２に滑り方向の力を加えると、両者の接触面１５に剪断力が発生する。剪断力が最大静止摩擦力以下のときには、固着状態が維持される。剪断力が最大静止摩擦力を越えると、第２部材１２が第１部材１１に対して滑り始める。

判定上限値ｕ_ｍａｘを、剪断力と最大静止摩擦力とが釣り合った状態の相対変位ベクトルの大きさと定義すると、判定上限値ｕ_ｍａｘは以下の式で表される。

ここで、Ｋは、式（５）の比例定数Ｋを表しており、Ｆ_Ｎ（ｉ）は、ｉ番目のタイムステップ後における接触面１５に対して垂直な方向に作用する垂直力である。一例として、式（６）を満たすように、判定上限値ｕ_ｍａｘを設定すればよい。

次に、上記実施例の優れた効果について説明する。
２つの部材の接触面に発生する摩擦力の方向は、２つの部材の相対速度の方向と平行及び反平行である。粒子法を用いて構造解析を行う際に、接触面に位置する粒子の相対速度に基づいて、各粒子に摩擦力を作用させると、摩擦力の向きが相対速度の揺らぎの影響を受けてしまう。この揺らぎの影響を受けると、各粒子に作用する摩擦力の向きが、２つの部材の接触面に発生する実際の摩擦力の向きからずれてしまう。

上記実施例では、図５Ａ、及び式（４）、式（５）に示したように、摩擦力Ｆｊ（ｉ）の向きを、相対速度ベクトルｖ（ｉ）ではなく、積算変位ベクトルｕ（ｉ）に基づいて決定している。積算変位ベクトルｕ（ｉ）には、過去の複数の相対速度ベクトルｖ（ｉ）が累積されている。すなわち、相対速度ベクトルｖ（ｉ）の揺らぎも累積されている。相対速度ベクトルｖ（ｉ）の揺らぎはランダムであるため、揺らぎを累積すると複数の揺らぎが打ち消し合って、揺らぎの影響が軽減される。このため、積算変位ベクトルｕ（ｉ）は、一時的な相対速度ベクトルｖ（ｉ）の揺らぎの影響を受けにくい。図４Ｂに示したように、第１部材１１の三角形要素３１の重心３３に作用する摩擦力Ｆｊ（ｉ）を合成することによって、粒子２１に作用する摩擦力Ｆを求めるため、粒子２１に作用する摩擦力Ｆは、粒子２１の速度ベクトルの揺らぎの影響を受けにくくなる。第２部材１２の粒子２２に作用する摩擦力Ｆについても同様である。したがって、各粒子２１、２２に作用する摩擦力Ｆ（ｉ）の向きと、第１部材１１と第２部材１２との接触面１５に発生する摩擦力の向きとの相違を小さくすることができる。これにより、構造解析の精度を高めることができる。

計算のタイムステップが進み、１タイムステップ分の相対変位ベクトルｖ（ｉ－１）ｄｔの大きさに比べて、積算変位ベクトルｕ（ｉ）が大きくなり過ぎると、現時点の相対速度の向きが摩擦力Ｆｊの向きに反映されにくくなる。上記実施例では、式（３）に示したように、積算変位ベクトルｕ（ｉ）の大きさが判定上限値ｕ_ｍａｘを超えたら、積算変位ベクトルｕ（ｉ）を、その大きさが判定上限値ｕ_ｍａｘに等しくなるように補正している。このため、現時点の相対速度の向きが摩擦力Ｆの向きに反映されにくくなるという不都合を回避することができる。

判定上限値ｕ_ｍａｘを大きくし過ぎると、積算変位ベクトルｕ（ｉ）の大きさを判定上限値ｕ_ｍａｘに制限することの効果が薄れてしまう。判定上限値ｕ_ｍａｘは、２つの部材の接触面１５に位置する粒子２１、２２に作用する剪断力と、最大静止摩擦力とが釣り合った時の相対変位ベクトルの大きさ以下とすることが好ましい。

逆に、判定上限値ｕ_ｍａｘを小さくし過ぎると、摩擦力Ｆｊの向きは相対速度ベクトルｖの揺らぎの影響を受けやすくなる。判定上限値ｕ_ｍａｘは、１タイムステップの相対変位ベクトルｖ（ｉ）ｄｔの揺らぎの大きさに比べて十分大きくすることが好ましい。

次に、図６～図７Ｂを参照して、上記実施例の効果を確認するために行った実際のシミュレーション及びその結果について説明する。

図６は、シミュレーションモデルを示す模式図である。板状の第１部材１１の上に板状の第２部材１２が載せられており、第１部材１１の上面と第２部材１２の下面とが接触している。第２部材１２の上面に垂直力Ｆ_Ｎを与える。第１部材１１の下面の高さ方向の位置は固定されている。第１部材１１と第２部材１２との接触面に摩擦力Ｆ_Ｓが発生する。

第１部材１１の１つの側面を、速度ｖで強制的に移動させる。第２部材１２の、速度ｖとは反対方向を向く側面に、バネ１７が取り付けられており、バネ１７が伸びると、速度ｖとは反対向きの復元力が第２部材１２の側面に作用する。

図７Ａは、第１部材１１と第２部材１２との接触面に発生する摩擦力の向きを、積算変位ベクトルｕ（図５Ａ）と平行にしたときの摩擦力の大きさの時間変化を示すグラフである。横軸は経過時間を単位「ｓ」で表し、縦軸は摩擦力の大きさを単位「Ｎ」で表す。なお、垂直力Ｆ_Ｎの大きさを１０００Ｎ、静止摩擦係数及び動摩擦係数を共に０．３とした。第１部材１１及び第２部材１２の物性値として、鉄の物性値を用いた。

第１部材１１の移動開始直後は、第２部材１２が第１部材１１に追随して移動する。第２部材１２の移動量が大きくなるにしたがって、バネ１７が伸ばされてバネ１７の復元力が大きくなる。バネ１７の復元力が最大静止摩擦力を越えると、第１部材１１に対して第２部材１２が滑り始める。このときの摩擦力は、垂直力Ｆ_Ｎの大きさ１０００Ｎに、静止摩擦係数０．３を乗じた３００Ｎである。

図７Ａに示したグラフにおいて、摩擦力の大きさが時間の経過とともに負の向きに大きくなっている期間は、第１部材１１に対して第２部材１２が固着している状態に対応する。第１部材１１の強制移動開始時点から約０．１３秒経過した時点で、摩擦力の大きさが最大静止摩擦力である３００Ｎに達する。この時点が、第２部材１２が第１部材１１に対して滑り始めた時点に対応する。本シミュレーションでは、動摩擦係数と静止摩擦係数とを同一に設定しているため、第２部材１２が滑り始めた後は、摩擦力の大きさが約３００Ｎで一定になっている。この結果から、実施例によるシミュレーション方法では、各粒子に作用させる摩擦力が、実際の摩擦力をほぼ正確に再現できていることがわかる。

図７Ｂは、摩擦力の向きを、相対速度ベクトルｖ（図５Ａ）の向きと平行にしたときの摩擦力の大きさの時間変化を示すグラフである。図７Ｂに示した例では、時間の経過とともに摩擦力の大きさが大きくなっているが、摩擦力の大きさは、最大静止摩擦力に達することなく低下し始める。これは、シミュレーションの各タイムステップにおいて算出された摩擦力が、本来の摩擦力を再現していないためである。

図６～図７Ｂに示したシミュレーション結果から、上記実施例によるシミュレーション方法において、クーロン摩擦力を十分正確に再現できていることが確認された。

次に、式（２）で示したように、積算変位ベクトルｕ（ｉ）の大きさを判定上限値ｕ_ｍａｘに制限することの効果を確認するために行ったシミュレーション及びその結果について、図８Ａ及び図Ｂを参照して説明する。本シミュレーションでは、支持面の上に滑り部材を載せて、支持面に対して滑り部材を強制的に移動させた。

図８Ａは、本シミュレーションにおいて行った滑り部材４０の移動の軌跡を示すグラフである。滑り部材４０は、ｘ軸方向に等速度で移動させ、ｙ軸方向に単振動させた。これにより、滑り部材４０の軌跡は正弦波の形状になる。ｘ軸から、滑り部材４０の相対速度ベクトルｖまでの回転角度をθ_ｖと表記する。ｘ軸から、支持面に作用する摩擦力Ｆまでの回転角度をθ_Ｆと表記する。

図８Ｂは、回転角度θ_ｖ、θ_Ｆの時間変化を示すグラフである。図８Ｂに示した実線が回転角度θ_ｖ、θ_Ｆを示す。両者はほぼ一致している。比較のために、積算変位ベクトルｕ（ｉ）の大きさに制限を設けないでシミュレーションを行った。比較例によるシミュレーションにおいて算出されたｘ軸から摩擦力Ｆまでの回転角度をθ_ＮＦと表記する。図８Ｂにおいて、回転角度θ_ＮＦを破線で示している。比較例においては、摩擦力Ｆの向きが相対速度ベクトルｖの向きから大きくずれていることがわかる。

図８Ａ及び図８Ｂに示したシミュレーションにより、積算変位ベクトルｕ（ｉ）の大きさが判定上限値ｕ_ｍａｘを超えないように、積算変位ベクトルｕ（ｉ）を式（２）に基づいて補正することの効果が確認された。

上記実施例では、積算変位ベクトルｕ（ｉ）が判定上限値ｕ_ｍａｘを超えたら、積算変位ベクトルｕ（ｉ）の大きさを判定上限値ｕ_ｍａｘに等しくしているが、求められた積算変位ベクトルｕ（ｉ）の大きさより小さくなるように補正してもよい。

上記実施例では、２つの部材の接触面に位置する粒子を頂点とする三角形要素を定義し、三角形要素の対ごとに摩擦力を求めたが、三角形要素以外の多角形要素を定義し、多角形要素の対ごとに摩擦力を求めてもよい。また、上記実施例では、三角形要素の対ごとに相対変位ベクトルを求める基準点として三角形要素の重心を採用したが、その他の点を基準点として採用してもよい。例えば、三角形要素の幾何学的中心位置等を基準点として採用してもよい。このように、複数の粒子の位置に基づいて基準点を定義すればよい。

上記実施例では、２つの部材の接触面が平面である解析モデルについてシミュレーションを行ったが、接触面は、必ずしも平面である必要はない。例えば、上記実施例によるシミュレーション方法は、滑り軸受の接触面や、ピストンとシリンダとの接触面のような柱面である部材の構造解析にも適用することが可能である。さらに、接触面が球面である部材の構造解析にも適用することが可能である。

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

１１第１部材
１２第２部材
１５接触面
１７バネ
２１第１部材を表す粒子
２１Ａ着目する粒子
２２第２部材を表す粒子
２４第１部材を表す２つの粒子を接続するバネ
２５第２部材を表す２つの粒子を接続するバネ
２６第１部材の接触面上の粒子と第２部材の接触面上の粒子とを接続するバネ
３１、３２三角形要素
３３、３４三角形要素の重心
４０滑り部材
５０入力装置
５１処理装置
５２出力装置
５３外部記憶装置

Claims

シミュレーション条件が入力される入力装置と、
前記入力装置に入力されたシミュレーション条件に基づいて、部材を複数の粒子の集まりで表して粒子法を用いて構造解析を行う処理装置と、
出力装置と
を有し、
前記処理装置は、
解析対象の２つの部材が接触している面に位置する複数の粒子にそれぞれ作用する摩擦力の向きを、一方の部材の複数の粒子で定義される基準点と、他方の部材の複数の粒子で定義される基準点との、１タイムステップごとの相対変位ベクトルを積算して得られる積算変位ベクトルに基づいて決定し、
決定された摩擦力に基づいて、前記複数の粒子について運動方程式を解き、
解析結果を前記出力装置に出力するシミュレーション装置。
前記基準点は、解析対象の２つの部材が接触している面に位置する複数の粒子の位置で定義される複数の三角形要素の各々の重心である請求項１に記載のシミュレーション装置。
前記処理装置は、
解析対象の一方の部材の複数の粒子で定義される前記三角形要素と、他方の部材の複数の粒子で定義される前記三角形要素との対ごとに、摩擦力を決定し、
前記三角形要素の対ごとに決定された摩擦力を、前記三角形要素の頂点に位置する粒子に分配して、粒子に作用する摩擦力を決定する請求項２に記載のシミュレーション装置。
前記処理装置は、
あるタイムステップで前記積算変位ベクトルの大きさが判定上限値を超えたら、前記積算変位ベクトルの大きさを、求められた前記積算変位ベクトルの大きさより小さく補正し、次のタイムステップにおいて、補正後の前記積算変位ベクトルに、当該タイムステップにおける相対変位ベクトルを加える請求項１乃至３のいずれか１項に記載のシミュレーション装置。
前記処理装置は、
あるタイムステップで前記積算変位ベクトルの大きさが前記判定上限値を超えたら、前記積算変位ベクトルの大きさを、前記判定上限値に一致するように補正する請求項４に記載のシミュレーション装置。
前記処理装置は、
解析対象の２つの部材が静止摩擦状態のときに、２つの部材の接触面に位置する粒子に作用する剪断力と、最大静止摩擦力とが釣り合った時の相対変位ベクトルの大きさを、前記判定上限値とする請求項５に記載のシミュレーション装置。
解析対象の２つの部材を複数の粒子の集まりで表した解析モデルにおいて、解析対象の２つの部材が接触している面に位置する複数の粒子にそれぞれ作用する摩擦力の向きを、一方の部材の複数の粒子で定義される基準点と、他方の部材の複数の粒子で定義される基
準点との、１タイムステップごとの相対変位ベクトルを積算して得られる積算変位ベクトルに基づいて決定する機能と、
決定された摩擦力に基づいて、前記複数の粒子について運動方程式を解く機能と
をコンピュータに実現させるプログラム。