JP5893333B2 - Particle behavior analysis method, particle behavior analysis apparatus, and analysis program - Google Patents
Particle behavior analysis method, particle behavior analysis apparatus, and analysis program Download PDFInfo
- Publication number
- JP5893333B2 JP5893333B2 JP2011230304A JP2011230304A JP5893333B2 JP 5893333 B2 JP5893333 B2 JP 5893333B2 JP 2011230304 A JP2011230304 A JP 2011230304A JP 2011230304 A JP2011230304 A JP 2011230304A JP 5893333 B2 JP5893333 B2 JP 5893333B2
- Authority
- JP
- Japan
- Prior art keywords
- particle
- contact
- particles
- center
- plane
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000002245 particle Substances 0.000 title claims description 488
- 238000004458 analytical method Methods 0.000 title claims description 24
- 239000013598 vector Substances 0.000 claims description 128
- 238000004364 calculation method Methods 0.000 claims description 71
- 238000000034 method Methods 0.000 claims description 51
- 230000010365 information processing Effects 0.000 claims description 25
- 239000000284 extract Substances 0.000 claims description 6
- 238000000605 extraction Methods 0.000 description 15
- 238000012545 processing Methods 0.000 description 14
- 238000007373 indentation Methods 0.000 description 12
- 238000010586 diagram Methods 0.000 description 9
- 230000035515 penetration Effects 0.000 description 8
- 230000008569 process Effects 0.000 description 8
- 238000004140 cleaning Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 230000006870 function Effects 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 230000005484 gravity Effects 0.000 description 4
- 239000000843 powder Substances 0.000 description 4
- 239000012798 spherical particle Substances 0.000 description 4
- 238000011439 discrete element method Methods 0.000 description 3
- 238000013519 translation Methods 0.000 description 3
- 230000001133 acceleration Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 239000000470 constituent Substances 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000013016 damping Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000009472 formulation Methods 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- 238000007792 addition Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 239000008187 granular material Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Description
本発明は、密集状態にある多数の粒子について、その形状を考慮して挙動を解析する技術に関するものであり、特に挙動解析計算における粒子同士の接触判定に関する。 The present invention relates to a technique for analyzing the behavior of a large number of particles in a dense state in consideration of the shape thereof, and more particularly to contact determination between particles in behavior analysis calculation.
粉体や粒体などの粒子を取り扱う分野では、個々の粒子の挙動を知ることが重要である。従来これらの分野では、実験と観察から粒子の挙動を把握することが多かった。しかし近年、計算機内でその運動を再現させ、その挙動を把握する、いわゆるシミュレーション技術の適用が活発に行われている。正確な粒子挙動のシミュレーションを行うには、粒子のもつ形状、機械的特性等の情報をできるだけ詳細に計算モデルに反映させ、それらの因子を考慮して解を求めることが重要である。また、土木分野における砕石、砂利等の比較的大きな物体も「粒子」と考えることにより、上記シミュレーション技術をこれらの物体の挙動を計算する場合へも適用可能である。 In the field of handling particles such as powder and granules, it is important to know the behavior of individual particles. Conventionally, in these fields, the behavior of particles is often grasped from experiments and observations. However, in recent years, so-called simulation technology has been actively applied to reproduce the motion in a computer and grasp its behavior. In order to accurately simulate particle behavior, it is important to reflect information such as the shape and mechanical characteristics of particles in the calculation model in as much detail as possible, and obtain a solution in consideration of those factors. Further, by considering relatively large objects such as crushed stones and gravel in the field of civil engineering as “particles”, the above simulation technique can be applied to the case of calculating the behavior of these objects.
一方、このような粒子を扱う装置の１つに、プリンタ、複写機、ファクシミリ等の電子写真技術を用いた画像形成装置がある。これらの装置では、直径数ミクロンの粉体粒子からなるトナーを用いて、帯電、露光、現像、転写、定着、クリーニングという６つのプロセスを経て、紙上に画像を形成する。このような画像形成装置では、トナー粒子の形状が、作成される画像に大きな影響を与えることがわかっている。そこでこのような画像形成装置の開発においても、シミュレーション技術を用いて、密集状態にある個々のトナー粒子の挙動を求める試みが活発に行われ、近年の傾向として特にトナー粒子の形状を考慮することが求められている。 On the other hand, there is an image forming apparatus using electrophotographic technology such as a printer, a copying machine, and a facsimile as one of apparatuses that handle such particles. In these apparatuses, an image is formed on paper through six processes of charging, exposure, development, transfer, fixing, and cleaning using toner composed of powder particles having a diameter of several microns. In such an image forming apparatus, it is known that the shape of toner particles has a great influence on an image to be created. Therefore, in the development of such an image forming apparatus, an attempt to obtain the behavior of individual toner particles in a dense state using a simulation technique has been actively carried out, and the shape of the toner particles is particularly considered as a recent trend. Is required.
粒子が密集した状態での各粒子の挙動をシミュレーションする代表的な方法として個別要素法（ＤＥＭ）がある。個別要素法の技術は非特許文献１に詳細に記載されているが、個々の粒子を球形で表現し、それらの接触に際してはそのめり込み量からヘルツの式とミンドリンの式を用いることで粒子の持つ機械特性を考慮した接触抗力を求める。そしてニュートンの運動方程式（並進方向）とオイラーの運動方程式（回転方向）を用いて、その接触抗力をもとに個々の粒子の時間進展に伴う運動を求めるものである。
As a typical method for simulating the behavior of each particle in a dense state, there is a discrete element method (DEM). The technique of the individual element method is described in detail in
個別要素法の特徴は個々の粒子を要素と捉え、それらが接触する部分にバネ、ダンパー、スライダを仮定して、その接触力を求め、運動方程式を解くところにある。この考え方を球形以外の形状をもつ粒子の挙動計算に適用したものとして、非特許文献２、非特許文献３がある。これらの文献では粒子の形状を多面体で表現し上述の個別要素法を用いてその挙動を計算している。そして、粒子の接触の検出に関して、Common-planeと呼ばれる面を想定して、接触の有無、接触力の作用点、接触力の作用方向、及び接触力の大きさを決定する。しかしながら、同文献の方法ではCommon-planeと粒子の接触断面積を求めることが必要となるが、断面積を求めるには粒子の各頂点、辺、面とCommon-planeとの接触状態を調べなければならないため、長い計算時間を要する。また本発明者らが同文献の方法を実際に試行しようとしたところ、同文献は接触力を作用させる位置についての情報開示が乏しく、再現することが困難であることが分かった。
The distinct element method is characterized by taking individual particles as elements, assuming springs, dampers, and sliders where they come in contact, finding the contact force, and solving the equation of motion.
本発明の目的は、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することにある。また本発明の目的は、密集状態にある多数の粒子の挙動を計算するに際し、粒子同士の接触判定を高速に行うことにある。 An object of the present invention is to calculate the behavior of a large number of particles in a dense state in consideration of their shapes at high speed. Another object of the present invention is to perform contact determination between particles at high speed when calculating the behavior of a large number of particles in a dense state.
上記課題を解決するために、本発明は以下の構成を採用する。
本発明の第一態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報
処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面を設定し、前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。
In order to solve the above problems, the present invention employs the following configuration.
A first aspect of the present invention is a method for calculating the motion of a plurality of particles in a dense state by an information processing device, and the information processing device acquires data defining the shape of each of the plurality of particles by a polyhedron. An acquisition step; a contact determination step for determining whether or not each particle is in contact with another object based on positional information of each of the plurality of particles; and an information processing device for determining the contact A calculation step of calculating a force received by each particle from another object determined to be in contact in the step, and calculating a motion of each particle by the force; and an information processing device, each of which is calculated in the calculation step An output step for outputting the motion of the particles, and in the contact determination step, when determining whether or not the particle i is in contact with the particle j, the center of the particle i and the center of the particle j are connected. Sets a first plane tangent to and particle j intersect a straight line, to determine the unit vector n1 facing the center side of and particle j perpendicular to said first plane, the position of each vertex of the unit vectors n1 and particles i The inner product n1 · Pi of the vector Pi is calculated, the inner product n1 · Pj of the unit vector n1 and the position vector Pj of each vertex of the particle j is calculated, and the maximum of the inner products n1 · Pi for all the vertices of the particle i is And the inner product n1 · Pj for all the vertices of the particle j, where the smallest one is D1min, the vertex of the particle i where the inner product n1 · Pi is D1min or more and D1max or less, and the inner product n1 · Pj is to define the set a consisting of the vertex of the particle j to be less D1max more D1min, extracts the face side and the particle j of particle i including vertices in the set a, the extract A particle behavior analysis method characterized by performing the contact determination in between the edges and faces of the particle j particle i.
本発明の第二態様は、情報処理装置により、密集状態にある複数の粒子の運動を計算する方法であって、情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析方法である。 The second aspect of the present invention is a method for calculating the motion of a plurality of particles in a dense state by an information processing device, wherein the information processing device acquires data defining the shape of each of the plurality of particles by a polyhedron. An acquisition step; a contact determination step for determining whether or not each particle is in contact with another object based on positional information of each of the plurality of particles; and an information processing device for determining the contact A calculation step of calculating a force received by each particle from another object determined to be in contact in the step, and calculating a motion of each particle by the force; and an information processing device, each of which is calculated in the calculation step An output step for outputting the motion of the particles, and in the contact determination step, when determining whether or not the particle i is in contact with the particle j, the center of the particle i and the center of the particle j are connected. More than a third plane that intersects the straight line and is closer to the center of the particle j than the first plane that contacts the particle j and that intersects the straight line and contacts the particle j at an angle different from the first plane. The side of the particle i located in the space on the center side of the particle j, the center side of the particle i from the second plane parallel to the first plane and in contact with the particle i, and the third The surface of the particle j located in the space on the center side of the particle i with respect to the fourth plane in contact with the particle i and between the side of the extracted particle i and the surface of the particle j is extracted. This is a particle behavior analysis method characterized in that contact determination is performed by using a method.
本発明の第三態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面を設定し、前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。 A third aspect of the present invention is a particle behavior analysis device that calculates the motion of a plurality of particles in a dense state, an acquisition unit that acquires data defining the shape of each of the plurality of particles by a polyhedron, and the plurality of the plurality of particles Based on the position information of each particle, a contact determination unit that determines whether or not each particle is in contact with another object, and each particle from another object that is determined to be in contact with the contact determination unit A calculation unit that calculates a force to be received and calculates a motion of each particle by the force; and a calculation result output unit that outputs a motion of each particle calculated by the calculation unit; When determining whether i is in contact with the particle j, a first plane that intersects the straight line connecting the center of the particle i and the center of the particle j and is in contact with the particle j is set. A unit vector n perpendicular to the center of the particle j , The inner product n1 · Pj of the unit vector n1 and the position vector Pi of each vertex of the particle i is calculated, the inner product n1 · Pj of the unit vector n1 and the position vector Pj of each vertex of the particle j is calculated, and the particle i Of the inner products n1 · Pi for all the vertices of the particle j, and the largest one of the inner products n1 · Pj for all the vertices of the particle j is D1min. Define a set A consisting of the vertices of particles i for which D1min is not less than D1max and the vertices of particles j for which the inner product n1 · Pj is not less than D1min and not more than D1max. The particle behavior analysis apparatus is characterized in that a surface of a particle j is extracted and contact determination is performed between the side of the extracted particle i and the surface of the particle j.
本発明の第四態様は、密集状態にある複数の粒子の運動を計算する粒子挙動解析装置であって、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行うことを特徴とする粒子挙動解析装置である。 A fourth aspect of the present invention is a particle behavior analyzer that calculates the motion of a plurality of particles in a dense state, an acquisition unit that acquires data defining the shape of each of the plurality of particles by a polyhedron, and the plurality of the plurality of particles Based on the position information of each particle, a contact determination unit that determines whether or not each particle is in contact with another object, and each particle from another object that is determined to be in contact with the contact determination unit A calculation unit that calculates a force to be received and calculates a motion of each particle by the force; and a calculation result output unit that outputs a motion of each particle calculated by the calculation unit; When determining whether or not i is in contact with the particle j, it intersects a straight line connecting the center of the particle i and the center of the particle j and is closer to the center of the particle j than the first plane in contact with the particle j. And at an angle different from the first plane The side of the particle i located in the space closer to the center of the particle j than the third plane intersecting the line and contacting the particle j, and the second plane parallel to the first plane and contacting the particle i Extracting the surface of the particle j located in the center side of the particle i and in the space on the center side of the particle i with respect to the fourth plane parallel to the third plane and in contact with the particle i; The particle behavior analysis apparatus is characterized in that contact determination is performed between the side of the extracted particle i and the surface of the particle j.
本発明の第五態様は、上述した粒子挙動解析方法の各ステップを情報処理装置（コンピュータ）に実行させる解析プログラムである。 A fifth aspect of the present invention is an analysis program that causes an information processing apparatus (computer) to execute each step of the particle behavior analysis method described above.
本発明によれば、粒子同士の接触判定に要する時間を短縮でき、密集状態にある多数の粒子の挙動を、その形状も考慮して高速に計算することができる。 According to the present invention, the time required for determining contact between particles can be shortened, and the behavior of many particles in a dense state can be calculated at high speed in consideration of their shapes.
本発明の実施形態は、密集状態にある多数の粒子の挙動をコンピュータ・シミュレーションにより計算するための粒子挙動解析方法及び装置を提供するものである。本実施形態に係る粒子挙動解析方法の特徴の一つは、多数の多角形から構成される多面体（ポリゴンモデル）により各粒子の形状を定義した点である。これにより任意の粒子形状を表現でき
るようになり、多様なケースのシミュレーションを精度良く行うことが可能になる。本方法のもう一つの特徴は、粒子同士の接触判定を行う際に、全ての頂点、辺、面について網羅的に接触判定計算を行うのではなく、予め接触可能性のある辺及び面を特定し、それらのあいだについてのみ接触判定を行う点である。これにより、接触判定に要する時間を短縮でき、粒子の挙動計算の高速化を図ることができる。以下、本実施形態の具体的な構成及び特徴について図に基づき詳しく説明する。
Embodiments of the present invention provide a particle behavior analysis method and apparatus for calculating the behavior of a large number of particles in a dense state by computer simulation. One of the features of the particle behavior analysis method according to the present embodiment is that the shape of each particle is defined by a polyhedron (polygon model) composed of a large number of polygons. As a result, an arbitrary particle shape can be expressed, and various cases of simulation can be accurately performed. Another feature of this method is that when performing contact determination between particles, instead of exhaustively performing contact determination calculations for all vertices, sides, and surfaces, the sides and surfaces that may be contacted are identified in advance. However, the contact determination is performed only between them. As a result, the time required for contact determination can be shortened, and the speed of particle behavior calculation can be increased. Hereinafter, specific configurations and features of the present embodiment will be described in detail with reference to the drawings.
（粒子挙動計算のモデル）
図１は、本実施形態における粒子挙動計算の力学モデルを説明するものであり、異形粒子（非球形粒子）が他の物体と接触した際に働く接触力を説明する図である。このモデルでは、粒子が他の物体と接触した際に発生するめり込みを、その接触面に対して法線方向と接線方向に分解して考える。粒子には、それぞれのめり込み量に対してスプリングで押し返す力が働く。また各スプリングと並列にダンパーを設定する。これはめり込み速度に比例して働く力であり、運動を減衰させるものである。またすべりを考慮するためのスライダを設定している。このように、スプリング、ダンパー、スライダを設定して衝突を考えることは従来の個別要素法と同じである。
(Particle behavior calculation model)
FIG. 1 is a diagram for explaining a dynamic model of particle behavior calculation in the present embodiment, and is a diagram for explaining a contact force that acts when an irregularly shaped particle (non-spherical particle) comes into contact with another object. In this model, the indentation that occurs when a particle comes into contact with another object is considered by decomposing it into a normal direction and a tangential direction with respect to the contact surface. The particles have a force that pushes them back with a spring for each amount of penetration. A damper is set in parallel with each spring. This is a force that works in proportion to the penetration speed and attenuates the movement. In addition, a slider is set to take slip into consideration. Thus, setting a spring, a damper, and a slider to consider a collision is the same as the conventional individual element method.
（解析装置のハードウエア構成）
図２は、本実施形態における解析装置が動作する計算装置を示すブロック図である。情報処理装置（コンピュータ）１０は、各種判断及び処理を行う中央処理装置（ＣＰＵ）１１と、処理プログラム及び固定データを格納するＲＯＭ１２と、処理データを記憶するＲＡＭ１３と、入出力回路（Ｉ／Ｏ）１４とから構成されている。また必要に応じて、情報処理装置１０には、Ｉ／Ｏ１４を介して記憶装置、入力装置、表示装置、他のコンピュータなどが接続される。Ｉ／Ｏ１４は、これら外部装置と情報処理装置１０との間でデータをやり取りするための回路である。この情報処理装置１０には、入力データ１５がＩ／Ｏ１４を介して入力される。計算結果は、Ｉ／Ｏ１４を介して出力データ１６として出力される。
(Hardware configuration of analysis device)
FIG. 2 is a block diagram illustrating a computing device on which the analysis device according to the present embodiment operates. An information processing apparatus (computer) 10 includes a central processing unit (CPU) 11 that performs various determinations and processes, a
（解析装置の機能）
図３は、粒子挙動解析を行う解析プログラムによって実現される機能の構成を示すブロック図である。解析プログラムは、ＲＯＭ１２や記憶装置などのコンピュータ読み取り可能な記録媒体に記憶されており、必要に応じて、ＣＰＵ１１によって読み込まれ実行されるソフトウエアである。この解析プログラムは、主な機能（モジュール）として、入力データ読込部（取得部）Ｂ１１０、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、粒子運動計算部Ｂ１４０、計算結果出力部Ｂ１５０を備えている。全体制御部Ｂ１００は、これらのモジュールの処理の順番やデータの引渡しなどの全体制御を担っている。
(Function of analysis device)
FIG. 3 is a block diagram illustrating a configuration of functions realized by an analysis program that performs particle behavior analysis. The analysis program is software that is stored in a computer-readable recording medium such as the
入力データ読込部Ｂ１１０は、粒子形状データ、壁形状データ、及び各種パラメータを読み込んでＲＡＭ１３に格納する。粒子形状データとは、粒子の表面を複数の多角形（ポリゴン）で表現したものである。ここで多角形としてはその頂点が一平面上にあることが望ましく、３角形を用いるのが最も単純である。また粒子は凸形状（凸多面体）に限定すると、他の粒子や壁との接触状態の判定が簡便になる。壁形状データも同様であり、壁の表面を複数の多角形で表現したものである。粒子形状データ、壁形状データともに、多角形の頂点座標と、各多角形を構成する頂点の番号とで構成されるデータである。各種パラメータには、粒子と壁の機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）、粒子に働く力（重力、付着力、電磁力、空気の粘性特性等）、及び計算条件としての計算刻み時間（Δｔ）、計算終了時刻がある。なお、壁とは、粒子が移動可能な空間を形成する（規定する）境界面である。
The input data reading unit B110 reads particle shape data, wall shape data, and various parameters and stores them in the
接触判定部Ｂ１２０は、粒子の位置情報（各頂点の座標、粒子の中心座標）に基づいて
、粒子が他の物体（周囲の粒子または壁）と接触しているかどうか、またどのように接触しているかを判定する。粒子受力計算部Ｂ１３０は、各粒子が他の物体から受ける力を計算する機能である。粒子は他の粒子または壁と接触している場合に、そのめり込み量に応じて接触抗力を受ける。粒子運動計算部Ｂ１４０は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。計算結果出力部Ｂ１５０は、得られた各粒子の位置と速度の結果を出力データ１６として出力するものである。
The contact determination unit B120 determines whether and how the particles are in contact with other objects (the surrounding particles or walls) based on the position information of the particles (the coordinates of each vertex and the center coordinates of the particles). Judge whether it is. The particle receiving force calculation unit B130 has a function of calculating the force that each particle receives from another object. When a particle is in contact with other particles or walls, it receives contact drag depending on the amount of penetration. The particle motion calculation unit B140 obtains the speed and position of each particle after the ticking time due to the force acting on the particles obtained by the particle force calculation unit B130. The calculation result output unit B150 outputs the obtained position and velocity results of each particle as
以下に、接触判定部Ｂ１２０、粒子受力計算部Ｂ１３０、及び粒子運動計算部Ｂ１４０の詳しい内容と計算方法を示す。 Below, the detailed content and calculation method of contact determination part B120, particle | grain receiving power calculation part B130, and particle motion calculation part B140 are shown.
（接触判定部Ｂ１２０）
まずは、接触判定部Ｂ１２０における粒子同士の接触判定の内容を図４を用いて説明する。図４の上段は４種類の接触パターンを３次元的に示したものであり、下段は各パターンにおける２つの粒子の接触によるめり込みの様子を示すものである。図４において、実線は６面体粒子ｉ、破線は６面体粒子ｊを示す。ｉ−ｖ０は粒子ｉを構成する多角形の頂点、ｉ−ｅ０は粒子ｉを構成する多角形の辺、ｊ−ｅ０とｊ−ｅ１は粒子ｊを構成する多角形の辺、ｊ−ｆ０は粒子ｊを構成する多角形の面、星印は粒子ｉと粒子ｊの接触点を示す。また図４の上段における矢印は粒子ｉが粒子ｊに衝突する方向を示し、また下段における黒丸は粒子ｊの辺を示す。粒子が凸形状である場合、粒子ｉと粒子ｊの接触は、図４に示す次の（ａ）〜（ｄ）の４つのパターンのいずれかとなる。なお、図４では、粒子同士の接触判定を例に挙げているが、粒子と壁の接触判定も同じように考えることができる。
(Contact determination unit B120)
First, the content of the contact determination between the particles in the contact determination unit B120 will be described with reference to FIG. The upper part of FIG. 4 shows four types of contact patterns in a three-dimensional manner, and the lower part shows a state of indentation due to contact of two particles in each pattern. In FIG. 4, the solid line indicates hexahedral particles i, and the broken line indicates hexahedral particles j. i-v0 is the vertex of the polygon constituting the particle i, i-e0 is the side of the polygon constituting the particle i, j-e0 and j-e1 are the sides of the polygon constituting the particle j, and j-f0 is A polygonal surface constituting the particle j, an asterisk, indicates a contact point between the particle i and the particle j. In addition, the arrow in the upper part of FIG. 4 indicates the direction in which the particle i collides with the particle j, and the black circle in the lower part indicates the side of the particle j. When the particles have a convex shape, the contact between the particles i and the particles j is any one of the following four patterns (a) to (d) shown in FIG. In addition, in FIG. 4, although the contact determination of particle | grains is mentioned as an example, the contact determination of particle | grains and a wall can be considered similarly.
（ａ）頂点のめり込み：
粒子ｉを構成する多角形の頂点が粒子ｊの１つの面にめり込む場合である。粒子ｉの頂点が粒子ｊの内部にあれば、接触していると判定する。図４（ａ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込んでいる。このパターンでは、頂点ｉ−ｖ０を、粒子ｉの粒子ｊに対する接触点として考える。
(A) Vertex penetration:
This is a case where the vertexes of the polygon constituting the particle i are embedded in one surface of the particle j. If the vertex of the particle i is inside the particle j, it is determined that they are in contact. In FIG. 4A, the vertex i-v0 of the particle i is embedded in the surface j-f0 of the particle j. In this pattern, the vertex i-v0 is considered as a contact point of the particle i with respect to the particle j.
（ｂ）辺のめり込み（ケース１）：
上記（ａ）に加えて、めり込んだ頂点を端点に持つ粒子ｉの辺が粒子ｊの辺にもめり込む場合である。（ａ）のめり込んだ頂点を端点に持つ粒子ｉの辺において、粒子ｊ構成面との交点数が１であり、かつその交点を持つ面が（ａ）におけるめり込み面でなければ、これに該当する。図４（ｂ）においては、粒子ｉの頂点ｉ−ｖ０が粒子ｊの面ｊ−ｆ０にめり込む（上記（ａ）に該当）とともに、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にもめり込んでいる。このパターンでは、頂点ｉ−ｖ０に加え、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点も、粒子ｉの粒子ｊに対する接触点として考える。
(B) Edge penetration (case 1):
In addition to the above (a), the side of the particle i having an indented vertex is also indented into the side of the particle j. This corresponds to the side of the particle i whose end point is the recessed vertex of (a) and the number of intersections with the surface constituting the particle j is 1, and the surface having the intersection is not the recessed surface in (a). . In FIG. 4B, the vertex i-v0 of the particle i sinks into the surface j-f0 of the particle j (corresponding to the above (a)), and the side i-e0 of the particle i is the side j-e0 of the particle j. I'm also into it. In this pattern, in addition to the vertex i-v0, a point on the side i-e0 closest to the side j-e0 is also considered as a contact point of the particle i with respect to the particle j.
（ｃ）辺のめり込み（ケース２）：
粒子ｉの辺が粒子ｊの１つの辺だけにめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、その交点を持つ２つの面が共有辺を持ち、かつその共有辺が粒子ｉの辺に近接していれば、これに該当する。図４（ｃ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０にめり込んでいる。このパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
(C) Edge indentation (Case 2):
In this case, the side of the particle i is sunk into only one side of the particle j. In the side i-e0 of the particle i, the number of intersections with the particle j constituting surface is 2, the two surfaces having the intersection have a shared side, and the shared side is close to the side of the particle i. This is the case. In FIG. 4C, the side i-e0 of the particle i is recessed into the side j-e0 of the particle j. In this pattern, a point on the side i-e0 closest to the side j-e0 is considered as a contact point of the particle i with respect to the particle j.
（ｄ）辺のめり込み（ケース３）：
粒子ｉの辺が粒子ｊの２つの辺にめり込む場合である。粒子ｉの辺ｉ−ｅ０において、粒子ｊ構成面との交点数が２であり、（ｃ）ではない場合がこれに該当する。図４（ｄ）においては、粒子ｉの辺ｉ−ｅ０が粒子ｊの辺ｊ−ｅ０とｊ−ｅ１にめり込んでいる。こ
のパターンでは、辺ｊ−ｅ０に最も近い辺ｉ−ｅ０上の点、及び、辺ｊ−ｅ１に最も近い辺ｉ−ｅ０上の点を、粒子ｉの粒子ｊに対する接触点として考える。
(D) Edge indentation (Case 3):
This is a case where the side of the particle i sinks into two sides of the particle j. This corresponds to the case where the number of intersections with the particle j constituting surface is 2 in the side i-e0 of the particle i and is not (c). In FIG. 4D, the side i-e0 of the particle i is embedded in the sides j-e0 and j-e1 of the particle j. In this pattern, a point on the side i-e0 closest to the side j-e0 and a point on the side i-e0 closest to the side j-e1 are considered as contact points of the particle i with respect to the particle j.
上記の接触判定を行うには、図４（ａ）の場合には、粒子ｉの頂点と粒子ｊの面の位置関係、図４（ｂ）〜（ｄ）においては、粒子ｉの辺と粒子ｊの面の交差を調べる必要がある。辺と面のめり込み量の計算を行うための接触判定の仕方を６面体の粒子を例に、図５を用いて説明する。 In the case of FIG. 4A, the above contact determination is performed by the positional relationship between the vertex of the particle i and the surface of the particle j, and in FIGS. 4B to 4D, the side of the particle i and the particle It is necessary to examine the intersection of j planes. A method of contact determination for calculating the amount of indentation between sides and faces will be described with reference to FIG. 5, taking hexahedral particles as an example.
図５の２０は粒子ｊを構成する１つの多角形の面を表し、この法線ベクトルをｎとする。２１は粒子ｉを構成する辺の１つｉ−ｅ０を表している。この辺２１の両端の頂点をそれぞれＰ０（ｐｘ０，ｐｙ０，ｐｚ０），Ｐ１（ｐｘ１，ｐｙ１，ｐｚ１）とする。また、頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）と辺２１の交差点を点Ｑとする。このとき、点Ｑが多角形面２０上にあるか否かの判定を行う。
5 in FIG. 5 represents one polygonal surface constituting the particle j, and this normal vector is n.
まず、粒子ｉの辺２１が粒子ｊの多角形面２０の４つの頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を含む平面（無限）を横切っているかどうかを調べる。
First, it is examined whether the
粒子ｊの多角形面２０上の一つの任意点Ａ０と粒子ｉの辺２１の頂点Ｐ０，Ｐ１をそれぞれ結び、ベクトルＰ０Ａ０，Ｐ１Ａ０とする。ベクトルＰ０Ａ０は、点Ｐ０から点Ａ０に向かうベクトルであり、ベクトルＰ１Ａ０は、点Ｐ１から点Ａ０に向かうベクトルである。
One arbitrary point A0 on the
このとき、２つのベクトルＰ０Ａ０，Ｐ１Ａ０と法線ベクトルｎの内積Ｌ０＝Ｐ０Ａ０・ｎとＬ１＝Ｐ１Ａ０・ｎが同じ符号ならば、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていないことがわかる。逆に、内積Ｌ０とＬ１の符号が異なる場合は、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っていることがわかる。
At this time, if the inner product L0 = P0A0 · n and L1 = P1A0 · n of the two vectors P0A0 and P1A0 and the normal vector n have the same sign, the
次に、粒子ｉの辺２１が粒子ｊの多角形面２０の頂点を含む平面（無限）を横切っている場合、粒子ｊの多角形面２０の頂点を含む平面（無限）と粒子ｉの辺２１との交点Ｑの座標（ｑｘ，ｑｙ，ｑｚ）を下記式で求める。
ここで、ｒは内分比であり、ｒ＝Ｌ０／（Ｌ０−Ｌ１）で求まる。
Next, when the
Here, r is an internal ratio, and is obtained by r = L0 / (L0−L1).
次に、交点Ｑと頂点ｊ−ｖ０，ｊ−ｖ１，ｊ−ｖ２，ｊ−ｖ３を結ぶ４つのベクトルＱｖ０，Ｑｖ１，Ｑｖ２，Ｑｖ３について、隣のベクトルとの外積を下記式を用いて求め、ベクトルＶ０，Ｖ１，Ｖ２，Ｖ３を得る。
ベクトルＶ０は、ベクトルＱｖ０とＱｖ１の両方に垂直なベクトル、すなわち多角形面２０に対し垂直なベクトルとなり、その向きはベクトルＱｖ０とＱｖ１のなす角が０〜πかπ〜２πかで反転する。他のベクトルＶ１〜Ｖ３も同様の性質をもつ。交点Ｑが多角形面２０の内部にある場合は、Ｑｖ０とＱｖ１、Ｑｖ１とＱｖ２、Ｑｖ２とＱｖ３、Ｑｖ３とＱｖ０のいずれのなす角も０〜πとなるから、ベクトルＶ０〜Ｖ３の向きは一致する。一方、交点Ｑが多角形面２０の外部にある場合は、なす角が０〜πとなる組と、π〜２πとなる組とがでるため、ベクトルＶ０〜Ｖ３の向きは一致しなくなる。このような性質を利用することで、交点Ｑが多角形面２０の内側に存在するか否か、つまり粒子ｉの辺２１と粒子ｊの多角形面２０とが交差するか否かを判定することができる。
The vector V0 is a vector perpendicular to both the vectors Qv0 and Qv1, that is, a vector perpendicular to the
ここでは、下記式のように内積Ｍ０〜Ｍ３を計算し、それらの符号を調べることで判定を行う。
Ｍ０〜Ｍ３の符号が全て同じ場合は、粒子ｊの多角形面２０と粒子ｉの辺２１は交差しており、一つでも符号が逆のものがあれば交差していないことがわかる。
Here, determination is performed by calculating inner products M0 to M3 as in the following equation and examining their signs.
When all the signs of M0 to M3 are the same, the
上記のような接触判定では、１つの四角形面に対して１つの辺が交差していることを判定するのに、４３個の加減算、４７個の乗算、５回の条件演算を要する。そのため、粒子の多面体の面数が多いほど、また粒子の数が多いほど、膨大な計算時間がかかってしまう。そこで、本実施形態では接触判定の前に接触可能頂点を絞り込み、その接触可能頂点を含む辺、面のみを接触の可能性がある辺、面として接触判定を行う。以下、接触の可能性がある辺及び面を抽出する方法として、２つの方法を説明する。 In the contact determination as described above, 43 additions / subtractions, 47 multiplications, and 5 conditional operations are required to determine that one side intersects one rectangular plane. For this reason, the larger the number of polyhedrons of particles and the larger the number of particles, the longer the calculation time. Therefore, in this embodiment, contactable vertices are narrowed down before contact determination, and contact determination is performed using only the sides and surfaces including the contactable vertices as possible sides and surfaces. Hereinafter, two methods will be described as methods for extracting sides and surfaces that may be contacted.
（第１の抽出方法）
図６は、辺及び面の第１の抽出方法について説明する図である。粒子ｉ、ｊは３次元形状であるが、簡単のため２次元で表示している。
(First extraction method)
FIG. 6 is a diagram for explaining a first extraction method of sides and faces. Particles i and j have a three-dimensional shape, but are displayed in two dimensions for simplicity.
第１の抽出方法は、粒子ｉの粒子ｊに対する接触を判定する際に、粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わる２つの平面６１，６２を設定し、平面６１，６２で挟まれた空間に位置する粒子ｉの辺及び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面６１（第１の平面）と平面６２（第２の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第１の抽出方法は、平面６１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。なお、辺又は面が「空間に位置する」とは、辺又は面の少なくとも一部分が空間内に存在すること（空間と交わりを有すること）を意味する。
In the first extraction method, when determining the contact of the particle i with the particle j, two
第１の抽出方法を実現する具体的なアルゴリズムとして、本実施形態の接触判定部Ｂ１２０では次のような計算を行う。粒子ｉの中心から粒子ｊの中心に向かう中心結線ベクトルを考え、原点を通り中心結線ベクトルに垂直な面を基準面Ｓとする。粒子ｉの各頂点と基準面Ｓとの距離Ｄ１を計算し、最大値をＤ１ｍａｘとする。粒子ｊの各頂点と基準面Ｓとの距離を計算し、最小値をＤ１ｍｉｎとする。基準面Ｓとの距離がＤ１ｍｉｎ以上Ｄ１
ｍａｘ以下となる粒子ｉ、ｊの各頂点を接触可能頂点として選ぶ。図６で塗りつぶし（黒丸）で示されている頂点が接触可能頂点である。そして、接触可能頂点を含む辺及び面が判定対象として抽出される。図６で実線で示されている辺が抽出された辺である。
As a specific algorithm for realizing the first extraction method, the contact determination unit B120 of the present embodiment performs the following calculation. A center connection vector from the center of the particle i to the center of the particle j is considered, and a plane passing through the origin and perpendicular to the center connection vector is defined as a reference plane S. A distance D1 between each vertex of the particle i and the reference plane S is calculated, and the maximum value is set to D1max. The distance between each vertex of the particle j and the reference plane S is calculated, and the minimum value is D1min. The distance from the reference surface S is D1 min or more D1
Each vertex of the particles i and j that is equal to or less than max is selected as a contactable vertex. In FIG. 6, the vertexes indicated by filling (black circles) are contactable vertices. Then, sides and surfaces including the contactable vertices are extracted as determination targets. The sides indicated by solid lines in FIG. 6 are the extracted sides.
Ｄ１ｍａｘ、Ｄ１ｍｉｎの算出方法をより詳しく説明する。中心結線ベクトルと向きが同じで原点を起点とする単位ベクトルをｎ１＝（ｎ１ｘ，ｎ１ｙ，ｎ１ｚ）とし、粒子ｉの各頂点の位置ベクトルをＰｉｋ＝（ｘｉｋ，ｙｉｋ，ｚｉｋ）とする。ｎ１のことを中心結線単位ベクトルとも称す。ここで、ｋは粒子ｉの各構成頂点の番号である。基準面Ｓと粒子ｉの各頂点との距離Ｄ１は、下記式のように、単位ベクトルｎ１と頂点の位置ベクトルＰｉｋの内積ｎ１・Ｐｉｋで算出される。
同様に、粒子ｊの各頂点の位置ベクトルをＰｊｋ＝（ｘｊｋ，ｙｊｋ，ｚｊｋ）とすると、基準面Ｓと粒子ｊの各頂点との距離Ｄ１は下記式により算出される。
Similarly, when the position vector of each vertex of the particle j is Pjk = (x jk , y jk , z jk ), the distance D1 between the reference plane S and each vertex of the particle j is calculated by the following equation.
粒子ｉ及び粒子ｊの全ての頂点について距離Ｄ１を求めた後、Ｄ１ｍａｘ及びＤ１ｍｉｎは下記式で算出される。
なお、説明を簡単にするためＤ１を距離と表現したが、実際にはＤ１は内積ｎ１・Ｐｉｋ、内積ｎ１・Ｐｊｋのようにベクトルの内積で定義されるものである。よって、ｎ１とＰｉｋ又はＰｊｋの位置関係によっては、各頂点に対するＤ１の一部またはすべてが負になることもあるが、その場合でも上記式によりＤ１ｍａｘ、Ｄ１ｍｉｎを求めることができる。また、ここでは簡単のためｎ１の起点を原点としたが、ベクトルの平行移動に対して内積の大小は変わらないため、起点をどの位置にとってもよい。また、ここでは単位ベクトルｎ１を中心結線ベクトルに平行にとったが、ベクトルｎ１の向きはこれに限られない。ベクトルｎ１が中心結線ベクトルに対して直交していなければ（つまり、平面６１，６２が粒子ｉ，ｊの中心を結んだ直線に交わっていれば）、同様の処理が可能である。
In order to simplify the description, D1 is expressed as a distance. Actually, D1 is defined as an inner product of vectors such as inner product n1 · Pik and inner product n1 · Pjk. Therefore, depending on the positional relationship between n1 and Pik or Pjk, part or all of D1 with respect to each vertex may be negative. Even in this case, D1max and D1min can be obtained by the above equations. For simplicity, the starting point of n1 is used as the origin. However, the magnitude of the inner product does not change with respect to the parallel movement of the vector, and the starting point may be at any position. Here, the unit vector n1 is parallel to the central connection vector, but the direction of the vector n1 is not limited to this. If the vector n1 is not orthogonal to the center connection vector (that is, if the
以上述べた方法により接触可能頂点を選択し、判定対象とする粒子ｉの辺と粒子ｊの面を絞り込むことで、接触判定を行う辺、面の組み合わせを減らし、計算時間を短縮することができる。 By selecting contactable vertices by the method described above and narrowing down the sides of the particle i and the particle j to be determined, the number of combinations of sides and surfaces for contact determination can be reduced, and the calculation time can be shortened. .
（第２の抽出方法）
上述した第１の抽出方法では、接触判定を行う範囲を平面６１と６２で挟まれた空間に限定している。しかし、図６で示されるように、平面６１と６２のあいだに、実際には他の粒子と接触していない頂点が多く含まれる場合もある。そこで、第２の抽出方法として、接触可能頂点をさらに絞り込む方法について説明する。
(Second extraction method)
In the first extraction method described above, the range in which the contact determination is performed is limited to a space sandwiched between the
第２の抽出方法は、図７に示すように、平面６１，６２とは異なる角度で平面７１，７２を設定し、４つの平面６１，６２，７１，７２で挟まれた空間に位置する粒子ｉの辺及
び粒子ｊの面のみを判定対象として抽出するものである。ここで、平面７１（第３の平面）と平面７２（第４の平面）はそれぞれ粒子ｊ，ｉに接し、且つ、互いに平行になるように設定される。言い換えると、第２の抽出方法は、平面６１及び平面７１よりも粒子ｊの中心側の空間に位置する粒子ｉの辺と、平面６２及び平面７２よりも粒子ｉの中心側の空間に位置する粒子ｊの面とを、判定対象として抽出する方法である。
In the second extraction method, as shown in FIG. 7, the
第２の抽出方法を実現する具体的なアルゴリズムを説明する。図７は２つのベクトルによる接触可能頂点の絞り込み方法を示す図である。図７左上は上記の中心結線単位ベクトルｎ１によって絞り込まれた接触可能頂点を示している。塗りつぶし（黒丸）で示されているのが接触可能頂点である。 A specific algorithm for realizing the second extraction method will be described. FIG. 7 is a diagram showing a method of narrowing down touchable vertices using two vectors. The upper left of FIG. 7 shows the contactable vertices narrowed down by the central connection unit vector n1. The vertices that can be touched are indicated by solid fills (black circles).
図７左下は、重なり最小単位ベクトルｎ２によって絞り込まれた接触可能頂点を表している。「重なり最小単位ベクトル」とは、粒子ｉと粒子ｊの重なりを最小とする向きの単位ベクトルをいい、以降の説明では単に「重なり最小ベクトル」とも称す。ｎ２による接触可能頂点の絞り込み方法はｎ１による接触可能頂点の絞り込み方法と同じである。つまり、ベクトルｎ２と粒子ｉの各頂点座標の位置ベクトルＰｉｋの内積ｎ２・Ｐｉｋの最大値Ｄ２ｍａｘと、ベクトルｎ２と粒子ｊの各頂点座標の位置ベクトルＰｊｋの内積ｎ２・Ｐｊｋの最小値Ｄ２ｍｉｎを求める。そして、内積がＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる頂点を接触可能頂点として抽出する。 The lower left of FIG. 7 represents the contactable vertices narrowed down by the overlapping minimum unit vector n2. The “minimum overlap unit vector” refers to a unit vector in a direction that minimizes the overlap between the particle i and the particle j, and is simply referred to as “minimum overlap vector” in the following description. The method for narrowing contactable vertices by n2 is the same as the method for narrowing contactable vertices by n1. That is, the maximum value D2max of the inner product n2 · Pik of the vector n2 and the position vector Pik of each vertex coordinate of the particle i and the minimum value D2min of the inner product n2 · Pjk of the position vector Pjk of the vertex coordinate of the vector n2 and particle j are obtained. . Then, a vertex whose inner product is D2min or more and D2max or less is extracted as a contactable vertex.
図７右は中心結線ベクトルによる接触可能頂点集合をＡ、重なり最小ベクトルによる接触可能頂点集合をＢとしたときに、共通部分である積集合Ａ∩Ｂで選択される頂点を表している。図７において、集合Ａに含まれる頂点数は１０、集合Ａ内の頂点を構成頂点とする辺数は１２、集合Ｂに含まれる頂点数は８、辺数は１０である。一方、積集合Ａ∩Ｂに含まれる頂点数は４、辺数は６となり、１つのベクトルのみの場合より判定対象を絞り込むことができる。よって、積集合Ａ∩Ｂに含まれる頂点を最終的な接触可能頂点と定義する。 The right side of FIG. 7 represents vertices selected by a product set A∩B, which is a common part, where A is a contactable vertex set based on a central connection vector and B is a contactable vertex set based on a minimum overlapping vector. In FIG. 7, the number of vertices included in the set A is 10, the number of sides having the vertices in the set A as constituent vertices is 12, the number of vertices included in the set B is 8, and the number of sides is 10. On the other hand, the number of vertices included in the product set A∩B is 4 and the number of sides is 6, so that the determination target can be narrowed down compared to the case of only one vector. Therefore, vertices included in the product set A∩B are defined as final contactable vertices.
積集合Ａ∩Ｂが空集合でも接触する場合がある。図８を用いて説明する。四角で表されている頂点が中心結線ベクトルにより選択された頂点（集合Ａ）、三角で表されているのが重なり最小ベクトルにより選択された頂点（集合Ｂ）である。塗りつぶしの円で表されているのが積集合Ａ∩Ｂに含まれる接触可能頂点である。粒子ｊは接触可能頂点を持たないが、図８に示すように粒子ｉと接触している。このような場合を含めるため、判定対象とする辺、面は次のように決定する。
・積集合Ａ∩Ｂに含まれる頂点（接触可能頂点）を持つ辺、面。
・集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を持つ辺、面。
There is a case where the product set A∩B is in contact with the empty set. This will be described with reference to FIG. Vertices represented by squares are vertices selected by the center connection vector (set A), and triangles are vertices selected by the overlapping minimum vector (set B). The contactable vertices included in the product set A∩B are represented by filled circles. The particle j does not have a contactable vertex, but is in contact with the particle i as shown in FIG. In order to include such a case, the sides and surfaces to be determined are determined as follows.
-Edges and faces having vertices (contactable vertices) included in the product set A∩B.
Edges and faces that have both vertices included in set A and vertices included in set B.
以上の２パターンによって、全ての接触の可能性がある辺と面を選択することができる。本方法により、接触判定を行う辺、面の組み合わせを大幅に減らすことができ、計算時間のより一層の短縮が可能となる。 By the above two patterns, it is possible to select a side and a face that can be contacted. With this method, the combination of sides and surfaces for performing contact determination can be greatly reduced, and the calculation time can be further shortened.
本方法で用いる重なり最小ベクトルｎ２の算出方法について説明する。図９は重なり最小ベクトルを求めるフローチャートである。図９に示す処理は、接触判定部Ｂ１２０が実行するものである。 A method for calculating the minimum overlap vector n2 used in this method will be described. FIG. 9 is a flowchart for obtaining the minimum overlap vector. The process shown in FIG. 9 is executed by the contact determination unit B120.
まず、重なり最小ベクトルの初期値を設定する（ステップＳ１）。初期値は中心結線単位ベクトルと平行な単位ベクトルとする。粒子ｉの中心位置ベクトルをＣｉ、粒子ｊの中心位置ベクトルをＣｊとすると、重なり最小ベクトルｎ２の初期値は下記式で表される。
以降の処理で、ベクトルｎ２の向きを微小に変動させながら、粒子ｉと粒子ｊの重なりが最小となるものを見つける。その処理で用いる微小量δの初期値δｍａｘ（＝０．１程度）をステップＳ１で設定する。
First, an initial value of a minimum overlap vector is set (step S1). The initial value is a unit vector parallel to the central connection unit vector. If the center position vector of the particle i is Ci and the center position vector of the particle j is Cj, the initial value of the overlap minimum vector n2 is expressed by the following equation.
In the subsequent processing, the one that minimizes the overlap between the particle i and the particle j is found while slightly changing the direction of the vector n2. The initial value δmax (= about 0.1) of the minute amount δ used in the processing is set in step S1.
次に、２粒子の重なり量Ｌ（＝Ｄ２ｍａｘ−Ｄ２ｍｉｎ）を算出する（ステップＳ２）。
次に、接触の可能性があるか否かを判定する（ステップＳ３）。Ｌが０以下の場合、２粒子は非接触であるので、接触判定から除外される（ステップＳ４）。Ｌが０より大きい場合は接触の可能性があるため、重なり最小ベクトルｎ２を更新（探索）するためのループに入る。
Next, the overlapping amount L (= D2max−D2min) of the two particles is calculated (step S2).
Next, it is determined whether there is a possibility of contact (step S3). When L is 0 or less, the two particles are non-contact and are excluded from the contact determination (step S4). If L is greater than 0, there is a possibility of contact, and therefore a loop for updating (searching) the overlapping minimum vector n2 is entered.
重なり最小ベクトルｎ２に直交し、たがいに直交するベクトルｕ、ｖを適当に定める（ステップＳ５）。そして、試行する重なり最小ベクトルｎ２（１）〜ｎ２（４）を下記式のように設定する（ステップＳ６）。
上記式の１行目で決定される試行重なり最小ベクトルｎ２（１）は、図９中に示されるように、直交ベクトルｕの方向へ微小量δで決まる量だけ回転させたものとなる。
Vectors u and v that are orthogonal to the overlapping minimum vector n2 and are orthogonal to each other are appropriately determined (step S5). Then, the minimum overlap vector n2 (1) to n2 ( 4) to be tried is set as shown in the following equation (step S6).
As shown in FIG. 9, the trial overlap minimum vector n2 (1) determined in the first row of the above equation is rotated in the direction of the orthogonal vector u by an amount determined by a minute amount δ.
このようにして決定された各試行重なり最小ベクトルｎ２（ｋ）に対して２粒子の重なり量Ｌｎｅｗ（ｋ）を新たに算出する（ｋ＝１，２，３，４）。算出された４つのＬｎｅｗ（ｋ）の最小値とＬを比較する（ステップＳ７）。ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより小さければ、重なり最小ベクトルｎ２をｎ２（ｋ）に更新する。また、ＬもＬｎｅｗ（ｋ）に更新する（ステップＳ８）。重なり量Ｌを最小とするループにおいて、ｍｉｎ｛Ｌｎｅｗ（ｋ）｝がＬより大きくなった場合は（ステップＳ９）、微小量δをδ／２と小さくした後、試行重なりベクトルの設定から繰り返す（ステップＳ１０）。 An overlap amount Lnew (k) of two particles is newly calculated for each trial overlap minimum vector n2 (k) determined in this way (k = 1, 2, 3, 4). The calculated minimum value of four Lnew (k) is compared with L (step S7). If min {Lnew (k) } is smaller than L, the overlapping minimum vector n2 is updated to n2 (k) . Also, L is updated to Lnew (k) (step S8). When min {Lnew (k) } is larger than L in the loop that minimizes the overlap amount L (step S9), the minute amount δ is reduced to δ / 2, and the process is repeated from the setting of the trial overlap vector ( Step S10).
上記のループを繰り返して、δがδｍｉｎ（＝０．００１程度）より小さくなったときのｎ２が最終的な重なり最小ベクトルとなる（ステップＳ１１）。このような処理により、粒子ｉとｊの重なり量Ｌを最も小さくするベクトルｎ２を得ることができる。 By repeating the above loop, n2 when δ becomes smaller than δmin (= 0.001) becomes the final overlap minimum vector (step S11). By such processing, a vector n2 that minimizes the overlapping amount L of the particles i and j can be obtained.
なお、図９のステップＳ５で示したように、ループを繰り返すたびにｕ，ｖを面内で回転させると、局所解への陥りを抑制し、重なり最小ベクトルの最適解を得やすいという効果が得られる。 As shown in step S5 of FIG. 9, if u and v are rotated in the plane every time the loop is repeated, the effect of suppressing the fall into the local solution and easily obtaining the optimum solution of the minimum overlap vector is obtained. can get.
図９で示した方法は重なり最小ベクトルｎ２を決定するための一手法にすぎない。他の公知の最適化手法を利用して、重なり量Ｌを最小とするベクトルｎ２を求めることも可能である。また、少なくともベクトルｎ１とｎ２の向きが異なっていれば、第１の抽出方法よりも判定対象を絞り込む効果は得られるため、ベクトルｎ１に対して任意の回転を行ったものをベクトルｎ２として用いるだけでもよい。 The method shown in FIG. 9 is only one method for determining the minimum overlap vector n2. It is also possible to obtain the vector n2 that minimizes the overlap amount L using another known optimization method. If at least the directions of the vectors n1 and n2 are different, an effect of narrowing down the determination target can be obtained as compared with the first extraction method. Therefore, a vector n1 obtained by arbitrarily rotating the vector n1 is used. But you can.
（接触判定処理）
接触判定部Ｂ１２０の処理フローを図１０で説明する。
まず、接触判定を行う２粒子ｉ、ｊを設定する（ステップＳ１００）。次に２粒子ｉ、ｊの中心座標から中心結線単位ベクトルｎ１を求める（ステップＳ１１０）。ｎ１を基に２粒子ｉ、ｊの重なり量Ｌを算出する（ステップＳ１２０）。重なり量Ｌを基に接触可能性を判定し（ステップＳ１３０）、Ｌ＞０でなければ、２粒子は非接触として接触判定を終了する（ステップＳ１３５）。次に、２粒子ｉ、ｊの位置座標情報から重なり最小単位ベクトルｎ２を求める（ステップＳ１４０）。ステップＳ１４０の処理内容は、図９で詳しく説明したとおりである。
(Contact judgment process)
A processing flow of the contact determination unit B120 will be described with reference to FIG.
First, two particles i and j for which contact determination is performed are set (step S100). Next, a center connection unit vector n1 is obtained from the center coordinates of the two particles i and j (step S110). The overlap amount L of the two particles i and j is calculated based on n1 (step S120). The contact possibility is determined based on the overlap amount L (step S130), and if L> 0, the two particles are not in contact and the contact determination is terminated (step S135). Next, the overlapping minimum unit vector n2 is obtained from the position coordinate information of the two particles i and j (step S140). The processing content of step S140 is as described in detail with reference to FIG.
ここで、ステップＳ１１０で設定したｎ１を用いて接触可能頂点集合Ａを設定し（ステップＳ１５０）、ステップＳ１４０で設定したｎ２を用いて接触可能頂点集合Ｂを設定する（ステップＳ１６０）。次に、接触可能頂点集合Ａ∩Ｂを設定する（ステップＳ１７０）。こうして設定した接触可能頂点を含む粒子ｉの辺と粒子ｊの面を抽出する（ステップＳ１８０）。以上で、接触判定に用いる判定対象が絞り込まれる。 Here, the contactable vertex set A is set using n1 set in step S110 (step S150), and the contactable vertex set B is set using n2 set in step S140 (step S160). Next, a contactable vertex set A∩B is set (step S170). The side of the particle i and the surface of the particle j including the contactable vertex set in this way are extracted (step S180). The determination target used for contact determination is narrowed down as described above.
次に、粒子ｉが粒子ｊに接触しているか否かの判定処理を行う。まず、接触判定を行う粒子ｉの辺を設定する（ステップＳ１９０）。接触判定を行う粒子ｊの面を設定する（ステップＳ２００）。ここで設定した粒子ｉの辺と粒子ｊの面との交差の有無を判定する（ステップＳ２１０）。ここで、粒子ｊの面に粒子ｉの辺が交差していれば、交差回数のカウントアップを行う（ステップＳ２２０）。この操作をステップＳ１８０で抽出された粒子ｊの面の全てに対して行う（ステップＳ２３０）。次に、粒子ｉの辺と粒子ｊの面の交差回数が２回かどうかの判定を行う（ステップＳ２４０）。交差回数が２回ならば、めり込み量を算出し、メモリに記憶しておく（ステップＳ２５０）。本操作をＳ１８０で抽出された粒子ｉの辺の全てに対して行う（ステップＳ２６０）。 Next, it is determined whether or not the particle i is in contact with the particle j. First, the side of the particle i for which contact determination is performed is set (step S190). The surface of the particle j for which contact determination is performed is set (step S200). The presence / absence of intersection between the side of the particle i set here and the surface of the particle j is determined (step S210). Here, if the side of the particle i intersects the surface of the particle j, the number of intersections is counted up (step S220). This operation is performed on all the surfaces of the particle j extracted in step S180 (step S230). Next, it is determined whether the number of intersections between the side of the particle i and the surface of the particle j is 2 (step S240). If the number of intersections is two, the amount of sinking is calculated and stored in the memory (step S250). This operation is performed on all the sides of the particle i extracted in S180 (step S260).
なお、図１０では、説明を簡単にするため、粒子ｉのある辺が粒子ｊの２つの面と交差した場合に、接触有りと判定し、粒子ｉとｊのあいだのめり込み量を算出する処理例を示した。しかし実際の処理では、図４に示したように、粒子ｉと粒子ｊの接触を４つのパターンで判定する。前述したようにパターンごとに判定条件が異なるが、いずれのパターンにおいても粒子ｉの辺と粒子ｊの面との交差判定が基本となるので、図１０のステップＳ１８０で抽出した辺及び面を利用すればよい。 In FIG. 10, for the sake of simplicity of explanation, when a certain side of the particle i intersects with two surfaces of the particle j, it is determined that there is a contact, and a processing example for calculating the amount of subtraction between the particles i and j showed that. However, in the actual processing, as shown in FIG. 4, the contact between the particle i and the particle j is determined by four patterns. As described above, although the determination conditions differ for each pattern, since the determination of the intersection between the side of the particle i and the surface of the particle j is fundamental in any pattern, the side and the surface extracted in step S180 in FIG. 10 are used. do it.
各パターン（ａ）〜（ｄ）におけるめり込み量の定義について、図４中の下段の図を用いて説明する。
（ａ）頂点のめり込み：粒子ｉの頂点ｉ−ｖ０と粒子ｊの面ｊ−ｆ０との距離をめり込み量とする。
（ｂ）辺のめり込み（ケース１）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｃ）辺のめり込み（ケース２）：粒子ｉの辺ｉ−ｅ０と粒子ｊの辺ｊ−ｅ０との距離をめり込み量とする。
（ｄ）辺のめり込み（ケース３）：粒子ｉの辺ｉ−ｅ０に対する、粒子ｊの辺ｊ−ｅ０、ｊ−ｅ１との距離をそれぞれの辺とのめり込み量とする。
The definition of the sinking amount in each of the patterns (a) to (d) will be described with reference to the lower diagram in FIG.
(A) Vertex dip: The distance between the vertex i-v0 of the particle i and the surface j-f0 of the particle j is defined as the dip amount.
(B) Side indentation (case 1): The distance between the side i-e0 of the particle i and the side j-e0 of the particle j is defined as the amount of indentation.
(C) Side indentation (case 2): The distance between the side i-e0 of the particle i and the side j-e0 of the particle j is defined as the amount of indentation.
(D) Side indentation (Case 3): The distance between the side i of the particle i and the sides j-e0 and j-e1 of the particle j is defined as the amount of indentation with each side.
（粒子受力計算部Ｂ１３０）
以上の４パターンにおいて、めり込み量（距離）の計算に用いた各粒子表面の２点を結ぶ方向を接触面法線方向と定める。この接触面法線方向は、接触点の（粒子ｊの面又は辺に対する）めり込み方向ということもできる。そして、接触判定部Ｂ１２０で求めためり込み量に基づいて、ヘルツの式を用いて、粒子ｉに働く接触面法線方向の接触力（バネ力）を求める。また各接触点に対して、接触が始まってからの接触面（上記接触面法線方向に垂直な面）内方向の移動量を求め、ミンドリンの式を用いて粒子ｉに働く接触面内方向の接触力（バネ力）を求める。また接触面法線方向の運動速度（めり込み速度）と面内方向の各運動速度（ずり速度）から、各方向に作用する粘性力（ダンパー力）を求める。そしてそれら接触力と粘性力の和に、粒子ｉに働くその他の力（重力、付着力、電磁力、空気の粘性抵抗等）を加えて、各粒子ｉが受ける並進力と回転力を求める。なお、本めり込み量から接触力を求める過程は、上記非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
(Particle force calculation unit B130)
In the above four patterns, the direction connecting two points on the surface of each particle used for calculating the amount of penetration (distance) is defined as the normal direction of the contact surface. This normal direction of the contact surface can also be referred to as a dip direction (with respect to the surface or side of the particle j) of the contact point. Then, the contact force in the normal direction of the contact surface (spring force) acting on the particle i is obtained using the Hertz formula based on the amount of stagnation obtained by the contact determination unit B120. In addition, for each contact point, the amount of movement in the contact surface (surface perpendicular to the contact surface normal direction) from the start of contact is obtained, and the contact surface inward direction acting on the particle i using the Mindlin equation The contact force (spring force) is obtained. Also, the viscous force (damper force) acting in each direction is obtained from the movement speed (sinking speed) in the normal direction of the contact surface and each movement speed (shear speed) in the in-plane direction. Then, other forces (gravity, adhesion force, electromagnetic force, air viscosity resistance, etc.) acting on the particle i are added to the sum of the contact force and the viscous force to obtain the translational force and the rotational force that each particle i receives. In addition, since the process of calculating | requiring contact force from this amount of indentation is the same as the method described in the said
（粒子運動計算部Ｂ１４０）
次に、粒子運動計算部Ｂ１４０における各粒子の刻み時間後の速度と位置の算出方法を説明する。この処理は、粒子受力計算部Ｂ１３０で求めた粒子に働く力による、各粒子の刻み時間後の速度と位置を求めるものである。すなわち、粒子運動計算部Ｂ１４０は、並進力からニュートンの運動方程式を用いて並進方向の加速度を求め、更に速度と位置を求める。また粒子運動計算部Ｂ１４０は、回転力からオイラーの運動方程式を用いて角加速度を求め、更に回転速度と回転位置を求める。本処理は、非特許文献１に記載されている手法と同じであるため、詳しい説明は割愛する。
(Particle motion calculator B140)
Next, a method of calculating the speed and position of each particle after the ticking time in the particle motion calculation unit B140 will be described. This process is to determine the speed and position of each particle after the ticking time due to the force acting on the particle obtained by the particle receiving force calculation unit B130. That is, the particle motion calculation unit B140 calculates the acceleration in the translation direction from the translation force using Newton's equation of motion, and further determines the velocity and position. Further, the particle motion calculation unit B140 calculates angular acceleration from the rotational force using Euler's equation of motion, and further determines the rotational speed and rotational position. Since this processing is the same as the method described in
（実施例）
以上の接触判定方法を用いて粒子挙動計算を実施した例を、図１１〜図１２を用いて説明する。図１１（ａ）は解析モデルを示す。解析モデルは電子写真技術を用いた画像形成装置におけるクリーニングプロセスである。０．１ｍ／ｓｅｃで矢印の方向に移動する感光体ドラムに対してクリーニングブレードをβ＝３０°で当接させておき、感光体ドラム上に乗って運ばれてきたトナーの運動を見る。各トナーの形状を粒子形状データとして与え、感光体ドラム及びクリーニングブレードの形状を壁形状データとして与える。
(Example)
An example in which the particle behavior calculation is performed using the above contact determination method will be described with reference to FIGS. FIG. 11A shows an analysis model. The analysis model is a cleaning process in an image forming apparatus using an electrophotographic technique. A cleaning blade is brought into contact with the photosensitive drum moving in the direction of the arrow at 0.1 m / sec at β = 30 °, and the movement of the toner carried on the photosensitive drum is observed. The shape of each toner is given as particle shape data, and the shapes of the photosensitive drum and the cleaning blade are given as wall shape data.
図１１（ｂ）は、解析条件として与えた、トナー、感光体ドラム、クリーニングブレードそれぞれの機械特性（ヤング率、ポアソン比、減衰係数、摩擦係数、比重）を示している。この条件で、刻み時間Δｔを１０ｎｓｅｃに設定し、約２百万ステップの計算を実行した。 FIG. 11B shows the mechanical properties (Young's modulus, Poisson's ratio, damping coefficient, friction coefficient, specific gravity) of the toner, the photosensitive drum, and the cleaning blade given as analysis conditions. Under these conditions, the step time Δt was set to 10 nsec, and a calculation of about 2 million steps was executed.
図１２は、（ａ）球形粒子と（ｂ）異形粒子の計算結果を比較したものである。球形粒子の計算は従来の計算方法で実施している。異形粒子は４面体の角をとった形状であり、その表面を８つの３角形と６つの４角形の合計１２の多角形で構成している。なお本計算においては、トナーがブレードに堰き止められた際の挙動をみるため、強制的にトナーがすり抜けないようにしている。球形粒子に対して異形粒子は回転しにくいことから、全体的に動き辛く、トナー粒子の密集度が上がる傾向にある。図１２の計算結果では、異形粒子の密集度が高くなっている様子が再現されている。 FIG. 12 compares the calculation results of (a) spherical particles and (b) deformed particles. Spherical particles are calculated by a conventional calculation method. The irregularly shaped particle has a tetrahedral corner shape, and its surface is composed of a total of 12 polygons including 8 triangles and 6 rectangles. In this calculation, the toner is forcedly prevented from slipping through in order to observe the behavior when the toner is blocked by the blade. Since the irregularly shaped particles are difficult to rotate with respect to the spherical particles, it is difficult to move as a whole, and the density of the toner particles tends to increase. In the calculation result of FIG. 12, a state in which the density of irregularly shaped particles is high is reproduced.
このトナーがブレードに堰き止められた際の挙動は、クリーニング不良であるトナーのすり抜けとの関係があり、流動的である程すり抜けやすいといわれている。このことから、製品設計を行う上で、トナーの形状によるクリーニング不良の傾向を粒子挙動計算で予測することが可能となる。 The behavior of the toner when it is dammed by the blade is related to the slipping of the toner that is poorly cleaned, and it is said that the fluid is more easily slipped. This makes it possible to predict the tendency of poor cleaning due to the shape of the toner by particle behavior calculation when designing the product.
また、図１３は面数の異なる異形粒子の自由落下の計算を示した図である。（ａ）は、
頂点数１２、面数１４、辺数２４の場合の異形粒子で、（ｂ）は頂点数９８、面数１０８、辺数２０４の異形粒子の場合である。（ａ）と（ｂ）の場合において、２５６個の異形粒子の自由落下の計算を刻み時間Δｔ＝１０ｎｓｅｃで８百万ステップ実行し、本発明を使用しない場合と使用した場合での計算時間の比較を行った。（ａ）では、従来の方法では計算に１９時間３５分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込み（第１の抽出方法）を用いた場合では、８時間０９分となり、高速化率は２．４倍となる。（ｂ）では、従来の方法では計算に５２４時間３２分かかったのに対し、中心結線ベクトルによる接触可能頂点の絞り込みを用いると２２時間２４分となり、高速化率は２３．４倍となった。このことから、粒子の面の数が多くなればなるほど、本発明を用いることによる計算時間の短縮効果が大きくなることがわかる。
FIG. 13 is a diagram showing calculation of free fall of deformed particles having different numbers of faces. (A)
The figure shows a deformed particle having 12 vertices, 14 faces, and 24 sides, and (b) shows a deformed particle having 98 vertices, 108 faces, and 204 edges. In the cases of (a) and (b), the calculation of the free fall of 256 irregularly shaped particles was executed in 8 million steps at a ticking time Δt = 10 nsec, and the calculation time of the case where the present invention is not used and the case where it is used is calculated. A comparison was made. In (a), in the conventional method, the calculation took 19 hours and 35 minutes, whereas in the case of using the center connection vector to narrow down the contactable vertices (first extraction method), the calculation time was 8 hours and 09 minutes. The conversion rate is 2.4 times. In (b), the calculation with the conventional method took 524 hours and 32 minutes, but using the center connection vector to narrow the vertices that can be contacted gave 22 hours and 24 minutes, and the speed-up rate was 23.4 times. . From this, it can be seen that as the number of particle surfaces increases, the effect of shortening the calculation time by using the present invention increases.
図１４は中心結線ベクトルと重なり最小ベクトルを用いた高速化の効果を示す図である。計算には頂点数１９８、面数３９２、辺数５８８の異形粒子を用い、図１３と同様の自由落下計算を３０００００ステップ、Δｔ＝１０ｎｓｅｃで行った。接触判定に中心結線ベクトルのみを用いた場合（第１の抽出方法）と中心結線ベクトルと重なり最小ベクトルの両方を用いた場合（第２の抽出方法）で計算時間を比較した。中心結線ベクトルのみを用いた場合は４７２６秒、中心結線ベクトルと重なり最小ベクトルの両方を用いた場合は２００９秒となり、高速化率は２．３５倍となった。このことから、中心結線ベクトルと重なり最小ベクトルの両方を用いた接触判定により、さらなる計算時間の短縮が実現できていることが分かる。 FIG. 14 is a diagram showing the effect of speeding up using the center connection vector and the overlapping minimum vector. For the calculation, irregular-shaped particles having a number of vertices of 198, a number of surfaces of 392, and a number of sides of 588 were used, and free fall calculation similar to FIG. Computation time was compared between the case where only the center connection vector was used for the contact determination (first extraction method) and the case where both the center connection vector and the minimum overlap vector were used (second extraction method). When only the central connection vector was used, it was 4726 seconds, and when both the central connection vector and the overlapped minimum vector were used, it was 2009 seconds, and the speed-up rate was 2.35 times. From this, it is understood that the calculation time can be further reduced by the contact determination using both the center connection vector and the minimum overlap vector.
１０：情報処理装置
Ｂ１００：全体制御部
Ｂ１１０：入力データ読込部
Ｂ１２０：接触判定部
Ｂ１３０：粒子受力計算部
Ｂ１４０：粒子運動計算部
Ｂ１５０：計算結果出力部
10: Information processing device B100: Overall control unit B110: Input data reading unit B120: Contact determination unit B130: Particle force calculation unit B140: Particle motion calculation unit B150: Calculation result output unit
Claims (8)
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面を設定し、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。 A method of calculating the motion of a plurality of particles in a dense state by an information processing device,
An information processing device obtains data defining a shape of each of a plurality of particles by a polyhedron, and
A contact determination step in which the information processing apparatus determines whether each particle is in contact with another object based on the position information of each of the plurality of particles;
An information processing device calculates a force received by each particle from another object determined to be in contact in the contact determination step, and calculates a motion of each particle by the force;
An information processing apparatus including an output step of outputting the motion of each particle calculated in the calculation step;
In the contact determination step, when determining whether or not the particle i is in contact with the particle j,
Setting a first plane that intersects the straight line connecting the center of the particle i and the center of the particle j and touches the particle j ;
Determining a unit vector n1 perpendicular to the first plane and pointing toward the center of the particle j;
Calculate the inner product n1 · Pi of the unit vector n1 and the position vector Pi of each vertex of the particle i,
The inner product n1 · Pj of the unit vector n1 and the position vector Pj of each vertex of the particle j is calculated,
Of the inner products n1 · Pi for all the vertices of the particle i, the largest one is D1max, and among the inner products n1 · Pj for all the vertices of the particle j, the smallest one is D1min, the inner product n1 Define a set A consisting of the vertices of particles i where Pi is D1min to D1max and the vertices of particles j whose inner product n1 · Pj is D1min to D1max,
Extract the sides of the particle i and the surface of the particle j including the vertices included in the set A,
A particle behavior analysis method, wherein contact determination is performed between the side of the extracted particle i and the surface of the particle j.
情報処理装置が、複数の粒子それぞれの形状を多面体により定義するデータを取得する取得ステップと、
情報処理装置が、前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に
接触しているか否かを判定する接触判定ステップと、
情報処理装置が、前記接触判定ステップで接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算ステップと、
情報処理装置が、前記計算ステップで計算された各粒子の運動を出力する出力ステップと、を含み、
前記接触判定ステップにおいて、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析方法。 A method of calculating the motion of a plurality of particles in a dense state by an information processing device,
An information processing device obtains data defining a shape of each of a plurality of particles by a polyhedron, and
A contact determination step in which the information processing apparatus determines whether each particle is in contact with another object based on the position information of each of the plurality of particles;
An information processing device calculates a force received by each particle from another object determined to be in contact in the contact determination step, and calculates a motion of each particle by the force;
An information processing apparatus including an output step of outputting the motion of each particle calculated in the calculation step;
In the contact determination step, when determining whether or not the particle i is in contact with the particle j,
Intersects a straight line connecting the center of the particle i and the center of the particle j and is closer to the center of the particle j than the first plane in contact with the particle j, and intersects the straight line at an angle different from the first plane. And the side of the particle i located in the space on the center side of the particle j from the third plane in contact with the particle j, and
Particles closer to the center side of the particle i than the second plane parallel to the first plane and in contact with the particle i, and parallel to the third plane and in contact with the particle i than the fourth plane extract the surface of the particle j located in the center space of i,
A particle behavior analysis method, wherein contact determination is performed between the side of the extracted particle i and the surface of the particle j.
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
前記第３の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ２を決定し、
単位ベクトルｎ２と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ２・Ｐｉを計算し、
単位ベクトルｎ２と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ２・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ２・Ｐｉのうち、最大となるものをＤ２ｍａｘ、粒子ｊの全ての頂点についての内積ｎ２・Ｐｊのうち、最小となるものをＤ２ｍｉｎとした場合に、内積ｎ２・ＰｉがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ２・ＰｊがＤ２ｍｉｎ以上Ｄ２ｍａｘ以下となる粒子ｊの頂点からなる集合Ｂを定義し、
集合Ａと集合Ｂの積集合に含まれる頂点を含む粒子ｉの辺及び粒子ｊの面と、集合Ａに含まれる頂点と集合Ｂに含まれる頂点の両方を含む粒子ｉの辺及び粒子ｊの面とを抽出するステップである
ことを特徴とする請求項２に記載の粒子挙動解析方法。 Extracting the sides of the particle i and the surface of the particle j
Determining a unit vector n1 perpendicular to the first plane and pointing toward the center of the particle j;
Calculate the inner product n1 · Pi of the unit vector n1 and the position vector Pi of each vertex of the particle i,
The inner product n1 · Pj of the unit vector n1 and the position vector Pj of each vertex of the particle j is calculated,
Of the inner products n1 · Pi for all the vertices of the particle i, the largest one is D1max, and among the inner products n1 · Pj for all the vertices of the particle j, the smallest one is D1min, the inner product n1 Define a set A consisting of the vertices of particles i where Pi is D1min to D1max and the vertices of particles j whose inner product n1 · Pj is D1min to D1max,
Determining a unit vector n2 perpendicular to the third plane and facing the center of the particle j;
Calculate the inner product n2 · Pi of the unit vector n2 and the position vector Pi of each vertex of the particle i,
The inner product n2 · Pj of the unit vector n2 and the position vector Pj of each vertex of the particle j is calculated,
Of the inner products n2 · Pi for all the vertices of the particle i, the largest one is D2max, and among the inner products n2 · Pj for all the vertices of the particle j, the smallest one is D2min, the inner product n2 Define a set B consisting of the vertices of particles i where Pi is D2min or more and D2max or less, and the vertices of particles j whose inner product n2 · Pj is D2min or more and D2max or less,
The side of particle i and the surface of particle j including the vertex included in the intersection set of set A and set B, the side of particle i including both the vertex included in set A and the vertex included in set B, and the particle j The particle behavior analysis method according to claim 2 , wherein the method is a step of extracting a surface.
ことを特徴とする請求項１又は３に記載の粒子挙動解析方法。 The particle behavior analysis method according to claim 1 or 3 , wherein the unit vector n1 is determined so as to be parallel to a straight line connecting the center of the particle i and the center of the particle j.
前記単位ベクトルｎ２を、Ｄ２ｍａｘ−Ｄ２ｍｉｎがＤ１ｍａｘ−Ｄ１ｍｉｎよりも小さくなるように決定する
ことを特徴とする請求項３に記載の粒子挙動解析方法。 The unit vector n1 is determined to be parallel to a straight line connecting the center of the particle i and the center of the particle j;
The particle behavior analysis method according to claim 3 , wherein the unit vector n2 is determined so that D2max-D2min is smaller than D1max-D1min.
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面を設定し、
前記第１の平面に垂直で且つ粒子ｊの中心側を向く単位ベクトルｎ１を決定し、
単位ベクトルｎ１と粒子ｉの各頂点の位置ベクトルＰｉの内積ｎ１・Ｐｉを計算し、
単位ベクトルｎ１と粒子ｊの各頂点の位置ベクトルＰｊの内積ｎ１・Ｐｊを計算し、
粒子ｉの全ての頂点についての内積ｎ１・Ｐｉのうち、最大となるものをＤ１ｍａｘ、粒子ｊの全ての頂点についての内積ｎ１・Ｐｊのうち、最小となるものをＤ１ｍｉｎとした場合に、内積ｎ１・ＰｉがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｉの頂点、及び、内積ｎ１・ＰｊがＤ１ｍｉｎ以上Ｄ１ｍａｘ以下となる粒子ｊの頂点からなる集合Ａを定義し、
集合Ａに含まれる頂点を含む粒子ｉの辺及び粒子ｊの面を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。 A particle behavior analyzer for calculating the motion of a plurality of particles in a dense state,
An acquisition unit for acquiring data defining the shape of each of a plurality of particles by a polyhedron;
A contact determination unit that determines whether or not each particle is in contact with another object based on position information of each of the plurality of particles;
Calculating a force received by each particle from another object determined to be in contact by the contact determination unit, and calculating a motion of each particle by the force;
A calculation result output unit for outputting the motion of each particle calculated by the calculation unit,
The contact determination unit determines whether or not the particle i is in contact with the particle j.
Setting a first plane that intersects the straight line connecting the center of the particle i and the center of the particle j and touches the particle j ;
Determining a unit vector n1 perpendicular to the first plane and pointing toward the center of the particle j;
Calculate the inner product n1 · Pi of the unit vector n1 and the position vector Pi of each vertex of the particle i,
The inner product n1 · Pj of the unit vector n1 and the position vector Pj of each vertex of the particle j is calculated,
Of the inner products n1 · Pi for all the vertices of the particle i, the largest one is D1max, and among the inner products n1 · Pj for all the vertices of the particle j, the smallest one is D1min, the inner product n1 Define a set A consisting of the vertices of particles i where Pi is D1min to D1max and the vertices of particles j whose inner product n1 · Pj is D1min to D1max,
Extract the sides of the particle i and the surface of the particle j including the vertices included in the set A,
A particle behavior analyzing apparatus characterized in that contact determination is performed between the side of the extracted particle i and the surface of the particle j.
複数の粒子それぞれの形状を多面体により定義するデータを取得する取得部と、
前記複数の粒子それぞれの位置情報に基づいて、各粒子が他の物体に接触しているか否かを判定する接触判定部と、
前記接触判定部で接触していると判定された他の物体から各粒子が受ける力を計算し、その力による各粒子の運動を計算する計算部と、
前記計算部で計算された各粒子の運動を出力する計算結果出力部と、を含み、
前記接触判定部は、粒子ｉが粒子ｊに接触しているか否かを判定する場合に、
粒子ｉの中心と粒子ｊの中心を結ぶ直線に交わり且つ粒子ｊに接する第１の平面よりも粒子ｊの中心側であって、且つ、前記第１の平面とは異なる角度で前記直線に交わり且つ粒子ｊに接する第３の平面よりも粒子ｊの中心側の空間に位置する粒子ｉの辺、及び、
前記第１の平面に平行で且つ粒子ｉに接する第２の平面よりも粒子ｉの中心側であって、且つ、前記第３の平面に平行で且つ粒子ｉに接する第４の平面よりも粒子ｉの中心側の空間に位置する粒子ｊの面、を抽出し、
前記抽出した粒子ｉの辺と粒子ｊの面とのあいだで接触判定を行う
ことを特徴とする粒子挙動解析装置。 A particle behavior analyzer for calculating the motion of a plurality of particles in a dense state,
An acquisition unit for acquiring data defining the shape of each of a plurality of particles by a polyhedron;
A contact determination unit that determines whether or not each particle is in contact with another object based on position information of each of the plurality of particles;
Calculating a force received by each particle from another object determined to be in contact by the contact determination unit, and calculating a motion of each particle by the force;
A calculation result output unit for outputting the motion of each particle calculated by the calculation unit,
The contact determination unit determines whether or not the particle i is in contact with the particle j.
Intersects a straight line connecting the center of the particle i and the center of the particle j and is closer to the center of the particle j than the first plane in contact with the particle j, and intersects the straight line at an angle different from the first plane. And the side of the particle i located in the space on the center side of the particle j from the third plane in contact with the particle j, and
Particles closer to the center side of the particle i than the second plane parallel to the first plane and in contact with the particle i, and parallel to the third plane and in contact with the particle i than the fourth plane extract the surface of the particle j located in the center space of i,
A particle behavior analyzing apparatus characterized in that contact determination is performed between the side of the extracted particle i and the surface of the particle j.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2011230304A JP5893333B2 (en) | 2011-10-20 | 2011-10-20 | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2011230304A JP5893333B2 (en) | 2011-10-20 | 2011-10-20 | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2013089100A JP2013089100A (en) | 2013-05-13 |
| JP5893333B2 true JP5893333B2 (en) | 2016-03-23 |
Family
ID=48532930
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2011230304A Active JP5893333B2 (en) | 2011-10-20 | 2011-10-20 | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP5893333B2 (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5901417B2 (en) * | 2012-05-11 | 2016-04-13 | キヤノン株式会社 | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program |
| KR102090456B1 (en) * | 2012-11-30 | 2020-03-18 | 엘지디스플레이 주식회사 | non-conductive type adhesive means and display device including the same |
| JP6520447B2 (en) * | 2015-06-18 | 2019-05-29 | 住友ベークライト株式会社 | Simulation method, simulation apparatus, and computer program |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPS6426981A (en) * | 1987-07-23 | 1989-01-30 | Toshiba Corp | Recognizing method for collision of objects |
| JP3391657B2 (en) * | 1996-05-22 | 2003-03-31 | 富士通株式会社 | Inter-object distance calculation device and inter-object distance calculation program storage medium |
| US5761391A (en) * | 1996-05-22 | 1998-06-02 | Fujitsu Ltd. | Arithmetic unit for calculating distance between objects |
| JP5183698B2 (en) * | 2010-01-22 | 2013-04-17 | キヤノン株式会社 | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program |
-
2011
- 2011-10-20 JP JP2011230304A patent/JP5893333B2/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| JP2013089100A (en) | 2013-05-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhao et al. | A poly‐superellipsoid‐based approach on particle morphology for DEM modeling of granular media | |
| Lim et al. | Granular element method for three‐dimensional discrete element calculations | |
| Lu et al. | Critical assessment of two approaches for evaluating contacts between super-quadric shaped particles in DEM simulations | |
| US9311745B2 (en) | Systems and methods of analysis of granular elements | |
| Nezamabadi et al. | Implicit frictional-contact model for soft particle systems | |
| Fu et al. | Polyarc discrete element for efficiently simulating arbitrarily shaped 2D particles | |
| Jin et al. | Probability-based contact algorithm for non-spherical particles in DEM | |
| Valera et al. | Modified algorithm for generating high volume fraction sphere packings | |
| Neto et al. | Discrete element model for general polyhedra | |
| Tang et al. | Energy conservative property of impulse‐based methods for collision resolution | |
| JP6671064B2 (en) | Particle simulation apparatus, particle simulation method, and particle simulation program | |
| JP5183698B2 (en) | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program | |
| Jansson et al. | A discrete mechanics model for deformable bodies | |
| JP5704246B2 (en) | Object motion analysis apparatus, object motion analysis method, and object motion analysis program | |
| He et al. | Simulations of realistic granular soils in oedometer tests using physics engine | |
| JP5893333B2 (en) | Particle behavior analysis method, particle behavior analysis apparatus, and analysis program | |
| Tschisgale et al. | A constraint-based collision model for Cosserat rods | |
| Zhang et al. | Metaball based discrete element method for general shaped particles with round features | |
| Zhang et al. | Investigation of particle shape and ambient fluid on sandpiles using a coupled micro-geomechanical model | |
| Lai et al. | Machine‐learning‐enabled discrete element method: Contact detection and resolution of irregular‐shaped particles | |
| Zhang et al. | Coupled metaball discrete element lattice Boltzmann method for fluid-particle systems with non-spherical particle shapes: A sharp interface coupling scheme | |
| He et al. | Simulating shearing behavior of realistic granular soils using physics engine | |
| US9020784B2 (en) | Methods for providing a bonded-particle model in computer aided engineering system | |
| He et al. | Physics engine based simulation of shear behavior of granular soils using hard and soft contact models | |
| JP5454693B2 (en) | Continuum motion analysis program, continuum motion analysis method, and continuum motion analysis apparatus |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20141017 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20150731 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20150804 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20151005 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20160126 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20160224 |
|
| R151 | Written notification of patent or utility model registration |
Ref document number: 5893333 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R151 |