WO2013038476A1 - 運動解析装置、運動解析方法及び運動解析プログラム - Google Patents

運動解析装置、運動解析方法及び運動解析プログラム Download PDF

Info

Publication number
WO2013038476A1
WO2013038476A1 PCT/JP2011/070764 JP2011070764W WO2013038476A1 WO 2013038476 A1 WO2013038476 A1 WO 2013038476A1 JP 2011070764 W JP2011070764 W JP 2011070764W WO 2013038476 A1 WO2013038476 A1 WO 2013038476A1
Authority
WO
WIPO (PCT)
Prior art keywords
particles
motion analysis
particle
motion
divided
Prior art date
Application number
PCT/JP2011/070764
Other languages
English (en)
French (fr)
Inventor
圭太 小笠原
正喜 風間
Original Assignee
富士通株式会社
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 富士通株式会社 filed Critical 富士通株式会社
Priority to JP2013533361A priority Critical patent/JP5761355B2/ja
Priority to PCT/JP2011/070764 priority patent/WO2013038476A1/ja
Publication of WO2013038476A1 publication Critical patent/WO2013038476A1/ja
Priority to US14/203,701 priority patent/US20140195212A1/en

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to a motion analysis device, a motion analysis method, and a motion analysis program.
  • the SPH (Smoothed Particle Hydrodynamics) method, the MPS (Moving Particle Semi-implicit) method, and the like are known (for example, see Non-Patent Documents 1 and 2).
  • SPH Smoothed Particle Hydrodynamics
  • MPS Moving Particle Semi-implicit
  • particles existing within a distance h called an influence radius from a certain particle are treated as neighboring particles, and the motion of the particles is analyzed assuming that a repulsive force acts between particles within the influence radius h.
  • Equation (1) is an example of an equation obtained by discretizing the equation of motion in the SPH method.
  • subscripts a and b are identifiers for identifying particles.
  • “#” Indicates that the following alphabet is a vector.
  • W is a kernel function.
  • the kernel function W is used to construct a continuous field from the distribution of particles.
  • a cubic spline function represented by the following equation (2) is used as the kernel function W.
  • h is an influence radius between particles, and about 2 to 3 times the average particle interval in the initial state is often used.
  • ⁇ is a value adjusted so that the total spatial integration amount of the kernel function is 1, and is determined to be 0.7 ⁇ h 2 in the case of two dimensions and ⁇ h 3 in the case of three dimensions.
  • An invention about a method of performing a flow analysis using a particle method is disclosed (for example, see Patent Document 1).
  • particles with different particle sizes are handled for each computer, and when each particle size of the particles existing in the analysis region in charge exceeds the upper limit value, each computer treats the particles as a plurality of particles. Divide into
  • the present invention aims to improve the accuracy of analysis of the motion of an object.
  • one embodiment of the present invention provides: A storage unit for storing data representing a plurality of particles representing an object; When a predetermined uneven distribution relationship exists between other particles existing within a predetermined range centered on the target particle included in the plurality of particles, data representing a plurality of divided particles obtained by dividing the target particle is stored in the memory.
  • a division processing unit stored in the unit;
  • a motion analysis unit that analyzes the motion of the object according to the type of the object, using data representing the plurality of divided particles stored in the storage unit; It is a motion analysis device characterized by having.
  • FIG. 6 is a diagram illustrating a state in which a particle a is divided into six particles a_1 to a_6. It is a figure which shows a mode that the particle
  • FIG. 6 is a diagram illustrating a state in which a particle a is divided into six particles a_1 to a_6. It is a figure which shows a mode that the particle
  • FIG. 4 is a diagram illustrating a state in which a particle a is divided into four particles a_1 to a_4. It is a figure which shows notionally density (rho) * represented by the kernel function W.
  • 3 is a flowchart showing a flow of overall processing executed by the motion analysis apparatus 1. It is a flowchart which shows the flow of the process at the time of determining whether a division
  • the motion analysis device 1 is a device that analyzes the motion of an object such as a fluid, an elastic body, or a rigid body in a virtual three-dimensional space or a two-dimensional plane based on a particle method.
  • FIG. 1 is a hardware configuration example of a motion analysis apparatus 1 according to the first embodiment of the present invention.
  • the motion analysis device 1 includes, for example, a CPU (Central Processing Unit) 10, a drive device 12, an auxiliary storage device 16, a memory device 18, an interface device 20, an input device 22, and a display device 24.
  • Information processing apparatus These components are connected via a bus, a serial line, or the like.
  • the CPU 10 is an arithmetic processing unit such as a processor having a program counter, an instruction decoder, various arithmetic units, an LSU (Load Store Unit), a general-purpose register, and the like.
  • arithmetic processing unit such as a processor having a program counter, an instruction decoder, various arithmetic units, an LSU (Load Store Unit), a general-purpose register, and the like.
  • the drive device 12 is a device that can read a program and data from the storage medium 14.
  • the program is installed from the storage medium 14 to the auxiliary storage device 16 via the drive device 12.
  • the storage medium 14 is a portable storage medium such as a CD-ROM, a DVD disk, or a USB memory.
  • the auxiliary storage device 16 is, for example, an HDD (Hard Disk Drive) or a flash memory.
  • the program can be installed by the interface device 20 being downloaded from another computer via a network and installed in the auxiliary storage device 16.
  • the network is the Internet, a LAN (Local Area Network), a wireless network, or the like.
  • the program may be stored in advance in the auxiliary storage device 16 or a ROM (Read Only Memory) when the information processing apparatus is shipped.
  • the information processing apparatus shown in FIG. 1 can function as the motion analysis apparatus 1 of the present embodiment.
  • the memory device 18 is, for example, a RAM (Random Access Memory) or an EEPROM (Electrically Erasable and Programmable Read Only Memory).
  • the interface device 20 controls connection with the network.
  • the input device 22 is, for example, a keyboard, a mouse, a touch pad, a touch panel, a microphone, or the like.
  • the display device 24 is a display device such as an LCD (Liquid Crystal Display) or a CRT (Cathode Ray Tube).
  • the motion analysis apparatus 1 may include other types of output devices such as a printer and a speaker in addition to the display device 24.
  • FIG. 2 is a functional configuration example of the motion analysis apparatus 1 according to an embodiment of the present invention.
  • the motion analysis apparatus 1 includes a data input unit 30, a division processing unit 32, and a motion analysis unit 34. These functional blocks function when the CPU 10 executes program software stored in the auxiliary storage device 16 or the like. Note that these functional blocks do not need to be realized by a clearly separated program, and may be called from other programs as subroutines or functions. Further, a part of the functional blocks may be hardware means such as an IC (Integrated Circuit) or an FPGA (Field Programmable Gate Array).
  • IC Integrated Circuit
  • FPGA Field Programmable Gate Array
  • the data input unit 30 receives input of data related to particles constituting the object and stores the data in the memory device 18.
  • the input data is the aggregate data of particles modeling fluids, elastic bodies, etc., and specifically includes particle center coordinates, velocity, influence radius, density, mass, deformation gradient tensor, material properties, temperature, etc. .
  • auxiliary storage device 16 may be stored in advance in the auxiliary storage device 16 or may be input by the user using the input device 22. These data may be installed from the storage medium 14 to the auxiliary storage device 16 or may be acquired from the network by the interface device 20.
  • the deformation gradient tensor F ij is generally expressed by the following equation (4).
  • X j is the j component of the distance between the positions of the particles
  • u i is the i component of the displacement of the particle
  • I is the identity tensor.
  • the i component and the j component represent three components of x, y, and z representing three-dimensional coordinates if they are three-dimensional.
  • the deformation gradient tensor of the particle a in the particle method is obtained by superimposing the influence from the particle b existing in the region where the distance from the particle a is within the influence radius of the particle a with the particle a as the center. expressed.
  • du represents the displacement of the particle a.
  • “present in the region where the distance from the particle a around the particle a is within the influence radius of the particle a” is simply referred to as “present around the particle a”.
  • the data input unit 30 may have a function of setting data using particles from object data in case the input data is not expressed using particles.
  • the division processing unit 32 considers each data such as the position of each particle and the radius of influence stored in the memory device 18, and performs the particle division determination and the division processing.
  • the particle division process is performed in a situation where a certain particle a does not accurately act on the repulsive force that should be received from particles existing around the particle a.
  • Predetermined uneven distribution relationship means, for example, (A) a situation in which the number of other particles b existing around the particle a is less than a predetermined number (for example, about three in three dimensions, about two in two dimensions), Or (B) The relationship in which all the other particles b existing around the particle a and the particle a exist on substantially the same plane or substantially the same straight line.
  • FIG. 3 is a diagram showing the situation (A) (two particles b), and FIGS. 4 and 5 are diagrams showing the situation (B).
  • f indicates the repulsive force that the particle a receives from the particle b.
  • the plane is represented as S, and the straight line is represented as L in FIG.
  • the situation where the particles a and other particles b exist on substantially the same plane is excluded from the “predetermined uneven distribution relationship”.
  • “exists on substantially the same plane” means, for example, “a plane including the particle a and any two particles b existing around the particle a is a reference plane S, and all the particles b are used as a reference. It can be defined that the distance D1 from the plane S is within a predetermined value Dref1. Further, the present invention is not limited to this, and “exists on substantially the same plane” is defined as “the angle ⁇ 1 formed by the vector from the particle a to all the particles b and the reference plane S is within a predetermined angle ⁇ ref1”. May be. 6 and 7 are diagrams illustrating the positional relationship between the particles a and b defined to exist on substantially the same plane.
  • “exists on substantially the same straight line” means, for example, “a straight line including the particle a and any one particle b existing around the particle a is set as a reference straight line L, and the reference straight line for all particles b. It can be defined that the distance D2 from L is within a predetermined value Dref2. Further, the present invention is not limited thereto, and “exists on substantially the same straight line” is defined as “the angle ⁇ 2 formed by the vector from the particle a to all particles b and the reference straight line L is within a predetermined angle ⁇ ref2”. May be. 8 and 9 are diagrams illustrating the positional relationship between the particles a and b that are defined to exist on substantially the same straight line.
  • the particles a do not receive a sufficient repulsive force from the surrounding particles b, and therefore can move relatively freely.
  • the particle a can move relatively freely in the direction orthogonal to the reference plane S or the reference straight line L. Since neither situation is appropriate as the actual behavior of the particles constituting the fluid or the elastic body, when the situation (A) or (B) above occurs for the particle a, the division processing unit 32 causes the particle a to be divided into a plurality of particles a. The data is divided into particles, and the data such as the position and mass of the new particles after the division are stored in the memory device 18.
  • the division processing unit 32 calculates a restandardized matrix L (#r) represented by the following equation (6), obtains a determinant detL (#r), and compares the determinant detL (#r) with the threshold Ref. By doing so, the particle division is determined. Specifically, the division processing unit 32 determines that the particle needs to be divided when the determinant detL (#r) exceeds the threshold Ref, and the determinant detL (#r) is equal to or less than the threshold Ref. In this case, it is determined that the particle division is unnecessary.
  • following Formula (6) has shown the restandardization matrix L (#r) in case the space which the motion analysis apparatus 1 handles is three-dimensional. #r is a position vector of the particle a, and (x, y, z) is a component.
  • the restandardized matrix L is a 2 ⁇ 2 matrix and is expressed by the following equation (7).
  • #r has (x, y) as a component.
  • detL (#r) is a value close to 1 if the particles b existing around the particles a are evenly distributed.
  • detL (#r) increases under the conditions (A) and (B) described above, and can take values up to infinity. Therefore, the division processing unit 32 performs the particle a division processing when the determinant detL (#r) is equal to or greater than the threshold value Ref.
  • the threshold value Ref is arbitrarily set within a range exceeding 1 and not exceeding 2 due to the property of the determinant detL (#r) as described above.
  • the threshold value Ref is a parameter that determines how much “substantially on the same plane” or “substantially on the same straight line” is included.
  • the larger the threshold value the lower the possibility of performing the division process. Therefore, “substantially on the same plane” and “substantially on the same straight line” include a narrower range.
  • the smaller the threshold value the higher the possibility of performing the division process. Therefore, “substantially on the same plane” and “substantially on the same straight line” include a wider range.
  • division processing unit 32 determines that the particle a division processing is necessary by the above-described method, the division processing unit 32 executes the particle a division processing according to the following procedure.
  • the number i for dividing the particle a is determined.
  • the number i may be arbitrarily determined, but it is preferable to previously determine a division rule according to the situation as follows.
  • FIG. 10 is a diagram illustrating a state in which the particle a is divided into six particles a_1 to a_6.
  • FIG. 11 is a diagram illustrating a state in which the particle a is divided into two particles a_1 and a_2.
  • FIG. 12 is a diagram illustrating a state in which the particle a is divided into four particles a_1 to a_4.
  • the uneven distribution direction that is, a direction or straight line parallel to the plane. It is preferable to divide the particles a in a direction perpendicular to the direction.
  • the influence radius of the particles after the division may be the same as the influence radius of the particles before the division.
  • each particle may not be divided independently, but two particles may be divided into a total of three particles.
  • the density distribution represented by the group of particles after the division is as close as possible to the density distribution of the particles before the division.
  • the density ⁇ * of the particle # x * before division is expressed by the following equation (8).
  • M is the mass of the particle # x * before division
  • W is the kernel function described above.
  • FIG. 13 is a diagram conceptually showing the density ⁇ * expressed by the kernel function W.
  • the density ⁇ represented by the divided particle group is represented by the following formula (9) by the superposition formula of the particle method.
  • Division processing unit 32 in order to find an optimal solution of the mass m i to approximate the density [rho density [rho * above constraint conditions, using Lagrange multipliers method (Lagrange Multiplier).
  • the division processing unit 32 defines energy J (m i ) represented by the following formula (10).
  • the division processing unit 32 a solution m i the defined energy J (m i) is the extremum following equation (11), the optimum mass of each particle after division.
  • division processing unit 32 the following equation (11) is any .delta.m i, and hold for [delta] [lambda], the above equation (10) is m i, for a formula of the secondary or less with respect to ⁇ From Equation (11), linear simultaneous equations relating to m i and ⁇ are obtained.
  • segmentation process part 32 calculates
  • the position of the divided particle group is first determined, and then the optimum mass of each divided particle is obtained. As described below, first, the divided particle group is first determined. It is also possible to determine the mass and then determine the optimal position of each divided particle.
  • the division processing unit 32 sets a constraint condition represented by the following expression (12), for example, that the total value of displacement from the position before the division of the particles becomes zero.
  • the division processing unit 32 defines energy J (#x i ) represented by the following equation (13).
  • a solution #x i the defined energy J (#x i) becomes an extreme value following equation (14), the optimum position of each particle after division.
  • the following expression (14) is established for arbitrary ⁇ x i and ⁇ , and the above expression (13) is an expression of quadratic or less for #x i and ⁇ . Therefore, linear simultaneous equations regarding #x i and ⁇ are obtained from the equation (14).
  • segmentation process part 32 calculates
  • requires the optimal mass after determining a position and the method which calculates
  • the dividing process may be performed using any technique.
  • the motion analysis unit 34 applies a physical model set in advance for each type of object such as a fluid, an elastic body, and a rigid body, and analyzes the motion of the object formed by the particles.
  • a physical model used by the motion analysis unit 34 will be exemplified.
  • the motion analysis unit 34 discretizes the motion equations of the incompressible fluid and the Navier-Stokes equations of the following equations (15) to (17) (1) By using the equation of motion, the time derivative of acceleration and temperature is obtained.
  • u is internal energy
  • ⁇ is a viscous stress tensor
  • ⁇ is a heat conduction coefficient
  • T is a temperature.
  • the motion analysis unit 34 calculates the Cauchy stress ⁇ by solving the following equations (18) to (20) from the deformation gradient tensor obtained by the above equation (4).
  • C is a Cauchy-Green deformation tensor
  • E corresponds to a Green-Lagrange strain tensor
  • j is a volume change rate
  • Q is an elastic potential.
  • the motion analysis unit 34 uses, as the elastic potential Q, a linear spring potential expressed by the following formula (21) or a Mooney-Rivlin potential expressed by the following formulas (22) to (24).
  • k is a spring constant
  • tr is a trace of matrix operation.
  • the motion analysis unit 34 obtains a motion equation by the following equation (25).
  • the motion analysis unit 34 can analyze the motion of an object by applying a well-known Newton equation of motion in addition to a physical model corresponding to the type of the object.
  • FIG. 14 is a flowchart showing the flow of overall processing executed by the motion analysis apparatus 1.
  • the data input unit 30 receives input of data relating to particles constituting the object, and stores the data in the memory device 18 (S100).
  • the input data includes the upper limit loop number N in addition to the data described above.
  • the upper limit loop number N is a value that designates how many times the loop processing of S110 to S170 is performed.
  • the number of loop processes is counted as a process loop number n every time the loop process is performed.
  • the initial value of the processing loop number n is 0.
  • the division processing unit 32 takes out each piece of particle data stored in the memory device 18, for example, in the order of storage (S110), and determines whether or not division processing is necessary (S120).
  • FIG. 15 is a flowchart showing the flow of processing when determining whether or not division processing is necessary.
  • the division processing unit 32 first generates a restandardized matrix L (#r) (S121), and calculates a determinant detL (#r) (S122). Next, the division processing unit 32 determines particle division by comparing the determinant detL (#r) with the threshold value Ref (S123). When the determinant detL (#r) exceeds the threshold value Ref, the division processing unit 32 determines that particle division is necessary (S124), and when the determinant detL (#r) is equal to or less than the threshold value Ref. Determines that the particle division is not required (S125).
  • the division processing unit 32 divides the particle and stores data such as the position and mass of the divided particle in the memory device 18 (S130, S140). ). On the other hand, the division processing unit 32 does not execute the division processing when it is determined in S120 to S125 that the division processing is unnecessary (S130).
  • FIG. 16 is a flowchart showing the flow of division processing.
  • the division processing unit 32 determines the number of divided particles based on the distribution of other particles present around the particles to be divided (S141). Next, the division processing unit 32 determines the position and influence radius of each particle after division (S142), and determines the optimum mass of each particle after division (S143). Then, the division processing unit 32 stores the determined data such as the position and mass of the divided particles in the memory device 18 (S144).
  • the division processing unit 32 determines whether or not the processing of S110 to S144 has been executed for all particles to be analyzed (S150). If the division processing unit 32 has not executed the processing of S110 to S144 for all the particles to be analyzed, the process returns to S110. When the division processing unit 32 executes the processing of S110 to S144 for all the particles to be analyzed, the process proceeds to S160.
  • the motion analysis unit 34 applies the physical model corresponding to the type of the object to each particle data stored in the memory device 18 to calculate the particle acceleration (S160).
  • the motion analysis unit 34 integrates the acceleration over time to calculate the position and velocity of the particles after the lapse of the specified time ⁇ t (S170), and outputs the calculation result as an external file as necessary (S180). Further, the motion analysis unit 34 may display the calculation result of each routine on the display device 24.
  • the motion analysis unit 34 increments the processing loop number n by 1 (S190), and determines whether the processing loop number n is equal to N (S200). When the processing loop number n is different from N, the process returns to S110, and when the processing loop number n is equal to N, this flow is finished.
  • the motion analysis apparatus 1 of the present embodiment even if the fluid or elastic body, which is the object to be analyzed, is greatly deformed by crushing or the like and the particles constituting the object move and become partially sparse, The number of particles can be increased by dividing the particles.
  • a part that is partially sparse is a part where the number of other particles existing around the particle is less than a predetermined number as described above, or a part where a predetermined uneven distribution relationship exists between the particle and the other particle. .
  • the motion analysis apparatus 1 can maintain an environment in which a sufficient number of other particles exist around the particle, and accurately calculate the force (repulsive force, shear force, etc.) to be applied to the particle. be able to. Therefore, the motion analysis apparatus 1 can analyze the motion of the object more accurately.
  • the processing for calculating acceleration and the like for each particle is increased by a factor of 6, which may result in excessive calculation load and calculation time.
  • deformation of a fluid or an elastic body may be difficult to predict, no matter how densely the particles are arranged, it cannot be said that there are no places where the particle distribution is sparse.
  • the motion analysis apparatus 1 since the motion analysis apparatus 1 according to the present embodiment performs the dividing process on necessary portions, the total calculation load can be reduced as compared with the case where many particles are arranged from the beginning. As a result, the calculation speed can be improved. Furthermore, since the particles can be divided as necessary, even if the fluid or the elastic body continues to be deformed, it is possible to continuously eliminate the sparse distribution of the particles.
  • the motion of the object can be analyzed more accurately.
  • the division processing unit 32 exemplifies performing the particle division determination by comparing the determinant detL (#r) of the restandardized matrix L (#r) with the threshold value Ref.
  • the method is not limited to this.
  • the division processing unit 32 determines whether or not the number of other particles existing within a predetermined range centered on the particle is less than a predetermined number, or the particle and the other particle are substantially on the same plane or substantially the same size.
  • the case of being on a line may be determined by vector calculation or the like, and particle division may be determined.
  • the present invention can be used in the construction industry, the manufacturing industry, the computer software industry, etc., by analyzing the motion of an object by the particle method. For example, by applying the present invention to fluid analysis of river or sea water, it can be used for planning disaster prevention. Further, by analyzing the casting process according to the present invention, it can be used for designing products such as metal products. In addition, by applying the present invention to the elastic body, the shape of the sealing gel can be appropriately determined during product design. Moreover, the program itself using the present invention can be sold.
  • motion analyzer 10
  • Reference Signs List 12
  • drive device 14
  • storage medium 16
  • auxiliary storage device 18
  • memory device 20
  • interface device 22
  • input device 24
  • display device 30
  • data input unit 32
  • division processing unit 34 motion analysis unit

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

　物体を表現する複数の粒子を表すデータを記憶する記憶部と、複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、着目粒子を分割した複数の分割粒子を表すデータを記憶部に記憶する分割処理部と、記憶部が記憶する複数の分割粒子を表すデータを用いて、物体の運動を物体の種類に応じて解析する運動解析部と、を有することを特徴とする運動解析装置。

Description

運動解析装置、運動解析方法及び運動解析プログラム
　本発明は、運動解析装置、運動解析方法及び運動解析プログラムに関する。
　従来、コンピュータ上で物体の運動を再現したシミュレーション結果を出力する装置が、幅広い分野で用いられている。これに関連し、物体を粒子の集合として認識し、各粒子の運動を解析することにより粒子が集合した物体の運動を解析するための粒子法という手法が知られている。
　代表的な粒子法としては、ＳＰＨ（Smoothed　Particle　Hydrodynamics）法やＭＰＳ（Moving　Particle　Semi-implicit）法等が知られている（例えば、非特許文献１、２参照）。これらの手法では、ある粒子から影響半径と呼ばれる距離ｈ以内に存在する粒子を近傍粒子として扱い、影響半径ｈ以内の粒子間には反発力が作用するものとして粒子の運動を解析する。
　次式（１）は、ＳＰＨ法における運動方程式を離散化した式の一例である。式中、下付添え字ａ、ｂは、粒子を識別するための識別子である。また、#ｒａ、#ｖａ、ρａ、Ｐａ、ｍａ、Пａｂはそれぞれ粒子ａの位置ベクトル、速度ベクトル、密度、圧力、質量、粘性応力テンソルである。なお、「#」は、後に続くアルファベットがベクトルであることを示している。
Figure JPOXMLDOC01-appb-M000002
　上式（１）において、Ｗはカーネル関数である。カーネル関数Ｗは、粒子の分布から連続場を構成するのに用いられる。カーネル関数Ｗは、具体的には、例えば次式（２）で示す三次のスプライン関数が用いられる。式中、ｈは粒子間の影響半径であり、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。βはカーネル関数の全空間積分量が１になるように調整された値であり、二次元の場合は０．７πｈ２、三次元の場合はπｈ３と決められる。
Figure JPOXMLDOC01-appb-M000003
　また、上式（１）は、流体の粘性を無視して完全流体と扱うことができるような場合には、粘性応力テンソルПａｂの項を無視することにより、次式（３）のように簡易に表すこともできる。
Figure JPOXMLDOC01-appb-M000004
　粒子法を用いた流動解析を行う方法等についての発明が開示されている（例えば、特許文献１参照）。この方法では、計算機毎に異なる粒子径の粒子を取り扱っており、各計算機は、担当する解析領域内に存在する粒子の粒子径が上限値を上回っている場合に、当該粒子を複数個の粒子に分割する。
特開２００７－１４０９９３号公報
J.J.Monaghan,　"Smoothed Particle Hydrodynamics",　Annu.　Rev.　Astron.　Astrophys.　30:　543-74（1992） S.Koshizuka,　A.Nobe　and　Y.Oka,　"Numerical　Analysis　of　Breaking　Waves　Using　The　Moving　Particle　Semi-Implicit　Method",　International　Journal　for　Numerical　Method　in　Fluids,　26:751-769（1998）
　ところで、粒子法を用いて流体や弾性体の運動を解析する際には、流体や弾性体が破砕やせん断等により大きく変形すると、物体を構成する粒子が移動して疎になる箇所が部分的に出現し得る。この結果、粒子の周囲に十分な数の他の粒子が存在せず、粒子に作用すべき力（反発力やせん断力等）が正確に算出されないため、運動解析結果が現実のものから解離してしまう場合がある。
　１つの側面では、本発明は、物体の運動の解析の正確性の向上を図る目的とする。
　上記目的を達成するための本発明の一態様は、
　物体を表現する複数の粒子を表すデータを記憶する記憶部と、
　前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、
　前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、
　を有することを特徴とする運動解析装置である。
　１実施態様によれば、物体の運動の解析の正確性の向上を図ることができる。
本発明の第１実施例に係る運動解析装置１のハードウエア構成例である。 本発明の一実施例に係る運動解析装置１の機能構成例である。 粒子ａの周囲に存在する他の粒子ｂが所定個数未満である状況を示す図である。 粒子ａと他の粒子ｂが略同一平面上に存在する状況を示す図である。 粒子ａと他の粒子ｂが略同一直線上に存在する状況を示す図である。 略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。 粒子ａを６個の粒子ａ＿１～ａ＿６に分割する様子を示す図である。 粒子ａを２個の粒子ａ＿１、ａ＿２に分割する様子を示す図である。 粒子ａを４個の粒子ａ＿１～ａ＿４に分割する様子を示す図である。 カーネル関数Ｗによって表現される密度ρ*を概念的に示す図である。 運動解析装置１によって実行される全体的な処理の流れを示すフローチャートである。 分割処理が必要か否かを判定する際の処理の流れを示すフローチャートである。 分割処理の流れを示すフローチャートである。
　以下、本発明を実施するための形態について、添付図面を参照しながら実施例を挙げて説明する。
　なお、図面及び数式内では原則としてベクトル表記を採用するが、文章中では、「#」が、後に続くアルファベットがベクトルであることを示すものとする。
　以下、本発明の一実施例に係る運動解析装置１及び運動解析方法、並びに運動解析プログラムについて説明する。運動解析装置１は、仮想的な三次元空間又は二次元平面上で、流体や弾性体、剛体等の物体の運動を粒子法に基づき解析する装置である。
　［ハードウエア構成］
　図１は、本発明の第１実施例に係る運動解析装置１のハードウエア構成例である。運動解析装置１は、例えば、ＣＰＵ（Central　Processing　Unit）１０と、ドライブ装置１２と、補助記憶装置１６と、メモリ装置１８と、インターフェース装置２０と、入力装置２２と、ディスプレイ装置２４と、を備える情報処理装置である。これらの構成要素は、バスやシリアル回線等を介して接続されている。
　ＣＰＵ１０は、例えば、プログラムカウンタや命令デコーダ、各種演算器、ＬＳＵ（Load　Store　Unit）、汎用レジスタ等を有するプロセッサ等の演算処理装置である。
　ドライブ装置１２は、記憶媒体１４からプログラムやデータを読み込み可能な装置である。プログラムを記憶した記憶媒体１４がドライブ装置１２に装着されると、プログラムが記憶媒体１４からドライブ装置１２を介して補助記憶装置１６にインストールされる。記憶媒体１４は、例えば、ＣＤ－ＲＯＭ、ＤＶＤディスク、ＵＳＢメモリ等の可搬型の記憶媒体である。また、補助記憶装置１６は、例えば、ＨＤＤ（Hard　Disk　Drive）やフラッシュメモリである。
　プログラムのインストールは、上記のように記憶媒体１４を用いる他、インターフェース装置２０がネットワークを介して他のコンピュータよりダウンロードし、補助記憶装置１６にインストールすることによって行うこともできる。ネットワークは、インターネット、ＬＡＮ（Local　Area　Network）、無線ネットワーク等である。また、プログラムは、情報処理装置の出荷時に、予め補助記憶装置１６やＲＯＭ（Read　Only　Memory）等に格納されていてもよい。
　このようにしてインストール又は予め格納されたプログラムをＣＰＵ１０が実行することにより、図１に示す態様の情報処理装置が、本実施例の運動解析装置１として機能することができる。
　メモリ装置１８は、例えば、ＲＡＭ（Random　Access　Memory）やＥＥＰＲＯＭ（Electrically　Erasable　and　Programmable　Read　Only　Memory）である。インターフェース装置２０は、上記ネットワークとの接続等を制御する。
　入力装置２２は、例えば、キーボード、マウス、タッチパッド、タッチパネルやマイク等である。また、ディスプレイ装置２４は、例えば、ＬＣＤ（Liquid　Crystal　Display）やＣＲＴ（Cathode　Ray　Tube）等の表示装置である。運動解析装置１は、ディスプレイ装置２４の他、プリンタ、スピーカ等、他の種類の出力装置を備えてもよい。
　［機能構成］
　図２は、本発明の一実施例に係る運動解析装置１の機能構成例である。運動解析装置１は、データ入力部３０と、分割処理部３２と、運動解析部３４と、を備える。これらの機能ブロックは、補助記憶装置１６等に格納されたプログラム・ソフトウエアをＣＰＵ１０が実行することにより機能する。なお、これらの機能ブロックが明確に分離したプログラムによって実現される必要はなく、サブルーチンや関数として他のプログラムから呼び出されるものであってもよい。また、機能ブロックの一部が、ＩＣ（Integrated　Circuit）やＦＰＧＡ（Field　Programmable　Gate　Array）等のハードウエア手段であっても構わない。
　〔データ入力〕
　データ入力部３０は、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置１８に格納する。入力データは、流体や弾性体等をモデル化した粒子の集合データであり、具体的には、粒子の中心座標、速度、影響半径、密度、質量、変形勾配テンソル、材料物性、温度等を含む。
　これらのデータは、予め補助記憶装置１６に格納されていたものであってもよいし、ユーザが入力装置２２を用いて入力したものであってもよい。また、これらのデータは、記憶媒体１４から補助記憶装置１６にインストールされたものであってもよいし、インターフェース装置２０によりネットワークから取得されたものであってもよい。
　ここで、変形勾配テンソルＦｉｊは、一般的には次式（４）で表される。式中、Ｘｊは粒子の位置間の距離のｊ成分、ｕｉは粒子の変位のｉ成分、Ｉは恒等テンソルである。ｉ成分及びｊ成分は、三次元であれば三次元座標を表すｘ、ｙ、ｚの３成分を表している。
Figure JPOXMLDOC01-appb-M000005
　粒子法における粒子ａの変形勾配テンソルは、粒子ａを中心として粒子ａとの距離が粒子ａの影響半径以内の領域に存在する粒子ｂからの影響を重ね合わせて、例えば次式（５）により表される。式中、ｄｕは粒子ａの変位を表している。以下、「粒子ａを中心として粒子ａとの距離が粒子ａの影響半径以内の領域に存在する」ことを、単に「粒子ａの周囲に存在する」と表記する。
Figure JPOXMLDOC01-appb-M000006
　なお、入力データが粒子を用いて表現されていない場合に備え、データ入力部３０は、物体のデータから粒子を用いたデータを設定する機能を有してもよい。
　分割処理部３２は、メモリ装置１８に格納された各粒子の位置、影響半径等の各データを考慮し、粒子の分割の判定及び分割処理を行う。
　〔分割の判定〕
　粒子の分割処理は、ある粒子ａが、当該粒子ａの周囲に存在する粒子から、本来受けるべき反発力が正確に作用しなくなる状況で行われる。
　係る状況は、粒子ａと他の粒子ｂの間に所定の偏在関係が存在する状況である。「所定の偏在関係」とは、例えば（Ａ）粒子ａの周囲に存在する他の粒子ｂが所定個数（例えば三次元ならば３個程度、二次元ならば２個程度）未満である状況、又は（Ｂ）粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが、略同一平面上又は略同一直線上に存在する関係である。
　図３は、上記（Ａ）の状況を示す図であり（粒子ｂが２個）、図４、５は、上記（Ｂ）の状況を示す図である。図３中、ｆは粒子ａが粒子ｂから受ける反発力を示している。また、図４において平面をＳと、図５において直線をＬとそれぞれ表記した。なお、運動解析装置１の扱う空間が二次元の場合、粒子ａと他の粒子ｂが略同一平面上に存在する状況は「所定の偏在関係」から除かれる。
　ここで、「略同一平面上に存在する」とは、例えば「粒子ａと、粒子ａの周辺に存在する任意の２個の粒子ｂを含む平面を基準平面Ｓとし、全ての粒子ｂについて基準平面Ｓからの距離Ｄ１が所定値Ｄｒｅｆ１以内である」ことと定義することができる。また、これに限らず、「略同一平面上に存在する」ことを、「粒子ａから全ての粒子ｂへのベクトルと基準平面Ｓのなす角度θ１が所定角度θｒｅｆ１以内である」ことと定義してもよい。図６、７は、略同一平面上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。
　また、「略同一直線上に存在する」とは、例えば「粒子ａと、粒子ａの周辺に存在する任意の１個の粒子ｂを含む直線を基準直線Ｌとし、全ての粒子ｂについて基準直線Ｌからの距離Ｄ２が所定値Ｄｒｅｆ２以内である」ことと定義することができる。また、これに限らず、「略同一直線上に存在する」ことを、「粒子ａから全ての粒子ｂへのベクトルと基準直線Ｌのなす角度θ２が所定角度θｒｅｆ２以内である」ことと定義してもよい。図８、９は、略同一直線上に存在すると定義される粒子ａと粒子ｂの位置関係を例示した図である。
　上記（Ａ）の状況では、粒子ａは周囲の粒子ｂから十分な反発力を受けないため、比較的自由に移動できる。また、上記（Ｂ）の状況では、粒子ａは基準平面Ｓ又は基準直線Ｌに直交する方向に関して比較的自由に移動できる。いずれの状況も流体や弾性体を構成する粒子の現実のふるまいとしては適切でないため、分割処理部３２は、粒子ａについて上記（Ａ）、（Ｂ）の状況が生じると、粒子ａを複数の粒子に分割して分割後の新たな粒子の位置及び質量等のデータをメモリ装置１８に格納する。
　分割処理部３２は、例えば次式（６）で表される再標準化行列Ｌ（#ｒ）を計算し、行列式ｄｅｔＬ（#ｒ）を求め、行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより、粒子の分割の判定を行う。具体的には、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆを超える場合には粒子の分割が必要であると判定し、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以下である場合には粒子の分割が不要であると判定する。なお、次式（６）は、運動解析装置１の扱う空間が三次元の場合の再標準化行列Ｌ（#ｒ）を示している。#ｒは粒子ａの位置ベクトルであり、（ｘ，ｙ，ｚ）を成分とする。
Figure JPOXMLDOC01-appb-M000007
　運動解析装置１の扱う空間が二次元の場合、再標準化行列Ｌは２×２の行列となり、次式（７）で表される。この場合、#ｒは（ｘ，ｙ）を成分とする。
Figure JPOXMLDOC01-appb-M000008
　上式（６）、（７）のいずれの場合にも、粒子ａの周囲に存在する粒子ｂが均等に分布しているならば、ｄｅｔＬ（#ｒ）は１に近い値となる。
　一方、前述した（Ａ）、（Ｂ）の状況下では、ｄｅｔＬ（#ｒ）は増加し、無限大までの値をとり得ることになる。従って、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以上である場合に、粒子ａの分割処理を行う。閾値Ｒｅｆは、上記のような行列式ｄｅｔＬ（#ｒ）の性質から、１を超え、２を超えない程度の範囲内で任意に設定される。
　なお、閾値Ｒｅｆは、「略同一平面上」や「略同一直線上」をどの程度まで含めるかを決定付けるパラメータである。閾値が大きい値であればあるほど、分割処理を行う可能性がより低くなるため、「略同一平面上」や「略同一直線上」はより狭い範囲を含むことになる。一方、閾値が小さい値であればあるほど、分割処理を行う可能性がより高くなるため、「略同一平面上」や「略同一直線上」はより広い範囲を含むことになる。
　〔分割処理〕
　分割処理部３２は、上記のような手法により粒子ａの分割処理が必要であると判定すると、以下のような手順で粒子ａの分割処理を実行する。
　まず、粒子ａを分割する個数ｉを決定する。個数ｉについては任意に定めてもよいが、以下のように、状況に応じた分割ルールを予め定めておくと好適である。
　粒子ａの周囲に存在する他の粒子ｂが所定個数未満である場合には、例えば粒子ａをｘ軸方向、ｙ軸方向、ｚ軸方向にそれぞれ分割して例えば６個の新たな粒子とする。図１０は、粒子ａを６個の粒子ａ＿１～ａ＿６に分割する様子を示す図である。
　また、粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが略同一平面上に存在する場合には、例えば粒子ａを前述した基準平面Ｓに直交する方向に分割して例えば２個の新たな粒子とする。図１１は、粒子ａを２個の粒子ａ＿１、ａ＿２に分割する様子を示す図である。
　また、粒子ａと粒子ａの周囲に存在する全ての他の粒子ｂが略同一直線上に存在する場合には、例えば粒子ａを前述した基準直線Ｌに直交する２方向に分割して例えば４個の新たな粒子とする。図１２は、粒子ａを４個の粒子ａ＿１～ａ＿４に分割する様子を示す図である。
　このように、粒子ａと粒子ａの周囲に存在する粒子ｂが、偏在関係すなわち略同一平面上や略同一直線状に存在するという関係を有する場合は、偏在方向すなわち平面に平行な方向や直線の方向に、直交する方向に粒子ａを分割すると好適である。
　分割処理部３２は、以下のように、分割された粒子間の距離、及び分割された粒子の影響範囲を決定する。分割処理部３２は、分割前の粒子の位置を#ｘ*とすると、当該粒子からの距離が影響半径以下の領域内で、前述した分割ルールに従って分割後の粒子の位置#ｘｉを決定する（ｉ＝１、２、…）。分割後の粒子の影響半径は、分割前の粒子の影響半径と同じものを用いてよい。
　なお、粒子の分割処理は、各粒子を独立に分割させるのではなく、２つの粒子を合計３つの粒子に分割するなどしても良い。
　ここで、粒子法に基づく解析の性質上、分割後の粒子群によって表される密度分布は、分割前の粒子の密度分布にできるだけ近づけることが好ましい。分割前の粒子#ｘ*の密度ρ*は、次式（８）で表される。式中、Ｍは分割前の粒子#ｘ*の質量であり、Ｗは前述のカーネル関数である。図１３は、カーネル関数Ｗによって表現される密度ρ*を概念的に示す図である。
Figure JPOXMLDOC01-appb-M000009
　一方、分割後の粒子群があらわす密度ρは、粒子法の重ねあわせの式により、次式（９）で表される。式中、ｍｉ（ｉ＝１、２、…）は分割後の各粒子の質量であり、制約条件としてＭ＝Σｉｍｉを満たす必要がある。
Figure JPOXMLDOC01-appb-M000010
　分割処理部３２は、上記制約条件下で密度ρを密度ρ*に近づけるための質量ｍｉの最適解を求めるために、ラグランジュの未定乗数法（Lagrange　Multiplier）を用いる。まず、分割処理部３２は、次式（１０）で表されるエネルギーＪ（ｍｉ）を定義する。
Figure JPOXMLDOC01-appb-M000011
　そして、分割処理部３２は、定義したエネルギーＪ（ｍｉ）が極値となる次式（１１）の解ｍｉを、分割後の各粒子の最適な質量とする。具体的には、分割処理部３２は、次式（１１）が任意のδｍｉ、δλに対して成立し、上式（１０）がｍｉ、λに対して二次以下の式であるため、式（１１）よりｍｉ、λに関する線形連立方程式を得る。そして、分割処理部３２は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。
Figure JPOXMLDOC01-appb-M000012
　以上説明した手法は、まず分割後の粒子群の位置を決定し、次に分割された各粒子の最適な質量を求めるものであるが、以下に説明するように、まず分割後の粒子群の質量を決定し、次に分割された各粒子の最適な位置を求めることもできる。
　この場合、分割処理部３２は、例えば粒子の分割前の位置からの変位の合計値が０になるという、次式（１２）で表される制約条件を設定する。
Figure JPOXMLDOC01-appb-M000013
　更に、分割処理部３２は、次式（１３）で表されるエネルギーＪ（#ｘｉ）を定義する。
Figure JPOXMLDOC01-appb-M000014
　そして、分割処理部３２は、定義したエネルギーＪ（#ｘｉ）が極値となる次式（１４）の解#ｘｉを、分割後の各粒子の最適な位置とする。具体的には、分割処理部３２は、次式（１４）が任意のδｘｉ、δλに対して成立し、上式（１３）が#ｘｉ、λに対して二次以下の式であるため、式（１４）より#ｘｉ、λに関する線形連立方程式を得る。そして、分割処理部３２は、この線形連立方程式を解くことにより、分割前の粒子により表される密度と、分割後の粒子群により表される密度の誤差を最小にする最適値を求める。
Figure JPOXMLDOC01-appb-M000015
　なお、位置を決定してから最適な質量を求める手法と、質量を決定してから最適な位置を求める手法を例示したが、分割前後の質量を一定にしつつ分割前後の密度の変化等を小さくする手法であれば、如何なる手法を用いて分割処理を行ってもよい。
　〔運動解析〕
　運動解析部３４は、流体、弾性体、剛体等の物体の種類毎に予め設定されている物理モデルを適用して、粒子が構成する物体の運動を解析する。以下、運動解析部３４によって用いられる物理モデルを例示する。
　例えば、運動解析部３４は、物体が流体である場合は、次式（１５）～（１７）の非圧縮性流体の運動方程式、及びナヴィエ・ストークス方程式を離散化して求められる上式（１）の運動方程式を用いることにより、加速度や温度の時間微分を求める。式中、ｕは内部エネルギーであり、Пは粘性応力テンソルであり、κは熱伝導係数であり、Ｔは温度である。
Figure JPOXMLDOC01-appb-M000016
　運動解析部３４は、物体が弾性体である場合は、上式（４）により求められる変形勾配テンソルから、次式（１８）～（２０）を解くことにより、コーシー応力σを算出する。式中、ＣはCauchy-Green変形テンソルであり、ＥはGreen-Lagrange歪テンソルに相当し、ｊは体積変化率であり、Ｑは弾性ポテンシャルである。
Figure JPOXMLDOC01-appb-M000017
　運動解析部３４は、上記弾性ポテンシャルＱとして、次式（２１）で表される線形バネのポテンシャルや、次式（２２）～（２４）で表されるMooney-Rivlinポテンシャルを用いる。式中、ｋはバネ定数であり、ｔｒは行列演算のトレースである。
Figure JPOXMLDOC01-appb-M000018
　運動解析部３４は、応力テンソルであるコーシー応力σを算出すると、次式（２５）によって運動方程式を得る。
Figure JPOXMLDOC01-appb-M000019
　なお、運動解析部３４は、これらの物体の種類に応じた物理モデルに加えて、周知のニュートンの運動方程式等を適用して物体の運動を解析することができる。
　［処理フロー］
　以下、運動解析装置１によって実行される処理の流れをフローチャートに則して説明する。図１４は、運動解析装置１によって実行される全体的な処理の流れを示すフローチャートである。
　まず、データ入力部３０が、物体を構成する粒子に関するデータの入力を受け付け、メモリ装置１８に格納する（Ｓ１００）。入力されるデータは、前述したデータに加え、上限ループ数Ｎが含まれる。上限ループ数Ｎは、Ｓ１１０～Ｓ１７０のループ処理を何回行うかを指定する値である。ループ処理の回数は、ループ処理を行なう毎に処理ループ数ｎとしてカウントされる。処理ループ数ｎの初期値は０である。
　分割処理部３２は、メモリ装置１８に格納された粒子の各データを、例えば格納順に一つ取り出し（Ｓ１１０）、分割処理が必要か否かを判定する（Ｓ１２０）。
　図１５は、分割処理が必要か否かを判定する際の処理の流れを示すフローチャートである。分割処理部３２は、まず、再標準化行列Ｌ（#ｒ）を生成し（Ｓ１２１）、行列式ｄｅｔＬ（#ｒ）を算出する（Ｓ１２２）。次に、分割処理部３２は、行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより、粒子の分割の判定を行う（Ｓ１２３）。分割処理部３２は、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆを超える場合には粒子の分割が必要であると判定し（Ｓ１２４）、行列式ｄｅｔＬ（#ｒ）が閾値Ｒｅｆ以下である場合には粒子の分割は不要であると判定する（Ｓ１２５）。
　分割処理部３２は、上記Ｓ１２０～Ｓ１２５において分割処理が必要と判定された場合は、当該粒子を分割して分割後の粒子の位置及び質量等のデータをメモリ装置１８に格納する（Ｓ１３０、Ｓ１４０）。一方、分割処理部３２は、Ｓ１２０～Ｓ１２５において分割処理が不要と判定された場合は、分割処理を実行しない（Ｓ１３０）。
　図１６は、分割処理の流れを示すフローチャートである。分割処理部３２は、分割対象となる粒子の周囲に存在する他の粒子の分布に基づき分割後の粒子の個数を決定する（Ｓ１４１）。次に、分割処理部３２は、分割後の各粒子の位置及び影響半径を決定し（Ｓ１４２）、分割後の各粒子の最適な質量を決定する（Ｓ１４３）。そして、分割処理部３２は、決定した分割後の粒子の位置及び質量等のデータをメモリ装置１８に格納する（Ｓ１４４）。
　次に、分割処理部３２は、解析対象の全ての粒子についてＳ１１０～Ｓ１４４の処理を実行したか否かを判定する（Ｓ１５０）。分割処理部３２が解析対象の全ての粒子についてＳ１１０～Ｓ１４４の処理を実行していない場合は、Ｓ１１０に戻る。分割処理部３２が解析対象の全ての粒子についてＳ１１０～Ｓ１４４の処理を実行すると、Ｓ１６０に進む。
　運動解析部３４は、メモリ装置１８に格納された粒子の各データに対して、物体の種類に応じた物理モデルを適用して、粒子の加速度を算出する（Ｓ１６０）。
　そして、運動解析部３４は、加速度を時間積分して規定時間Δｔ経過後の粒子の位置、及び速度を算出し（Ｓ１７０）、必要に応じて算出結果を外部ファイルとして出力する（Ｓ１８０）。また、運動解析部３４は、各ルーチンの算出結果をディスプレイ装置２４に表示させてもよい。
　次に、運動解析部３４は、処理ループ数ｎを１インクリメントし（Ｓ１９０）、処理ループ数ｎがＮと等しいか否かを判定する（Ｓ２００）。処理ループ数ｎがＮと異なる場合はＳ１１０に戻り、処理ループ数ｎがＮと等しい場合は本フローを終了する。
　［まとめ］
　本実施例の運動解析装置１は、解析対象の物体である流体や弾性体が破砕等により大きく変形し、物体を構成する粒子が移動して疎になる箇所が部分的に出現したとしても、粒子を分割して粒子数を増やすことができる。部分的に疎になる箇所とは、前述のように粒子の周囲に存在する他の粒子が所定個数未満である箇所や、粒子と他の粒子の間に所定の偏在関係が存在する箇所である。
　この結果、運動解析装置１は、粒子の周囲に十分な数の他の粒子が存在する環境を維持することができ、粒子に作用すべき力（反発力やせん断力等）を正確に算出することができる。従って、運動解析装置１は、より正確に物体の運動を解析することができる。
　なお、流体や弾性体が大きく変形しても粒子の分布が疎になる箇所が出現しないように、多くの粒子を最初から配置することも考えられる。しかしながら、この場合、演算対象となる粒子の数が増加することにより、演算負荷や演算時間が過大となってしまう場合がある。例えば、前述のように粒子をｘ軸方向、ｙ軸方向、ｚ軸方向にそれぞれ分割して例えば６個の新たな粒子とする処理を、最初から各粒子に対して実行したとすると、粒子の分布が密となり、部分的に疎になる箇所は生じにくくなる。反面、各粒子について加速度等を演算する処理が６倍に増加するため、演算負荷や演算時間が過大となる可能性がある。また、流体や弾性体の変形は予測が困難な場合があるため、どれだけ密に粒子を配置したとしても、粒子の分布が疎になる箇所が出現しないとは言い切れない。
　これに対し、本実施例の運動解析装置１は、必要な箇所について分割処理を行うため、多くの粒子を最初から配置する場合に比して、トータルの演算負荷を軽減することができる。この結果、演算速度を向上させることができる。更に、必要に応じて粒子を分割することができるため、流体や弾性体が変形し続けたとしても、粒子の分布が疎になるのを持続的に解消することができる。
　以上説明した本実施例の運動解析装置１によれば、より正確に物体の運動を解析することができる。
　以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。
　例えば、分割処理部３２は、再標準化行列Ｌ（#ｒ）の行列式ｄｅｔＬ（#ｒ）を閾値Ｒｅｆと比較することにより粒子の分割の判定を行うことを例示したが、粒子の分割の判定手法は、これに限られない。分割処理部３２は、粒子の配置から、当該粒子を中心とする所定範囲内に存在する他の粒子が所定個数未満である場合や、当該粒子と他の粒子が略同一平面上又は略同一直線上に存在する場合をベクトル演算等によって判別し、粒子の分割の判定を行ってもよい。
　本発明は、物体の運動を粒子法により解析することを通じて、建築業や製造業、コンピュータソフトウエア産業等に利用することができる。例えば、本発明を河川や海の水の流体解析に対して適用することで、防災の計画立案に役立てることができる。また、本発明により鋳造過程を解析することで金属製品等の製品設計に用いることが出来る。また、本発明を弾性体に対して適用することで、製品設計の際に封止ゲルの形状などを適切に決定することが出来る。また、本発明を利用したプログラム自体を販売することが出来る。
　　１　運動解析装置
　１０　ＣＰＵ
　１２　ドライブ装置
　１４　記憶媒体
　１６　補助記憶装置
　１８　メモリ装置
　２０　インターフェース装置
　２２　入力装置
　２４　ディスプレイ装置
　３０　データ入力部
　３２　分割処理部
　３４　運動解析部

Claims (9)

  1. 　物体を表現する複数の粒子を表すデータを記憶する記憶部と、
    　前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶する分割処理部と、
    　前記記憶部が記憶する前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析する運動解析部と、
    　を有することを特徴とする運動解析装置。
  2. 　前記運動解析装置において、
    　前記所定の偏在関係は、
    　前記着目粒子を中心とした所定範囲内に存在する他の粒子が所定個数未満の関係であることを特徴とする請求項１記載の運動解析装置。
  3. 　前記運動解析装置において、
    　前記所定の偏在関係は、
    　前記着目粒子と前記他の粒子が略同一平面上又は略同一直線上に存在する関係であることを特徴とする請求項１記載の運動解析装置。
  4. 　前記運動解析装置において、
    　前記分割処理部は、
    　前記着目粒子と前記他の粒子との分布を表す関数の値を成分とする再標準化行列の行列式の値が所定値以上である場合、前記着目粒子を分割した複数の分割粒子を表すデータを前記記憶部に記憶することを特徴とする請求項１記載の運動解析装置。
  5. 　前記運動解析装置において、
    　前記再標準化行列は、
    Figure JPOXMLDOC01-appb-M000001
     　であり、
    　前記数式中、Ｗはカーネル関数、ｈは前記着目粒子の影響半径、ｒ（ベクトル表記）は前記着目粒子の位置ベクトル、（ｘ，ｙ，ｚ）は前記着目粒子のベクトル成分、ｒｂ（ベクトル表記）は前記他の粒子の位置ベクトル、（ｘｂ，ｙｂ，ｚｂ）は前記他の粒子のベクトル成分であることを特徴とする請求項４記載の運動解析装置。
  6. 　前記運動解析装置において、
    　前記分割処理部は、
    　前記着目粒子と前記他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子と前記他の粒子の偏在方向と直交する方向に、前記着目粒子を分割することを特徴とする請求項１記載の運動解析装置。
  7. 　前記運動解析装置において、
    　前記分割処理部は、
    　前記複数の分割粒子の位置を決定した後、前記複数の分割粒子の分割前後の密度関数の差分が小さくなるように、前記複数の分割粒子の質量を決定することを特徴とする請求項１記載の運動解析装置。
  8. 　物体を表現する複数の粒子の運動を解析する運動解析方法において、
    　コンピュータが、
    　前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割し、
    　前記複数の分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析することを特徴とする運動解析方法。
  9. 　物体を表現する複数の粒子の運動を解析する運動解析プログラムにおいて、
    　コンピュータに、
    　前記複数の粒子に含まれる着目粒子を中心とした所定範囲内に存在する他の粒子の間に所定の偏在関係が存在する場合、前記着目粒子を表すデータを複数の分割粒子を表すデータに分割させ、
    　前記分割粒子を表すデータを用いて、前記物体の運動を前記物体の種類に応じて解析させることを特徴とする運動解析プログラム。
PCT/JP2011/070764 2011-09-12 2011-09-12 運動解析装置、運動解析方法及び運動解析プログラム WO2013038476A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2013533361A JP5761355B2 (ja) 2011-09-12 2011-09-12 運動解析装置、運動解析方法及び運動解析プログラム
PCT/JP2011/070764 WO2013038476A1 (ja) 2011-09-12 2011-09-12 運動解析装置、運動解析方法及び運動解析プログラム
US14/203,701 US20140195212A1 (en) 2011-09-12 2014-03-11 Motion analysis apparatus and motion analysis method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2011/070764 WO2013038476A1 (ja) 2011-09-12 2011-09-12 運動解析装置、運動解析方法及び運動解析プログラム

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/203,701 Continuation US20140195212A1 (en) 2011-09-12 2014-03-11 Motion analysis apparatus and motion analysis method

Publications (1)

Publication Number Publication Date
WO2013038476A1 true WO2013038476A1 (ja) 2013-03-21

Family

ID=47882744

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2011/070764 WO2013038476A1 (ja) 2011-09-12 2011-09-12 運動解析装置、運動解析方法及び運動解析プログラム

Country Status (3)

Country Link
US (1) US20140195212A1 (ja)
JP (1) JP5761355B2 (ja)
WO (1) WO2013038476A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102459848B1 (ko) * 2014-10-24 2022-10-27 삼성전자주식회사 입자에 기반하여 대상 객체를 고속으로 모델링하는 방법 및 장치
CN108171800B (zh) * 2016-12-07 2021-05-25 北京乌有园艺术有限公司 创造具有无穷造型的制品的方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009134504A (ja) * 2007-11-30 2009-06-18 Fuji Xerox Co Ltd 粒子挙動解析装置、プログラム
JP2010049561A (ja) * 2008-08-22 2010-03-04 Toshiba Corp 流動解析方法、流動解析装置、及び流動解析プログラム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7755765B2 (en) * 2003-03-17 2010-07-13 Massachusetts Institute Of Technology Method and apparatus for inertial sensing via measurement of trapped orbit dynamics

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009134504A (ja) * 2007-11-30 2009-06-18 Fuji Xerox Co Ltd 粒子挙動解析装置、プログラム
JP2010049561A (ja) * 2008-08-22 2010-03-04 Toshiba Corp 流動解析方法、流動解析装置、及び流動解析プログラム

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
MASAYUKI TANAKA: "multi-resolution MPS method", PROCEEDINGS OF THE CONFERENCE ON COMPUTATIONAL ENGINEERING AND SCIENCE, vol. 2009, 28 January 2009 (2009-01-28) *
YOSHIKAZU SHIMIZU: "Ryushi-ho o Mochiita Asshukusei Nensei Ryutai no Suchi Kaiseki ni Okeru Seido no Kenkyu", THE JAPAN SOCIETY OF MECHANICAL ENGINEERS HEISEI 13 NENDO PREPRINTS OF GRADUATION RESEARCH THESIS OF KANSAI STUDENT SECTION OF THE JAPAN, 21 March 2002 (2002-03-21), pages 7 - 10 *

Also Published As

Publication number Publication date
JP5761355B2 (ja) 2015-08-12
JPWO2013038476A1 (ja) 2015-03-23
US20140195212A1 (en) 2014-07-10

Similar Documents

Publication Publication Date Title
Shojaei et al. An adaptive multi-grid peridynamic method for dynamic fracture analysis
Qu et al. An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids
US11645433B2 (en) Computer simulation of physical fluids on irregular spatial grids stabilized for explicit numerical diffusion problems
Rodríguez et al. Simulation of metal cutting using the particle finite-element method and a physically based plasticity model
Cornejo et al. Combination of an adaptive remeshing technique with a coupled FEM–DEM approach for analysis of crack propagation problems
US9250259B2 (en) Object motion analysis apparatus, object motion analysis method, and storage medium
Chadil et al. Accurate estimate of drag forces using particle-resolved direct numerical simulations
EP2535828A1 (en) Method for simulating the loss tangent of rubber compound
Daxini et al. Parametric shape optimization techniques based on Meshless methods: A review
De Corato et al. Finite element formulation of fluctuating hydrodynamics for fluids filled with rigid particles using boundary fitted meshes
Musehane et al. Multi-scale simulation of droplet–droplet interaction and coalescence
JP5761355B2 (ja) 運動解析装置、運動解析方法及び運動解析プログラム
Hashemi et al. Direct numerical simulation of particle–fluid interactions: a review
Dulikravich et al. Inverse problems in aerodynamics, heat transfer, elasticity and materials design
Garg et al. Accelerated element-free Galerkin method for analysis of fracture problems
Roden et al. MultiShape: A spectral element method, with applications to dynamic density functional theory and PDE-constrained optimization
Poe et al. A low-dissipation second-order upwind flux formulation for simulation of complex turbulent flows
Amani et al. On estimating the interface normal and curvature in piecewise linear interface calculation‐volume of fluid approach for three‐dimensional arbitrary meshes
Thornburg Overview of the pettt workshop on mesh quality/resolution, practice, current research, and future directions
Kumaran Dense shallow granular flows
Esteban et al. A contact line force model for the simulation of drop impacts on solid surfaces using volume of fluid methods
Zhang et al. An integrated coupling framework for highly nonlinear fluid-structure problems
Gao et al. A coupled isogeometric/multi-sphere discrete element approach for the contact interaction between irregular particles and structures
Kor et al. Extension of the unified interpolation stencil for immersed boundary method for moving boundary problems
Saberi et al. The Influence of Nitsche Stabilization on Geometric Multigrid for the Finite Cell Method

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11872337

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2013533361

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 11872337

Country of ref document: EP

Kind code of ref document: A1