JP7248487B2 - FEM analysis method, system and program for structures - Google Patents

Description

The present invention relates to an FEM analysis method, system and program for a structure.

A finite element method (FEM) is used for dynamic analysis of structures. In FEM model data used for FEM, a structure to be analyzed is expressed by dividing it into a plurality of elements, and a closed region connecting a plurality of nodes is expressed as one element. In FEM, the deformation of a structure is calculated by solving an equation of motion using the physical property values (Young's modulus, Poisson's ratio, mass, etc.) of the structure and the acting force for each node. Non-Patent Document 1 is available as a document describing FEM.

Although it is not a structure simulation, the molecular dynamics solver (software) for executing molecular simulation has a parallel computing function for solving the equation of motion of a multi-particle system with high efficiency, and LAMMPS, gromacs, etc. Freely available software is known. Non-Patent Document 2 is a document related to molecular dynamics.

Hutton, D. V. (2017). Fundamentals of finite element analysis, Chapter 10, Structural Dynamics. McGraw-Hill. Frenkel, D., Smit, B., Tobochnik, J., McKay, S. R., & Christian, W. (1997). Understanding molecular simulation. Computers in Physics, 11(4), 351-354.

In order to improve analysis accuracy in FEM analysis, a model with smaller element sizes may be used. Since most FEM analysis software performs matrix calculations that treat all elements as one, there are many cases where there is no parallel calculation function. .

SUMMARY OF THE INVENTION An object of the present invention is to provide an FEM analysis method, system and program for a structure capable of performing parallel calculations at a reasonable cost.

The FEM analysis method of the structure of the present invention comprises:
A method, performed by one or more processors, comprising:
(a) Acquiring FEM model data representing a structure with a plurality of elements, wherein the FEM model data is set with a plurality of nodes and physical property values that partition each element;
(b) obtaining a lumped mass matrix expressing the mass of the structure distributed to each node and the initial position and initial velocity of each node based on the FEM model data;
(c) calculating an element stiffness matrix for each element based on the initial position of each node and the physical property value;
(d) calculating, element by element, a displacement matrix representing the displacement of the position at time (t) with respect to the initial position of each node;
(e) calculating, for each element, a force acting on each node based on the element stiffness matrix and the displacement matrix;
(f) based on the position and velocity of each node at time (t), the force acting on each node, and the lumped mass matrix, the position and position of each node at the next time (t+unit time) and calculating the velocity;
including
The processes (d), (e) and (f) are repeatedly executed with the time (t) passing by a unit time until the target time is reached from the initial time.

In this way, the node in the FEM model is treated as one particle, and the force acting on the node is calculated based on the element stiffness matrix based on the FEM model and the positional displacement with respect to the initial position. ), the position of each node at the next time (t+unit time) can be calculated from the force and the position of the node, and FEM calculation can be executed using a molecular dynamics solver. In addition, molecular dynamics solvers, including software that can be used free of charge, have parallel computing functions, so FEM analysis of structures can be realized by parallel computing at a reasonable cost.

Figure showing the FEM analysis system of the structure of the present invention Flowchart of processing performed by the system Diagram showing FEM model A diagram showing the initial position of each node of the FEM model of the embodiment. A diagram showing the position of each node of the FEM model when it collides with the wall of the embodiment.

An embodiment of the present invention will be described below with reference to the drawings.

[FEM analysis system for structures]
A system 1 of this embodiment is a system for analyzing a structure using the finite element method.

As shown in FIG. 1, the system 1 includes an FEM model acquisition unit 10, a lumped mass matrix acquisition unit 11, an initial position/velocity acquisition unit 12, an element stiffness matrix calculation unit 13, a displacement matrix calculation unit 14, a force It has a calculator 15 and a molecular dynamics solver 16 . These units 10 to 15 are implemented by cooperation of software and hardware when the processor executes the processing routine shown in FIG. In this embodiment, the processor in one device executes the processing of each unit, but the present invention is not limited to this. For example, it may be distributed using a network so that a plurality of processors execute the processing of each unit. That is, one or more processors perform processing.

Generally, FEM analysis uses an FEM solver (software), transfers various data such as FEM model data, loads applied to the FEM model, and boundary conditions to the FEM solver to operate the FEM solver. FEM analyzes a structure by dividing it into multiple elements, which have multiple nodes that make up the elements. The inventor recognized that it is calculating the force to act. In other words, we thought that a molecular dynamics solver could be used by implementing FEM-like force calculation.

Therefore, in the FEM analysis method of this embodiment, instead of using the FEM solver, one node constituting the element of the FEM model is treated as one particle, and a molecular simulation of a material composed of a plurality of particles is performed. It uses a molecular dynamics solver that Molecular dynamics solvers include LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), gromacs, etc. as free software, and LAMMPS is used in this embodiment. Data to be passed to the molecular dynamics solver include initial positions and initial velocities of a plurality of particles that make up the material, and interactions (forces and potentials) acting on each particle. The interaction is determined according to the interparticle distance. The molecular dynamics solver has a potential calculator and a particle position updater. As an overview of the operation flow of a general molecular dynamics solver, the potential calculation unit calculates the force (potential) acting on each particle at time (t) based on the position and velocity of each particle at time (t). Calculated according to distance. Next, the particle position update unit calculates the position and velocity of each particle at the next time (t+1) based on the position and velocity of each particle at time (t) and the force (potential) acting on each particle. In this way, the calculation of the force (potential) acting on the particles, the position and velocity of each particle at time (t), and the position and velocity of each particle at the next time (t+1) based on the force acting on each particle Update and are executed repeatedly. This iterative process is performed until the target point in time. In this embodiment, since the potential calculation unit of the molecular dynamics solver 16 cannot be used as it is in the FEM analysis, the force calculation unit 15 is provided instead. The particle position updater of the molecular dynamics solver 16 is used as is.
A specific description will be given below.

The FEM model acquisition unit 10 shown in FIG. 1 acquires FEM model data D1 representing a structure with a plurality of elements. The FEM model data D1 indicates the three-dimensional shape of the structure by a plurality of nodes that partition each element, and physical property values (Young's modulus E, Poisson's ratio v, mass density ρ) are set for the elements. The illustration in FIG. 3 is a tetrahedron in which one element is represented by four nodes, and two elements indicated by elements e=1 and 2 are illustrated. Nodes i=1 to 4 forming an element e=1, and nodes i=2 to 5 forming an element e=2. Each node i has xyz coordinate values and information indicating to which node it is connected and to which element it belongs. The FEM model data D1 is frequently used in general finite element method simulations. In addition, in the present embodiment shown in FIG. 3, the element of the FEM model data D1 is a tetrahedron, but it is not limited to this. In principle any polyhedral element can be used.

ただし、ρは構造物の質量密度を示し、Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示す。 The lumped mass matrix acquisition unit 11 shown in FIG. 1 acquires a lumped mass matrix M _e that indicates the mass of the structure distributed to each node. The lumped mass matrix acquisition unit 11 acquires the lumped mass matrix M _e as it is when the lumped mass matrix M _e is defined in the FEM model data D1, and calculates it when the lumped mass matrix M _e is not defined. . The lumped mass matrix M _e of the element e is expressed by the following equation (6). The lumped mass matrix acquisition unit 11 acquires or calculates the lumped mass matrix M _e for all elements. In the example shown in FIG. 3, since the number of elements is 5070, 5070 lumped mass matrices M _e are obtained. The lumped mass matrix M _e has masses only on the diagonal. A standard FEM analysis may use a consistent mass matrix with off-diagonal elements, but we use a lumped mass matrix to avoid drastic changes to the molecular dynamics solver.

where ρ indicates the mass density of the structure, V _e indicates the volume of the element e, and L _e indicates the number of nodes in the element e.

The initial position/velocity acquisition unit 12 shown in FIG. 1 acquires the initial position and initial velocity of each node based on the FEM model data D1, and stores these values in the working memory D2. The initial position of each node is the coordinates (xyz coordinates) of each node in the FEM model data D1. The initial velocity at each node is zero. In the example shown in FIG. 3, the number of elements is 5070 and the total number of nodes is 1204, so the coordinates of 1204 nodes are acquired.

Ｖ_ｅは要素ｅの体積を示し、Ｌ_ｅは要素ｅにおける節点数を示し、Ｎ^ｅ _ｉ（ｒ）は要素ｅの節点ｉ（ｉ＝１，…，Ｌ_ｅ）の形状関数を示し、Ｅは構造物のヤング率を示し、ｖは構造物のポアソン比を示す。 The element stiffness matrix calculator 13 shown in FIG. 1 calculates the element stiffness matrix _Ke for each element based on the initial position and physical property value of each node. Element stiffness matrix K _e of element e is expressed by equations (1) to (3). The element stiffness matrix calculator 13 calculates the element stiffness matrix _Ke for all the elements. In the embodiment shown in FIG. 3, since the number of elements is 5070, 5070 element stiffness matrices K _e are calculated.

V _e denotes the volume of element e, L _e denotes the number of nodes in element e, N ^e _i (r) denotes the shape function of node i (i=1, . . . , L _e ) of element e, and E denotes the Young's modulus of the structure and v denotes the Poisson's ratio of the structure.

ただし、Ｓ^ｅ _ｉは要素ｅの節点ｉにおけるｘ成分の変位を示し、Ｕ^ｅ _ｉは要素ｅの節点ｉにおけるｙ成分の変位を示し、Ｗ^ｅ _ｉは要素ｅの節点ｉにおけるｚ成分の変位を示す。 The displacement matrix calculator 14 shown in FIG. 1 calculates a displacement matrix d _e indicating the displacement of the position at time (t) with respect to the initial position of each node for each element. The displacement matrix d _e is expressed by Equation (5). The displacement matrix calculator 14 calculates a displacement matrix d _e indicating displacements of all nodes forming each element for all elements. In the example shown in FIG. 3, since the number of elements is 5070, 5070 displacement matrices d _e are calculated.

where S ^e _i denotes the x-component displacement at node i of element e, U ^e _i denotes the y-component displacement at node i of element e, and W ^e _i denotes the z-component displacement at node i of element e. indicates

ただし、Ｆ^ｅ _ｉ，ｘは要素ｅの節点ｉに作用する力のｘ成分を示し、Ｆ^ｅ _ｉ，ｙは要素ｅの節点ｉに作用する力のｙ成分を示し、Ｆ^ｅ _ｉ，ｚは要素ｅの節点ｉに作用する力のｚ成分を示す。 The force calculator 15 shown in FIG _. 1 calculates the node points _i (i = 1 _{, ...} , _{L e} ) is ^calculated for ^each element ^. A force F _e acting on each node i constituting the element e is expressed by Equation (4).

where F ^e _i,x denotes the x component of the force acting on the node i of the element e, F ^e _i,y denotes the y component of the force acting on the node i of the element e, and F ^e _i,z The z component of the force acting on node i of element e is shown.

The molecular dynamics solver 16 shown in FIG. Based on _Me , the position and velocity of each node at the next time (t+unit time) are calculated. The calculated position and velocity of each node are stored in the working memory D2.

[FEM analysis method for structures]
A structure FEM analysis method executed by one or more processors in the system 1 shown in FIG. 1 will be described with reference to FIG.

First, in step ST101, the FEM model acquisition unit 10 acquires FEM model data D1 representing a structure with a plurality of elements, and a plurality of nodes and physical property values (Young's modulus E, Poisson's ratio v, FEM model data D1 in which mass density ρ, etc.) is set is obtained.

In the next step ST102, the lumped mass matrix acquisition unit 11 acquires the lumped mass matrix M _e representing the mass of the structure distributed to each node based on the FEM model data D1 for each element, and the initial position/velocity acquisition unit 12 obtains the initial position and initial velocity of each node based on the FEM model data D1.

Steps ST103, ST104, and ST105 are executed for each element (e=1 to 5070) for all elements (5070 elements in the example of FIG. 3). When the processes of steps ST103, ST104, and ST105 are executed for all elements, the process proceeds to step ST106.

In step ST103, the element stiffness matrix calculator 13 calculates the element stiffness matrix based on the initial positions and physical property values of each node of the element e. Element stiffness matrix K _e of element e is expressed by equations (1) to (3).

In step ST104, the displacement matrix calculator 14 calculates the displacement matrix d _e of the element e at time (t). The displacement matrix d _e indicates the displacement of the current position with respect to the initial position of each node i (i=1 to L _e ) that constitutes the element e. The displacement matrix d _e is expressed by Equation (5).

In step ST105, the force calculator 15 calculates the force acting on each node of the element e based on the element stiffness matrix K _e and the displacement matrix d _e of the element e. Forces (F ^e _i,x , F ^e _i,y , F ^e _i,z ) acting on _each node i (i=1, . )

After executing the processes of steps ST103, ST104, and ST105 for all elements, in step ST106, the molecular dynamics solver 16 calculates the position and velocity of each node at time (t) and the force and concentration acting on each node. The position and velocity of each node at the next time (t+unit time) are calculated based on the mass matrix _Me .

In the next step ST107, it is determined whether or not the processing of ST103 to ST106 has been repeatedly executed from the initial time until reaching the target time with the time (t) having elapsed by the unit time. If the target time has not been reached (ST107: NO), the time (t) is allowed to elapse by a unit time, and the process returns to ST103. If the target time has been reached (ST107: YES), the FEM analysis process ends.

In this embodiment, the calculation of the element stiffness matrix K _e in step ST103 is executed each time a unit time elapses. This is because the element stiffness matrix K _e has 5070 elements, so memory consumption cannot be ignored when stored in the memory, and the calculation cost is low even if the calculation is repeated every unit time. Of course, the element stiffness matrix _Ke for all elements may be stored in the memory at the time of calculation, and may not be calculated each time a unit time elapses.

The results of the simulation method of the present invention will be explained.
4 and 5 show the results of simulating the falling and repulsive motion of an elastic sphere with a radius of 10 mm. FIG. 4 shows the position of each node of the FEM model in the initial position, and FIG. 5 shows the position of each node of the FEM model at the time of collision with the wall. The physical properties are Young's modulus E of 1 MPa, Poisson's ratio v of 0.49, and mass density ρ of 970 kg/m ³ , and the structure is equivalent to a rubber ball. The element type is a tetrahedral primary element with four nodes per element, and there are 5070 elements and 1204 nodes in total. The height of the center of the elastic sphere is set to 20 mm, and the ball is allowed to fall naturally by gravity, and a rigid wall is set at a height of 0 mm. The number of time steps is 30000 (corresponding to about 1.0 second).

In order to show that the parallel computing function has been realized, the number of processes is set to 1, 2, 4, 8, and 16, and the computing time is shown below. As shown below, the calculation time is shortened according to the number of processes, so it can be seen that the parallel calculation function can be reasonably implemented in the FEM analysis.
Number of processes: 1 Calculation time (seconds): 1106.8
Number of processes: 2 Computation time (seconds): 560.4
Number of processes: 4 Computation time (seconds): 282.6
Number of processes: 8 Computation time (seconds): 150.2
Number of processes: 16 Computation time (seconds): 83.0

As described above, the FEM analysis method for the structure of the present embodiment is
A method, performed by one or more processors, comprising:
(a) FEM model data D1 representing a structure with a plurality of elements e, in which a plurality of nodes i and physical property values (Young's modulus E, Poisson's ratio v, mass density ρ) that partition each element e are set obtaining FEM model data D1 (ST101);
(b) Acquiring a lumped mass matrix M _e representing the mass of the structure distributed to each node and the initial position and initial velocity of each node based on the FEM model data D1 (ST102);
(c) calculating an element stiffness matrix K _e for each element based on the initial position and physical property value of each node (ST103);
(d) calculating, element by element, a displacement matrix _de representing the displacement of the position at time (t) with respect to the initial position of each node (ST104);
(e) calculating the force F _e acting on each node based on the element stiffness matrix K _e and the displacement matrix d _e for each element (ST105);
(f) based on the position and velocity of each node at time (t), the force acting on each node, and the lumped mass matrix, the position and position of each node at the next time (t+unit time) and calculating the speed (ST106);
including
The processes of (d), (e) and (f) are repeatedly executed with the unit time elapsed from the initial time point until the target time point is reached.

The FEM analysis system for the structure of this embodiment is
FEM model data D1 representing a structure with a plurality of elements e, in which a plurality of nodes i partitioning each element e and physical property values (Young's modulus E, Poisson's ratio v, mass density ρ) are set. an FEM model acquisition unit 10 that acquires data D1;
a lumped mass matrix acquisition unit 11 that acquires a lumped mass matrix M _e that expresses the mass of the structure by distributing it to each node based on the FEM model data D1;
an initial position/velocity acquisition unit 12 that acquires the initial position and initial velocity of each node based on the FEM model data D1;
an element stiffness matrix calculator 13 for calculating an element stiffness matrix K _e for each element based on the initial position and physical property value of each node;
a displacement matrix calculator 14 for calculating, for each element, a displacement matrix _de representing the displacement of the position at time (t) with respect to the initial position of each node;
a force calculation unit 15 for calculating, for each element, a force F _e acting on each node based on the element stiffness matrix K _e and the displacement matrix d _e ;
A molecular dynamics solver that calculates the position and velocity of each node at the next time (t + unit time) based on the position and velocity of each node at time (t), the force acting on each node, and the lumped mass matrix. 16 and
with
The process by the displacement matrix calculator 14, the force calculator 15, and the molecular dynamics solver 16 is repeatedly executed from the initial time point until the target time point is reached, with the time point (t) passing by a unit time. .

In this way, the node in the FEM model is treated as one particle, and the force acting on the node is calculated based on the element stiffness matrix based on the FEM model and the positional displacement with respect to the initial position. ), the position of each node at the next time (t+unit time) can be calculated from the force and the position of the node, and FEM calculation can be executed using a molecular dynamics solver. In addition, molecular dynamics solvers, including software that can be used free of charge, have parallel computing functions, so FEM analysis of structures can be realized by parallel computing at a reasonable cost.

As in the present embodiment, the calculation of the element stiffness matrix _Ke by the element stiffness matrix calculator 13 in (c) is preferably performed each time a unit time elapses.

In this way, even if the element stiffness matrix K _e is calculated each time the unit time elapses, the calculation cost is not large, and the memory consumption caused by storing the element stiffness matrix K _e of all the elements in the memory can be reduced. can be suppressed.

As in this embodiment, the element stiffness matrix K _e of the element e is preferably represented by equations (1) to (3). A preferred embodiment.

As in this embodiment, the forces ( ^Fei _,x, Fei ^, _y , ^Fei _{,z )} acting on each node i (i=1, . . . , L e ) constituting the element e are Expressions (4) and (5) are preferred. A preferred embodiment.

A program according to the present embodiment is a program that causes a computer to execute the above method. By executing this program, it is possible to obtain the effects of the above method.

Although the embodiments of the present invention have been described above with reference to the drawings, it should be considered that the specific configuration is not limited to these embodiments. The scope of the present invention is indicated not only by the description of the above embodiments but also by the scope of claims, and includes all modifications within the scope and meaning equivalent to the scope of claims.

For example, each unit 10 to 15 shown in FIG. 1 is realized by executing a predetermined program by a processor of a computer, but each unit may be configured by a dedicated circuit. In addition, in the present embodiment, the processors in one computer implement the units 10 to 15, but they may be installed in at least one or a plurality of processors in a distributed manner.

It is possible to adopt the structure adopted in each of the above embodiments in any other embodiment. The specific configuration of each part is not limited to the above-described embodiment, and various modifications are possible without departing from the scope of the present invention.

REFERENCE SIGNS LIST 10 FEM model acquisition unit 11 lumped mass matrix acquisition unit 12 initial position/velocity acquisition unit 13 element stiffness matrix calculation unit 14 displacement matrix calculation unit 15 force calculation unit 16 molecular dynamics solver

Claims (9)

A method, performed by one or more processors, comprising:
(a) Acquiring FEM model data representing a structure with a plurality of elements, wherein the FEM model data is set with a plurality of nodes and physical property values that partition each element;
(b) obtaining a lumped mass matrix expressing the mass of the structure distributed to each node and the initial position and initial velocity of each node based on the FEM model data;
(c) calculating an element stiffness matrix for each element based on the initial position of each node and the physical property value;
(d) calculating, element by element, a displacement matrix representing the displacement of the position at time (t) with respect to the initial position of each node;
(e) calculating, for each element, a force acting on each node based on the element stiffness matrix and the displacement matrix;
(f) based on the position and velocity of each node at time (t), the force acting on each node, and the lumped mass matrix, the position and position of each node at the next time (t+unit time) and calculating the velocity;
including
An FEM analysis method for a structure, wherein the processes of (d), (e) and (f) are repeatedly executed while the time (t) elapses as a unit time until the target time is reached from the initial time. 2. The method according to claim 1, wherein the calculation of said element stiffness matrix in said (c) is calculated each time said unit time elapses. The method according to claim 1 or 2, wherein the element stiffness matrix K _e of the element e is represented by equations (1)-(3).

V _e denotes the volume of element e, L _e denotes the number of nodes in element e, N ^e _i (r) denotes the shape function of node i (i=1, . . . , L _e ) of element e, and E denotes the Young's modulus of the structure and v denotes the Poisson's ratio of the structure. Forces (F ^e _i,x , F ^e _i,y , F ^e _i,z ) acting on _each node i (i=1, . ), the method of claim 3 .

where S ^e _i denotes the x-component displacement at node i of element e, U ^e _i denotes the y-component displacement at node i of element e, and W ^e _i denotes the z-component displacement at node i of element e. , F ^e _i,x is the x component of the force acting on the node i of the element e, F ^e _i,y is the y component of the force acting on the node i of the element e, and F ^e _i,z indicates the z component of the force acting on the node i of the element e. an FEM model acquisition unit that acquires FEM model data representing a structure with a plurality of elements and in which a plurality of nodes that partition each element and physical property values are set;
a lumped mass matrix acquisition unit that acquires a lumped mass matrix representing the mass of the structure distributed to each node based on the FEM model data;
an initial position/velocity acquisition unit that acquires the initial position and initial velocity of each node based on the FEM model data;
an element stiffness matrix calculation unit that calculates an element stiffness matrix for each element based on the initial position of each node and the physical property value;
a displacement matrix calculation unit that calculates, for each element, a displacement matrix representing the displacement of the position at time (t) with respect to the initial position of each node;
a force calculation unit that calculates, for each element, a force acting on each node based on the element stiffness matrix and the displacement matrix;
A molecular dynamics solver that calculates the position and velocity of each node at the next time (t + unit time) based on the position and velocity of each node at time (t), the force acting on each node, and the lumped mass matrix. and,
with
The process by the displacement matrix calculator, the force calculator, and the molecular dynamics solver is repeatedly executed from the initial time point until the target time point is reached, with the time point (t) passing by a unit time. FEM analysis system for structures. 6. The system according to claim 5, wherein the calculation of the element stiffness matrix by the element stiffness matrix calculator is performed each time the unit time elapses. The system according to claim 5 or 6, wherein the element stiffness matrix K _e of the element e is expressed by equations (1)-(3).

V _e denotes the volume of element e, L _e denotes the number of nodes in element e, N ^e _i (r) denotes the shape function of node i (i=1, . . . , L _e ) of element e, and E denotes the Young's modulus of the structure and v denotes the Poisson's ratio of the structure. Forces (F ^e _i,x , F ^e _i,y , F ^e _i,z ) acting on _each node i (i=1, . 8. The system of claim 7, expressed as .

where S ^e _i denotes the x-component displacement at node i of element e, U ^e _i denotes the y-component displacement at node i of element e, and W ^e _i denotes the z-component displacement at node i of element e. , F ^e _i,x is the x component of the force acting on the node i of the element e, F ^e _i,y is the y component of the force acting on the node i of the element e, and F ^e _i,z indicates the z component of the force acting on the node i of the element e. A program that causes one or more processors to execute the method according to any one of claims 1 to 4.

Images