SIMULATION METHOD, SIMULATION APPARATUS, AND NON-TRANSITORY COMPUTER READABLE MEDIUM STORING PROGRAM

Abstract

Provided is a simulation method that represents an elastic analysis object as an aggregate of a plurality of virtual particles, and the method includes, when solving the equation of motion for each particle, regarding the analysis object as a rigid body, and calculating a force acting on the analysis object and a velocity of the analysis object, based on the force acting on each particle obtained at each time step and the velocity and the position of each particle, applying FIRE to a motion of the analysis object to calculate a force to be applied to the analysis object, obtaining a force to be additionally applied to each particle, by distributing the force to be applied to the analysis object to the plurality of particles, and solving the equation of motion, by adding the force to be additionally applied, to the force acting on each particle.

Claims

1. A simulation method that represents an elastic analysis object as an aggregate of a plurality of virtual particles, obtains a force acting on each particle at each time step, and updates a velocity and a position of each particle by solving an equation of motion for each particle, based on the force acting on each particle, the simulation method comprising: when solving the equation of motion for each particle, regarding the analysis object as a rigid body, and calculating a force acting on the analysis object and a velocity of the analysis object, based on the force acting on each particle obtained at each time step and the velocity and the position of each particle; applying FIRE to a motion of the analysis object to calculate a force to be applied to the analysis object; obtaining a force to be additionally applied to each particle, by distributing the force to be applied to the analysis object to the plurality of particles; and solving the equation of motion, by adding the force to be additionally applied, to the force acting on each particle.

2. The simulation method according to claim 1, wherein a translational force and a torque are calculated as the force acting on the analysis object, a translational velocity and an angular velocity are calculated as the velocity of the analysis object, and when applying the FIRE to the motion of the analysis object, the FIRE is applied to each of a translational motion and a rotational motion of the analysis object to calculate the translational force and the torque as the force to be applied to the analysis object.

3. A simulation apparatus that simulates a behavior of an elastic analysis object using a molecular dynamics method, the simulation apparatus comprising: an input unit to which a simulation condition is input; and a processing unit that performs analysis, based on the simulation condition input to the input unit, wherein the processing unit represents the analysis object as an aggregate of a plurality of virtual particles, obtains a force acting on each particle at each time step, updates a velocity and a position of each particle by solving an equation of motion for each particle, based on the force acting on each particle, and when solving the equation of motion for each particle, regards the analysis object as a rigid body, and calculates a force acting on the analysis object and a velocity of the analysis object, based on the force acting on each particle obtained at each time step and the velocity and the position of each particle; applies FIRE to a motion of the analysis object to calculate a force to be applied to the analysis object, obtains a force to be additionally applied to each particle, by distributing the force to be applied to the analysis object to the plurality of particles, and solves the equation of motion, by adding the force to be additionally applied, to the force acting on each particle.

4. The simulation apparatus according to claim 3, wherein the processing unit calculates a translational force and a torque, as the force acting on the analysis object, calculates a translational velocity and an angular velocity, as the velocity of the analysis object, and when applying the FIRE to the motion of the analysis object, applies the FIRE to each of a translational motion and a rotational motion of the analysis object to calculate the translational force and the torque as the force to be applied to the analysis object.

5. A non-transitory computer readable medium storing a program causing a computer to execute a process comprising: acquiring a simulation condition; representing an elastic analysis object as an aggregate of a plurality of virtual particles; obtaining a force acting on each particle at each time step; updating a velocity and a position of each particle by solving an equation of motion for each particle, based on the force acting on each particle; and when solving the equation of motion for each particle, regarding the analysis object as a rigid body, and calculating a force acting on the analysis object and a velocity of the analysis object, based on the force acting on each particle obtained at each time step and the velocity and the position of each particle; applying FIRE to a motion of the analysis object to calculate a force to be applied to the analysis object; obtaining a force to be additionally applied to each particle, by distributing the force to be applied to the analysis object to the plurality of particles; and solving the equation of motion, by adding the force to be additionally applied, to the force acting on each particle.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0036] FIG. 1A is a schematic diagram showing an analysis model in which an elastic analysis object is represented by an aggregate of a plurality of particles, and FIG. 1B is a schematic diagram showing a physical quantity defined for the analysis object regarded as a rigid body.

[0037] FIG. 2 is a block diagram of a simulation apparatus according to a present embodiment.

[0038] FIG. 3 is a flowchart of a simulation method according to the present embodiment.

[0039] FIG. 4 is a graph showing a temporal change of a torque applied to a fixing member that fixes a part of a region of one end surface of a rotating body.

[0040] FIG. 5 is a graph showing a temporal change of an angular velocity co when the rotating body is regarded as a rigid body.

[0041] FIGS. 6A and 6B are schematic diagrams of an analysis model to be analyzed by applying FIRE.

DETAILED DESCRIPTION

[0042] FIRE is applied to a single particle and cannot be applied to an aggregate of particles representing an elastic analysis object. It is desirable to provide a simulation method, a simulation apparatus, and a program that represent an elastic analysis object as an aggregate of a plurality of particles, and capable of obtaining a solution at high speed when solving the equation of motion for each particle to simulate the behavior of the analysis object.

[0043] An elastic analysis object is represented as an aggregate of a plurality of particles, and when solving the equation of motion for each particle to simulate the behavior of the analysis object, a solution can be obtained at high speed.

General FIRE

[0044] Before explaining the embodiment, FIRE, which is a method for obtaining a solution at high speed in a simulation using a molecular dynamics, will be described with reference to FIGS. 6A and 6B.

[0045] FIGS. 6A and 6B are schematic diagrams of an analysis model analyzed by applying FIRE. A particle 20 is disposed in the analysis space. The velocity of the particle 20 is marked as v, and the force acting on the particle 20 is marked as f. When applying FIRE, the following determination value P(t) is calculated for each time step.

P(t)=f(t).Math.v(t)   (1)

[0046] When the determination value P(t) is negative or zero, that is, when the angle between the velocity v and the force f is 90° or more as shown in FIG. 6A, the velocity v(t+Δt) of the particle 20 in the next state is determined as follows.

v(t+Δt)=0, (P(t)≤0)   (2)

[0047] In other words, when the force acts in the direction of decelerating the particle 20, the velocity of the particle 20 in the next state is set to zero.

[0048] When the determination value P(t) is positive, that is, when the angle between the velocity v and the force f is less than 90° as shown in FIG. 6B, the following equation of motion is solved for the particle 20 to obtain the velocity v(t+Δt) in the next state.

[00001]$\begin{array}{cc}\stackrel{.}{v}\left(t\right)=\frac{f\left(t\right)}{m}-\gamma \left(t\right).Math.v\left(t\right).Math..Math.\left\{\hat{v}\left(t\right)-\hat{f}\left(t\right)\right\},\left(P\left(t\right)>0\right)& \left(3\right)\end{array}$

[0049] Here, m is the mass of the particle 20, and γ is the attenuation constant. The hat attached to the vector means that it is a unit vector. Expression (3) means that the particle 20 is further accelerated in the direction of the force f.

FIRE by Reference Example

[0050] Next, FIRE by a reference example in which an elastic analysis object is applied to an analysis model represented by an aggregate of particles will be described.

[0051] In the FIRE according to the reference example, the determination value P(t) is defined by the following expression.

[00002]$\begin{array}{cc}P\left(t\right)=\underset{i}{\overset{N}{.Math.}}{f}_i\left(t\right).Math.v_i\left(t\right)& \left(4\right)\end{array}$

[0052] Here, f.sub.i is the force acting on the i-th particle, and v.sub.i is the velocity of the i-th particle. N is the number of particles, and Σ on the right side of Expression (4) means that all the particles constituting the analysis object are totaled.

[0053] When the determination value P(t) is negative or zero, the velocity v.sub.i(t+Δt) of the i-th particle 20 in the next state is determined as follows.

v.sub.i(t+Δt)=0, (P(t)≤0)   (5)

[0054] When the determination value P(t) is positive, the velocity v.sub.i(t+Δt) in the next state is obtained by solving the following equation of motion for the i-th particle 20.

[00003]$\begin{array}{cc}{\stackrel{.}{v}}_i\left(t\right)=\frac{f_i\left(t\right)}{m_i}-\gamma \left(t\right).Math.v_i\left(t\right).Math..Math.\left\{{\hat{v}}_i\left(t\right)-{\hat{f}}_i\left(t\right)\right\},\left(P\left(t\right)>0\right)& \left(6\right)\end{array}$

[0055] Here, m.sub.i is the mass of the i-th particle.

[0056] The present inventors have recognized that the FIRE according to the above reference example has the follows.

[0057] In the FIRE according to the reference example, as shown in Expression (4), the total value of the inner product of the force acting on each of the particles and the velocity of the particles is used as the determination value P(t). Therefore, even when the velocity of the center of gravity of the analysis object is zero, the determination value P(t) is not zero due to the internal oscillation, but has a certain value. When the temperature of the analysis object is high or the velocity of the analysis object approaches convergence, the condition that the determination value P(t) is negative is often satisfied. When this condition is satisfied, the velocity of each particle becomes zero, which causes a delay in convergence.

[0058] Further, when the determination value P(t) is positive, as shown in Expression (6), the equation of motion to which the force accelerating in the direction of the force is additionally applied to each particle is solved. The attenuation constant γ (t) needs to be set for each of the particles 20. Since the mass and the interaction with the surroundings are different between the particles 20, it is difficult to set the attenuation constant γ, and the calculation may fail.

Embodiment

[0059] The simulation method, simulation apparatus, and program according to one embodiment of the present invention will be described with reference to FIGS. 1A to 5.

[0060] FIG. 1A is a schematic diagram showing an analysis model in which an elastic analysis object 10 is represented by an aggregate of a plurality of particles 11. Note that FIG. 1A shows only a part of the particles 11 constituting the analysis object 10. The xyz orthogonal coordinate system is defined for the analysis space, the position vector of the i-th particle 11 is marked as r.sub.i, and the velocity of the i-th particle 11 is marked as v.sub.i.

[0061] In the embodiment, the elastic analysis object 10 is assumed to be one rigid body, and the physical quantity is defined for the analysis object 10 regarded as a rigid body.

[0062] FIG. 1B is a schematic diagram showing a physical quantity defined for the analysis object 10 regarded as a rigid body. The position vector of the center of gravity of the analysis object 10 is marked as R. The translational force acting on the center of gravity of the analysis object 10 is marked as F, and the torque when the center of gravity is the center of rotation is marked as T. The translational velocity of the center of gravity of the analysis object 10 is marked as V, and the angular velocity is marked as ω. For example, when the analysis object such as a metal member is deformed only slightly, the analysis object is regarded as a rigid body, and the translational force F, torque T, translational velocity V, and angular velocity ω applied to the analysis object can be defined.

[0063] FIG. 2 is a block diagram of the simulation apparatus according to the present embodiment. The simulation apparatus according to the present embodiment includes an input unit 40, a processing unit 41, an output unit 42, and a storage unit 43. Simulation conditions and the like are input from the input unit 40 to the processing unit 41. Further, various commands are input from the user to the input unit 40. Examples of the input unit 40 include a communication device, a removable media reading device, a keyboard, a pointing device, and the like.

[0064] The processing unit 41 executes a simulation by the molecular dynamics, based on the input simulation conditions and commands. Further, the simulation result is output to the output unit 42. Examples of the simulation result include information representing temporal changes of the positions of the particles 11 (FIG. 1A) constituting the simulation object. Examples of the processing unit 41 include a central processing unit (CPU) of a computer. A program for causing the computer to execute the simulation by the molecular dynamics is stored in the storage unit 43. Examples of the output unit 42 include a communication device, a removable media writing device, a display, and the like.

[0065] FIG. 3 is a flowchart of the simulation method according to the present embodiment. By executing the program stored in the storage unit 43 (FIG. 2) by the processing unit 41 (FIG. 2) , the process of each step of the simulation method is executed.

[0066] First, the processing unit 41 generates an analysis model in which the elastic analysis object 10 (FIG. 1A) is represented by an aggregate of a plurality of particles, based on the simulation conditions input from the input unit 40 (step S1). The simulation conditions include the shape of the analysis object 10, information that depends on the physical property value of the analysis object 10, the initial condition and boundary condition of simulation, a time step width, an analysis end condition, and the like. The information that depends on the physical property value of the analysis object 10 includes information that defines the interaction potential between particles in the analysis model.

[0067] Next, the equation of motion is solved for each particle 11 (FIG. 1A), and the velocity and position of each particle 11 are updated by one time step.

[0068] Next, the analysis object 10 (FIG. 1B) is regarded as a rigid body, and the force acting on the analysis object 10 and the velocity of the analysis object 10 are calculated based on the position and velocity of the particles 11 (step S3). The velocity of the analysis object 10 includes a translational velocity of the center of gravity and an angular velocity with the center of gravity as the center.

[0069] The translational force F acting on the analysis object 10 is calculated by the following expression.

[00004]$\begin{array}{cc}F=\underset{i}{\overset{N}{.Math.}}{f}_i& \left(7\right)\end{array}$

[0070] That is, the translational force F acting on the analysis object 10 is the vector sum of the forces f.sub.i acting on each of the particles 11.

[0071] The translational velocity V of the analysis object 10 is calculated by the following expression.

[00005]$\begin{array}{cc}V=\frac{1}{M}\underset{i}{\overset{N}{.Math.}}{m}_i{v}_i& \left(8\right)\end{array}$

[0072] Here, M is the sum of the masses of all the particles 11 constituting the analysis object 10. The translational velocity V of the analysis object 10 is a value obtained by averaging the velocity v.sub.i of the particles 11 constituting the analysis object 10 weighted by the mass m.sub.i of the particles 11.

[0073] The torque T=(T.sub.x, T.sub.y, T.sub.z) acting on the analysis object 10 is calculated by the following expression.

[00006]$\begin{array}{cc}{T}_x=\stackrel{N}{\underset{i}{.Math.}}\left(y_i{f}_{i,z}-z_i{f}_{i,y}\right)T_y=\underset{i}{\overset{N}{.Math.}}\left(z_i{f}_{i,x}-x_i{f}_{i,x}\right)T_z=\underset{i}{\overset{N}{.Math.}}\left(x_i{f}_{i,y}-y_i{f}_{i,x}\right)& \left(9\right)\end{array}$

[0074] Here, T.sub.x, T.sub.y, and T.sub.z are the x component, the y component, and the z component of the torque T acting on the analysis object 10, respectively. x.sub.i, y.sub.i, and z.sub.i are the x-coordinate, the y-coordinate, and the z-coordinate of the i-th particle 11 when the center of gravity of the analysis object 10 is the origin, respectively. f.sub.i,x, f.sub.i,y, and f.sub.i,z are the x component, the y component, and the z component of the force f.sub.i acting on the i-th particle 11, respectively.

[0075] The angular velocity ω=(ω.sub.x, ω.sub.y, and ω.sub.z) of the analysis object 10 is calculated by the following expression.

[00007]$\begin{array}{cc}{\omega }_x=\frac{1}{I_x}\underset{i}{\overset{N}{.Math.}}{m}_i\left(y_i^2+z_i^2\right)\left(y_i{v}_{i,z}-z_i{v}_{i,y}\right){\omega }_y=\frac{1}{I_y}\underset{i}{\overset{N}{.Math.}}{m}_i\left(z_i^2+x_i^2\right)\left(z_i{v}_{i,x}-x_i{v}_{i,z}\right){\omega }_z=\frac{1}{I_z}\underset{i}{\overset{N}{.Math.}}{m}_i\left(x_i^2+y_i^2\right)\left(x_i{v}_{i,y}-y_i{v}_{i,x}\right)& \left(10\right)\end{array}$

[0076] Here, ω.sub.x, ω.sub.y, and ω.sub.z are the x component, the y component, and the z component of the angular velocity ω, respectively. v.sub.i,x, v.sub.i,y, and v.sub.i,z are the x component, the y component, and the z component of the velocity of the i-th particle 11, respectively.

[0077] I.sub.x, I.sub.y, and I.sub.z are moments of inertia around the x-axis, y-axis, and z-axis, respectively, and are defined by the following expression.

[00008]$\begin{array}{cc}{I}_x=\underset{i}{\overset{N}{.Math.}}{m}_i\left(y_i^2+z_i^2\right)I_y=\underset{i}{\overset{N}{.Math.}}{m}_i\left(z_i^2+x_i^2\right)I_z=\underset{i}{\overset{N}{.Math.}}{m}_i\left(x_i^2+y_i^2\right)& \left(11\right)\end{array}$

[0078] After step S3 (FIG. 3), the processing unit 41 (FIG. 2) applies FIRE to the motion of the analysis object 10 so as to alleviate the oscillation of the physical quantity defining the motion of the analysis object 10 regarded as a rigid body. Thus, the force to be applied to the analysis object 10 regarded as a rigid body is calculated. More specifically, the translational force and the torque to be applied to the analysis object 10 are calculated (step S4).

[0079] Hereinafter, the process of step S4 will be described in detail.

[0080] The motion of the analysis object 10 is separated into translational motion and rotational motion, and FIRE is applied to each motion. First, the FIRE applied to the translational motion will be described.

[0081] The determination value P.sub.t(t) of FIRE applied to the translational motion is defined as follows.

P.sub.t(t)=F(t).Math.V(t)   (12)

[0082] This expression is obtained by replacing the force f acting on the particles in the determination value (Expression (1)) of FIRE applied to the individual particles with the translational force F acting on the analysis object 10, and replacing the velocity v of the particles with the translational velocity V of the analysis object 10.

[0083] In a case where FIRE is applied to the translational motion of the analysis object 10, when the determination value P.sub.t(t) is zero or negative, the translational velocity V(t+Δt) of the analysis object 10 in the next state may be determined as follows.

V(t+Δt)=0, (P(t)≤0)   (13)

[0084] Expression (13) corresponds to Expression (2) in which the velocity v of the particles is set to zero.

[0085] When the determination value P.sub.t(t) is positive, the translational velocity V(t+Δt) in the next state may be obtained by solving the following equation of motion for the analysis object 10.

[00009]$\begin{array}{cc}\hat{V}\left(t\right)=\frac{F\left(t\right)}{M}-\alpha .Math.V\left(t\right).Math..Math.\left\{\hat{V}\left(t\right)-\hat{F}\left(t\right)\right\},\left(P_t\left(t\right)>0\right)& \left(14\right)\end{array}$

[0086] Expression (14) corresponds to the equation of motion for particles (Expression (3)). Here, α is an attenuation constant.

[0087] The determination value P.sub.r(t) of FIRE applied to the rotational motion is defined as follows.

P.sub.r(t)=T(t).Math.ω(t)   (15)

[0088] Expression (15) is obtained by replacing the translational force F and the translational velocity V of Expression (12) with respect to the translational motion with the torque T and the angular velocity ω, respectively.

[0089] Ina case where FIRE is applied to the rotational motion of the analysis object 10, when the determination value P.sub.r(t) is zero or negative, the angular velocity ω(t+Δt) of the analysis object 10 in the next state may be determined as follows.

ω(t+Δt)=0, (P.sub.r(t)≤0)   (16)

[0090] Expression (16) corresponds to Expression (13) in which the translational velocity V of the analysis object 10 is set to zero.

[0091] When the determination value P.sub.r(t) is positive, the angular velocity ω(t+Δt) in the next state may be obtained by solving the following equation of motion for the analysis object 10.

[00010]$\begin{array}{cc}\stackrel{.}{\omega }\left(t\right)=\frac{T\left(t\right)}{I}-\beta .Math.\omega \left(t\right).Math..Math.\left\{\hat{T}\left(t\right)-\hat{\omega }\left(t\right)\right\},\left(P_r\left(t\right)>0\right)& \left(17\right)\end{array}$

[0092] Expression (17) corresponds to the equation of motion for translational motion (Expression (14)). Here, β is an attenuation constant. I is the moment of inertia of the analysis object 10 and is defined by Expression (11). Here, T(t)/I means a vector (Tx/Ix, Ty/Iy, Tz/Iz).

[0093] From Expression (14), the translational force F.sub.add to be additionally applied to the analysis object 10 by FIRE is described by the following expression.

F.sub.add=−αM|V(t)|.Math.{{circumflex over (V)}(t)−{circumflex over (F)}(t)}  (18)

[0094] From Expression (17), the x component of the torque T.sub.add to be additionally applied to the analysis object 10 by FIRE is described by the following expression. The same applies to the y component and the z component of the torque T.sub.add.

T.sub.add,x=−βI.sub.x|ω(t)|.Math.[{circumflex over (T)}(t)−{circumflex over (ω)}(t)].sub.x   (19)

[0095] After step S4, the processing unit 41 distributes the translational force F.sub.add and the torque T.sub.add to be additionally applied to the analysis object 10 to each particle 11, and the force to be additionally applied to the particles 11 is calculated (step S5).

[0096] The translational force F.sub.add to be additionally applied is distributed to each particle 11 based on the following expression.

[00011]$\begin{array}{cc}{f}_i=\frac{F}{N}& \left(20\right)\end{array}$

[0097] The x component t.sub.i,x, the y component t.sub.i,y, and the z component t.sub.i,z of the force to be applied to the i-th particle 11 when the torque T.sub.add to be additionally applied is distributed to each particle 11 are described by the following expression.

[00012]$\begin{array}{cc}{t}_{i,x}=\frac{m_i\left(y^2+z^2\right)}{I_x}{T}_{\mathrm{add},x}{t}_{i,y}=\frac{m_i\left(z^2+x^2\right)}{I_x}{T}_{\mathrm{add},y}{t}_{i,x}=\frac{m_i\left(x^2+y^2\right)}{I_x}{T}_{\mathrm{add},z}& \left(21\right)\end{array}$

[0098] When the determination value P.sub.t(t) for the translational motion of the analysis object 10 is zero or negative, the translational velocity V defined by Expression (8) becomes zero, by subtracting the translational velocity V of the analysis object 10 at the present time from the velocity v.sub.i of each of the particles 11. When the determination value P.sub.r(t) for the rotational motion is zero or negative, ω.sub.x, ω.sub.y, and ω.sub.z of Expression (10) become zero, by subtracting the velocity r.sub.iω of each of the particles 11 corresponding to the rotational velocity of the analysis object 10, from the velocity v.sub.i of each of the particles 11.

[0099] The procedure from step S2 to step S5 is repeated until the analysis end condition is satisfied (step S6). When solving the equation of motion in step S2, the force to be additionally applied, obtained in step S5, is added to each of the particles 11.

[0100] Next, the excellent effects of the above embodiment will be described.

[0101] In the above embodiment, in step S4 (FIG. 3), since FIRE is applied to the motion of the analysis object 10 regarded as a rigid body, the calculation time until the analysis object 10 as a whole reaches an equilibrium state can be shortened.

[0102] In order to confirm the excellent effect of the above embodiment, a static analysis of a cylindrical rotating body was performed using the simulation method according to the above embodiment. Next, the simulation results will be described with reference to FIGS. 4 and 5. In the simulation, a part of the region of one end surface of the rotating body is attached to the fixing member and fixed, and a constant torque is applied to the outer peripheral surface of the rotating body.

[0103] FIG. 4 is a graph showing a temporal change of a torque applied to a fixing member that fixes a part of a region of one end surface of a rotating body. The horizontal axis represents the elapsed time from the time when a certain time has elapsed from the start of analysis, and the vertical axis represents the torque normalized with the applied torque. The thick solid line in FIG. 4 indicates the torque obtained by using the simulation method according to the above embodiment, the thin solid line indicates the torque obtained by using the simulation method according to the reference example using Expressions (4) to (6), and the broken line indicates the torque obtained by simulation without applying FIRE. In the equilibrium state, the torque applied to the fixing member is equal to the applied torque. That is, the normalized torque becomes 1.

[0104] When FIRE is not applied, the calculated value of torque oscillates around an equilibrium value. The calculated value of torque according to the reference example gradually approaches the equilibrium value after overshooting. On the other hand, the calculated value of torque according to the above embodiment is a substantially equilibrium value, the calculated value of torque gradually increases, and the calculated value of torque reaches a substantially equilibrium state in a short time.

[0105] FIG. 5 is a graph showing a temporal change of an angular velocity co when the rotating body is regarded as a rigid body. The horizontal axis represents the elapsed time from the time when a certain time has elapsed from the start of analysis, and the vertical axis represents the angular velocity normalized with the maximum value set to 1. The thick solid line in FIG. 4 indicates the angular velocity obtained by using the simulation method according to the above embodiment, the thin solid line indicates the angular velocity obtained by using the simulation method according to the reference example using Expressions (4) to (6), and the broken line indicates the angular velocity obtained by simulation without applying FIRE. In the equilibrium state, the angular velocity becomes zero.

[0106] When FIRE is not applied, the calculated value of the angular velocity oscillates around zero. The calculated value of the angular velocity according to the reference example decreases after reaching the maximum, and becomes zero discontinuously immediately before the equilibrium state. When the calculated value of the angular velocity becomes zero discontinuously, the determination value P(t) of Expression (4) is zero or negative. On the other hand, in the calculated value of the angular velocity according to the above embodiment, the calculated value of the angular velocity becomes zero discontinuously immediately after the calculated value of the angular velocity becomes a substantially maximum value. When the calculated value of the angular velocity becomes zero discontinuously, the determination value P.sub.t(t) of Expression (12) or the determination value P.sub.r(t) of Expression (15) is zero or negative. In the present simulation, the determination value P.sub.r(t) is zero or negative.

[0107] As shown in FIGS. 4 and 5, by using the simulation method according to the above embodiment, it can be seen that the time required to reach the equilibrium state is shortened, as compared to the case where FIRE is not applied and the case where the simulation method according to the reference example is used.

[0108] Next, a modification example of the above embodiment will be described.

[0109] In the above embodiment, in step S3 (FIG. 3), the translational velocity and the angular velocity are obtained as the velocity of the analysis object 10 regarded as a rigid body, but only one of the translational velocity and the angular velocity may be obtained, and FIRE may be applied to only one of the translational motion and the rotational motion of the analysis object 10.

[0110] The above-described embodiment is an example, and the present invention is not limited to the above-described embodiment. For example, it will be apparent to those skilled in the art that various modifications, improvements, combinations, and the like can be made.

[0111] It should be understood that the invention is not limited to the above-described embodiment, but may be modified into various forms on the basis of the spirit of the invention. Additionally, the modifications are included in the scope of the invention.