CN102184298B - Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming - Google Patents

Abstract

The invention relates to a method for storing and generating stiffness matrix of finite element analysis in metal volume plastic forming, which can significantly save the storage space of stiffness matrix and improve the efficiency of solving stiffness equation. It comprises the following steps: (1.1) according to the finite element network of volume plastic forming workpiece The grid model is used to determine the "generalized adjacent node pairs" of grid nodes and their correspondence with the non-zero submatrix in the overall stiffness matrix of the grid model; (1.2) Based on the "generalized adjacent node pairs", determine the overall grid model The position and quantity of the non-zero sub-matrix in the stiffness matrix, store the non-zero sub-matrix; (1.3) Determine the non-zero sub-matrix where each element stiffness matrix element of the grid model is located in the overall stiffness matrix and its storage array in the overall stiffness matrix The position in the grid model is assembled into the overall stiffness matrix storage array of the grid model; (1.4) The overall stiffness equation of the grid model is solved by using the symmetric positive definite coefficient matrix relaxation preconditioning co-value gradient method.

Description

Storage and Generation Method of Stiffness Matrix for Finite Element Analysis in Metal Volume Plastic Forming

technical field

The invention relates to a finite element analysis stiffness matrix storage and generation method in metal volume plastic forming, which is suitable for finite element numerical simulation of volume forming processes such as forging, extrusion, rolling, and rotary forging.

Background technique

In industrial production, forging, extrusion, rolling and other metal volume plastic forming manufacturing technologies are one of the main technologies supporting national economic development and national defense construction. use. The volume plastic forming process is to plastically deform the metal material (blank) at room temperature or high temperature (cold, warm, hot forming) under the action of the forging equipment and the mold installed on it, and form it into parts of the required shape and performance The bulk plastic forming material is generally assumed to be a rigid (viscosity) plastic material. Therefore, the blank undergoes large and severe plastic deformation in the mold cavity, and there are material nonlinearity, geometric nonlinearity and complex boundary conditions. The forming process Solving is very complicated. Forming process and mold design have an important impact on the quality of formed parts and whether they can be formed smoothly, and are also the key to energy saving and precise forming. Mastering the plastic deformation behavior and flow law of metal materials in the mold cavity is the key theoretical basis for designing the forming process and mold. As an important method for the numerical analysis of engineering problems, the finite element method has played an important role in the engineering field. The finite element numerical simulation method of the metal volume plastic forming process is the main method for studying the plastic deformation behavior and flow law of metal materials. It can obtain the flow of metal materials. Regularity and distribution results of various field variables, verify the feasibility of process design, parameter selection and mold design scheme, and then realize the optimization of forming process and mold, improve the speed and quality of product development and reduce development costs.

For the finite element analysis of the metal volume plastic forming problem, firstly, the mathematical model of the metal volume plastic deformation mechanics problem is established through the plastic deformation control equation, and the geometric model (mold and denatured material) is discretized by the finite element method, and the finite element stiffness equation is established. To realize the numerical simulation of the forming process, numerical analysis methods such as finite element have played an important role in the field of engineering, especially in the fields of metal forming manufacturing and structural analysis.

The steps of finite element numerical solution to the metal volume plastic forming problem are as follows: establish the CAD geometric model of the problem (geometric shape of the deformed material), discretize the CAD geometric model with finite elements, and obtain the finite element mesh model of the geometric model, the mesh model It is composed of a finite number of grid units; firstly, the stiffness equation of each unit in the grid model is established, and the stiffness equations of all units are assembled to form a nonlinear overall stiffness equation system for all node velocities of the grid model; by The nonlinear stiffness equation of node velocity is linearized to obtain the linear stiffness equation system of node velocity increment; the Newton-Raphson iterative method is used to iteratively solve the velocity field. Since the node velocity increment is solved, for a metal volume plastic forming process, the forming time needs to be decomposed into hundreds or even hundreds of time steps (long), and the element stiffness matrix needs to be calculated in each time step and nodal force vectors, assembled into the overall stiffness equation of the nodal velocity increment, applying velocity boundary conditions, solving the overall stiffness equation of the nodal velocity increment (large linear equations), obtaining the velocity correction, and then using the velocity correction to the velocity field Make corrections, and then repeat the above steps to iteratively solve the velocity field until the velocity field converges. Generally, at least 7 to 10 iterations are required to obtain the converged velocity field; according to the obtained convergent velocity field, refresh the shape of the deformed workpiece and die position, the numerical calculation within one time step can be completed. Then carry out the numerical calculation in the next step until the required degree of deformation is reached.

It can be seen that for the numerical simulation of the metal volume plastic forming process, it is necessary to frequently solve the overall stiffness equation (large linear equation system). For example, for medium and complex three-dimensional metal volume plastic forming problems, more than 200 simulation steps are generally required, and the velocity field needs to be iterated 9 times in each simulation step, and each iteration of the velocity field needs to solve a large linear equation set. If the number of finite element nodes in a metal volume plastic forming problem is 50,000 (a moderately complex three-dimensional problem), it is necessary to solve a large linear square system consisting of 150,000 equations 1800 times. If more complex problems are analyzed, the computational workload will be even greater.

In recent years, with the needs of large-scale projects such as major equipment, ships, automobiles and wind power, the demand for large and heavy metal volume plastic forming parts is increasing, and the forming problems that need to be analyzed are becoming more and more complex. The number of units and units is large, the amount of data is huge, the efficiency of finite element calculation is low, and the required computer memory increases sharply. Even when analyzing large and super-large problems, the problem of insufficient storage space often occurs, which greatly limits the use of finite element numerical methods in metal Application in forming process analysis.

Not only the problem of metal volume plastic forming, but also the finite element analysis of many other engineering problems will eventually come down to solving large linear equations. The efficiency of finite element solution mainly depends on the storage and generation of the stiffness matrix of the linear equations and the efficiency of solving the stiffness equation. .

The finite element stiffness matrix of the metal volume plastic forming problem is a large, strip-shaped sparse matrix, and the non-zero elements of the matrix are generally concentrated in the diagonal area. In view of this situation and in order to improve the calculation efficiency, many storage methods for saving storage space have been proposed , for example, half-bandwidth storage method, variable bandwidth storage method, sparse equations preprocessing iterative solution method, etc. Although these methods reduce the storage space of the stiffness matrix, they inevitably store some non-zero elements and waste some At the same time, these methods all depend on the node number, and the order of the node number needs to be optimized, which also wastes calculation time. The existing one-dimensional variable bandwidth compression storage method uses a one-dimensional array to store the non-zero elements of the stiffness matrix in the order of rows, and uses two auxiliary arrays to record the column number and the position of the starting element of each column in the stiffness matrix array, saving storage space, avoiding the problem of node number optimization, but the auxiliary array used has the same length as the stiffness matrix storage array, which still consumes storage space and wastes computing time.

Therefore, when finite element analysis is used to analyze engineering problems such as metal volume plastic forming, the efficient compression storage and generation method of finite element stiffness matrix elements has become an urgent problem to be solved.

Contents of the invention

The purpose of the present invention is to provide a new finite element analysis stiffness matrix storage in metal volume plastic forming for the problems of large storage space of stiffness matrix data and low calculation efficiency of stiffness equation in the finite element numerical simulation process of metal volume plastic forming Compared with the generation method, it not only has the characteristics of significantly saving the storage space of the stiffness matrix and improving the calculation efficiency, but also has the more significant space saving and higher solution efficiency with the increase of the number of nodes and elements of the finite element mesh model of the metal volume plastic forming problem specialty.

The present invention is achieved through the following technical solutions:

A method for storing and generating a stiffness matrix for finite element analysis in metal volume plastic forming, comprising the following steps:

(1.1) First, the existing CAD design software or three-coordinate measuring instrument and other reverse equipment are used to obtain the three-dimensional geometric model of the metal volume plastic forming problem, and the existing grid generation method is used to generate the finite element mesh model of the formed workpiece. Form the topology of the finite element mesh model of the workpiece, determine the "generalized adjacent node pairs" of its mesh nodes, and establish the relationship between the "generalized adjacent node pairs" and the non-zero submatrix of the overall stiffness matrix of the finite element mesh model of the workpiece Correspondence; the "generalized adjacent node pair" refers to the node pair formed by any two nodes in each unit in the finite element mesh model of the volume forming problem, that is, the node itself and the node formed by itself right;

(1.2) Based on the "generalized adjacent node pair", store the non-zero sub-matrix, open up a non-zero sub-matrix column number storage array J, which is an integer index array, and record the "generalized adjacent node pair" in the finite element mesh model of the workpiece "Corresponding to the column number of the non-zero sub-matrix of the overall stiffness matrix, open up another integer index array K, record the position of the first non-zero sub-matrix in each row in the stiffness matrix of the workpiece finite element mesh model, and determine the workpiece The position and quantity of the non-zero sub-matrix in the stiffness matrix of the finite element grid model, according to the corresponding relationship between the non-zero sub-matrix and the non-zero elements, determine the position and quantity of the non-zero elements of the stiffness matrix of the finite element grid model of the workpiece;

(1.3) Calculate the non-zero sub-matrix of each unit stiffness matrix element in the workpiece finite element grid model in the overall stiffness matrix of the workpiece finite element grid model, and use the data in the index arrays J and K to determine the element stiffness matrix elements The location in the storage array A of the overall stiffness matrix of the grid model (see 4.1 to 4.9 for the specific determination steps); according to the generation rules of the overall stiffness matrix of the grid model, the stiffness matrix elements of each unit of the grid model are assembled into the grid The overall stiffness matrix of the lattice model is stored in the array A; after traversing all the cells in the grid model, the overall stiffness matrix storage array of the metal volume forming problem grid model is finally generated;

(1.4) Using the existing symmetric positive definite coefficient matrix relaxation pretreatment co-value gradient method, the overall stiffness equation (large linear equation system) of the mesh model is solved to obtain the velocity increase of all nodes in the mesh model during the metal volume forming process. Numerical solution of the quantity field.

In the described step (1.1), the method for determining the "generalized adjacent node pair" of the grid node is:

若不考虑节点的顺序，则构成一个“无序广义相邻节点对”B_IJ，即认为和

是同一个节点对，记为B_IJ或B_JI。为了应用和表示的方便，可采取下标由小到大的形式来表示“无序广义相邻节点对”，即统一记为B_IJ(I≤J)；(2.1) In the finite element mesh model of the volume forming problem, the node pair formed by any two nodes in each element (including the node itself and the node pair formed by itself) is called "generalized adjacent node pair"; The node pair with the order of adjacent nodes is called "ordered generalized adjacent node pair"; on the contrary, it is called "disordered generalized adjacent node pair". If the global node number of two nodes of a cell in the grid model is (I, J), then constitute two "ordered generalized adjacent node pairs" and

If the order of the nodes is not considered, then a "disordered generalized adjacent node pair" B _IJ is formed, that is, and

are the same node pair, denoted as B _IJ or B _JI . For the convenience of application and expression, the form of subscripts from small to large can be used to represent "disordered generalized adjacent node pairs", which is uniformly recorded as B _IJ (I≤J);

对应网格模型的总体刚度矩阵中第I行、第J列的非零子矩阵；(2.2) In the grid model, each "ordered generalized adjacent node pair" corresponds to a non-zero sub-matrix in the overall stiffness matrix of the grid model that has the same dimension as the analyzed problem, and the position of the non-zero sub-matrix Determined by the global node number of the "Ordered Generalized Adjacent Node Pair", e.g. "Ordered Generalized Adjacent Node Pair" in the grid model

The non-zero submatrix of row I and column J of the overall stiffness matrix corresponding to the grid model;

(2.3) In the grid model, each "disordered generalized adjacent node pair" corresponds to a non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix of the grid model that has the same dimension as the analyzed problem, non-zero The position of the sub-matrix is determined by the global node number of the "unordered generalized adjacent node pair"; for example, the "unordered generalized adjacent node pair" B _IJ (I≤J) in the grid model corresponds to the overall The non-zero submatrix of the I row and J column in the upper triangular matrix of the stiffness matrix;

(2.4) Since the overall stiffness matrix of the finite element analysis of the metal volume plastic forming process is a symmetrical matrix, in order to save space, only the upper triangular matrix can be stored, so the overall stiffness matrix storage of the grid model mentioned in the present invention refers to the upper The storage of the triangular matrix, because the "disordered generalized adjacent node pair" in the grid model corresponds to the non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix of the grid model, and the global node number of the node pair is the same as the non-zero sub-matrix The row number and column number of the zero submatrix correspond, so through this correspondence, the upper limit of the stiffness matrix of the grid model can be determined by the number of all "disordered generalized adjacent node pairs" in the grid model and the global node number. The number and position (row number and column number) of the non-zero sub-matrix in the triangular matrix, and can further determine the number and position of the non-zero elements that make up the non-zero sub-matrix.

In described step (1.2), also comprise the following steps:

(3.1) Use a floating-point array A to store the non-zero elements in the upper triangular matrix of the overall stiffness matrix of the mesh model in row order;

(3.2) Use the integer array J to record the column number of the non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix corresponding to the "disordered generalized adjacent node pair";

(3.3) Use the integer array K to record the number of the first non-zero sub-matrix in the array J of each row in the upper triangular matrix of the overall stiffness matrix, and use this to determine the row number of the non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix, Because K[i] is the number of the first non-zero sub-matrix in the array J of the i-th row in the upper triangular matrix of the overall stiffness matrix, and K[i+1] is the i+1th row in the upper triangular matrix of the overall stiffness matrix The number of the first non-zero sub-matrix in array J, it can be seen that all non-zero sub-moments corresponding to numbers from K[i] to K[i+1]-1 in array J are in the i-th row.

The finite element mesh model after the geometric model division of the analyzed problem is a three-dimensional hexahedral unit, tetrahedral unit or a mixed unit model of the two; or the finite element mesh model after the geometric model division of the analyzed problem is a two-dimensional The method of the present invention can be used to store the overall stiffness matrix of the mesh model.

In the step (1.3), for a unit in the grid model, if the number of nodes of the unit is n and the dimension of the physical quantity is m, record the overall number of the n nodes of the unit as I[i], i=1 ,2,Λ,n, the element stiffness matrix is U(m×n,m×n), then any element of the element stiffness matrix is

U _ijkl = U[(i-1)×m+k][(j-1)×m+l]

Wherein, i, j=1, 2, Λ, n are the local numbers of the unit nodes, k, l=1, 2, Λ, m are the physical quantity dimension numbers;

The process of determining the serial number of U _ijkl in the storage array A of the overall stiffness matrix of the grid model is as follows:

(4.1) Determine the overall node number of U _ijkl : I=I[i], J=I[j], where I and J are the non-zero sub-matrix A _IJ where U _ijkl is located on the overall stiffness matrix of the grid model Row and column numbers in the triangular matrix; the location of the nonzero submatrix is shown in Figure 1.

(4.2) Determine the number of non-zero sub-matrices before the non-zero sub-matrix A _II : N _I =K[I]-1;

(4.3) Determine the number of nonzero elements preceding the nonzero submatrix A _II :

NN _I ＝N _I ×m ² -(I-1)×m×(m-1)/2

Among them, m ² is the number of non-zero elements of each non-zero sub-matrix, and m(m-1)/2 is the reduction term of non-zero elements caused by the diagonal non-zero sub-matrix of each row; The zero sub-matrix A _II is a three-dimensional sub-matrix with m=3, and only 6 non-zero elements are stored, which is 3 less than other sub-matrixes.

(4.4) Determine the number of non-zero elements before the k-th row in the non-zero submatrix of the I-th row:

NN _IK ＝(K[I+1]-K[I])×m×(k-1)

(4.5) Determine the number of non-zero sub-matrices in row I before column J. Search the array J from the label K[I] to the array label II of J[II]=J in K[I+1]. In the I row, the number of non-zero sub-matrixes before the J column is N _J =II-K [I];

(4.6) Determine the number of U _ijkl in the row of the overall stiffness matrix:

NN _Jl =N _J ×m-TT+l

Wherein, TT is the reduction term of the non-zero elements due to the diagonal non-zero sub-matrix, TT=k(k-1)/2;

(4.7) Determine the number of U _ijkl in the storage array A of the overall stiffness matrix of the grid model:

NN _IJkl = NN _I +NN _Ik +NN _Jl

(4.8) According to the generation rules of the overall stiffness matrix of the grid model and the position number NN _IJkl of U _ijkl in the overall stiffness matrix storage array A of the grid model, assemble or accumulate U _ijkl into the overall stiffness matrix storage array A of the grid model middle;

(4.9) Repeat steps (4.1) to (4.8) to traverse each element in the stiffness matrix of all elements in the grid model to obtain the one-dimensional storage array A of the overall stiffness matrix of the grid model.

The numerical analysis method used to solve the grid model ultimately boils down to the storage and generation of the sparse and symmetrical overall stiffness matrix and the solution of its stiffness equation (ie, the linear equation system). The parameters for the solution of the overall stiffness equation of the analyzed grid model are the full amount or incremental.

The beneficial effects of the present invention are: according to the "generalized adjacent node pair" proposed by the present invention and its corresponding relationship with the non-zero sub-matrix of the overall stiffness matrix of the grid model, the non-zero sub-matrix in the overall stiffness matrix of the grid model is stored The matrix and its position have changed the previous method of storing non-zero elements and their positions in the overall stiffness matrix of the grid model; the storage space required by the compressed storage method of the present invention increases as a linear function with the increase of the number of grid model nodes, different Because the storage space required by the existing compressed storage method increases as a quadratic function with the increase of the number of grid model nodes, the method of the present invention not only significantly saves the storage space of the overall stiffness matrix of the grid model and improves the calculation efficiency, Moreover, as the number of grid model nodes and elements increases, its ability to save storage space becomes more significant; the relaxation preconditioning co-gradient method based on the overall stiffness matrix equation of the grid model storing non-zero sub-matrices does not need to perform symmetric stepwise over-relaxation The iterative splitting matrix is inverted and another array is opened to store the splitting matrix, and the non-zero elements of the splitting matrix can be directly operated, and the solution and calculation efficiency are high; therefore, it is especially suitable for metal volume forming that requires repeated incremental solutions and calculations. Numerical analysis of large complex engineering problems solved iteratively.

Description of drawings

Figure 1 is a schematic diagram of the distribution of non-zero sub-matrices in the overall stiffness matrix;

Fig. 2 is a schematic diagram of the positions of non-zero elements in the overall stiffness matrix;

Fig. 3 is a flow chart of the compression storage and generation algorithm for the overall stiffness matrix equation of the metal rigid (viscosity) plastic forming mesh model;

Fig. 4 is the grid model that is made of an eight-node hexahedron element;

Fig. 5 is the grid model that is made of two eight-node hexahedron elements;

Figure 6a is the finite element numerical simulation of the forging process;

Figure 6b is the finite element numerical simulation of the forging process;

Figure 6c is the finite element numerical simulation of the forging process;

Figure 6d is the finite element numerical simulation of the forging process.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

Figure 3 is a flow chart of the compression storage and generation algorithm for the overall stiffness matrix equation of the metal rigid (viscosity) plastic forming mesh model. As shown in Figure 3, the overall stiffness matrix equation compression storage and generation algorithm flow of the metal rigid (viscosity) plastic forming mesh model is as follows: According to the generated geometric model and finite element mesh model of the metal plastic forming problem, determine the entire mesh Calculate the number of "generalized adjacent node pairs" of the model; open up the overall stiffness matrix storage array A of the grid model, open up the storage array J for the column number of the non-zero sub-matrix, and open up the storage array K for the position of the first non-zero sub-matrix in each row; Loop through all the nodes in the grid model, for any node i, search all the cells that contain node i in the grid model; search all the nodes that contain all the cells that contain node i, and determine the "generalized neighbor" of node i node pair"; according to the number of "generalized adjacent node pairs" containing node i, an index array J is generated; according to the number of "generalized adjacent node pairs" containing node i, an index array K is generated; to determine the number of all nodes in the grid model Whether the loop is over, if not, repeat the above steps for the next node until all the nodes in the grid model are processed. Cycle through all the units in the grid model to generate the stiffness matrix U of each unit, and determine the position of each element of the unit stiffness matrix in the triangular matrix of the overall stiffness matrix of the grid model. The number of non-zero sub-matrices and non-zero elements before the row where the zero sub-matrix is located; calculate the number of all non-zero elements before the row where the elements of the triangular matrix are located in the overall stiffness matrix of the grid model corresponding to the elements of the element stiffness matrix; calculate the element stiffness The sequence number of the row where the triangular matrix element of the overall stiffness matrix of the grid model corresponding to the matrix element is located; the position of the upper triangular matrix element of the overall stiffness matrix of the grid model corresponding to the element stiffness matrix element in the storage array A of the overall stiffness matrix of the grid model According to the generation rules of the overall stiffness matrix of the grid model, the element stiffness matrix elements are loaded into the upper triangular matrix storage array of the overall stiffness matrix of the grid model. It is judged whether all the elements in the unit matrix have been calculated, and if not, the above processing is performed on the next element in the unit until all the elements in the unit are processed. Then judge whether there are still units that have not generated the stiffness matrix, and if there are still, perform the cyclic calculation of the next unit until the stiffness matrices of all units are generated according to the above steps, so as to obtain the upper triangular matrix of the overall stiffness matrix of the entire grid model and its stiffness equation. Finally, the overall stiffness equation of the entire grid model is solved by using the relaxation preconditioning co-gradient method to obtain the incremental solution of the velocity field, and then use this incremental solution to refresh the velocity field of the grid model at this moment.

The number of "generalized adjacent node pairs" in the grid model is determined by the topology of the grid model. Any two nodes belonging to the same unit can form a "generalized adjacent node pair", so the "generalized adjacent node pair The determination of the "quantity" is to calculate the number of non-repeated node pairs belonging to the same unit in the entire grid model. For the grid models of different unit types, the principle of determining "generalized adjacent node pairs" is the same, but the specific methods are as follows: different.

Since the "generalized adjacent node pairs" correspond to the non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix of the mesh model, the non-zero sub-matrix in the upper triangular matrix of the overall stiffness matrix can be determined by the number of "generalized adjacent node pairs". The number of matrices, which in turn determines the number of nonzero elements. However, since the overall stiffness matrix of the grid model is stored in an upper triangular matrix, each non-zero sub-matrix on the diagonal only needs to store an upper triangular sub-matrix. If m is the dimension of the problem to be solved, the non-zero elements of the diagonal non-zero sub-matrix are m(m+1)/2, and the non-diagonal non-zero sub-matrix data is m ² , so it can be determined The number of non-zero elements of the triangular matrix in the overall stiffness matrix of the mesh model.

The above calculation process is illustrated by taking a mesh model composed of a finite number of eight-node hexahedral elements as an example.

For the mesh model composed of eight-node hexahedral elements, its "generalized adjacent node pairs" are divided into three types, one is edge adjacent node pairs; the other is surface adjacent node pairs; the third is volume adjacent node pairs. The edge adjacent node pairs are determined by the number of adjacent node connections in the unit, and each connection corresponds to a "generalized adjacent node pair"; the number of surface adjacent node pairs is determined by the number of cell faces, and each face has two diagonals, and each face diagonal corresponds to a "generalized adjacent node pair"; for body adjacent node pairs, each unit has four body diagonals, and each body diagonal corresponds to a "generalized adjacent node pair". node pair". Let N ₁ be the number of nodes (that is, the number of node pairs generated by the node and itself), N ₂ is the number of adjacent node connections in the grid model (that is, the number of edge-adjacent node pairs), and N ₃ is the grid The number of faces in the model (each face has two face-adjacent node pairs, 2N ₃ represents the number of face-adjacent node pairs), N ₄ is the number of cells in the mesh model (each cell has two or four volume phases Adjacent node pairs, 4N ₄ represents the number of volume adjacent node pairs), then the number of "generalized adjacent node pairs" in the grid model is:

S _SM ＝N ₁ +N ₂ +2×N ₃ +4×N ₄

Among them, N ₁ is the number of node pairs formed by the node itself. The row number of the non-zero sub-matrix corresponding to this type of node pair is the same as the column number, and is located on the diagonal, so N ₁ is the number of non-zero sub-matrix on the diagonal.

For three-dimensional problems (m=3), the data to be stored in each diagonal non-zero sub-matrix is m(m+1)/2=6, and the data to be stored in off-diagonal non-zero sub-matrix is m ² =9. Therefore, the number of non-zero elements in the triangular matrix of the overall stiffness matrix of the mesh model is:

S _NZ ＝6×N ₁ +9×(N ₂ +2×N ₃ +4×N ₄ )

In the following, two specific grid models are taken as examples to further illustrate the above calculation process.

For the grid model composed of a hexahedral unit shown in Figure 4, N ₁ =8, N ₂ =12, N ₃ =6, N ₄ =1, when the upper triangular matrix is used for storage, the number of non-zero sub-matrices for:

S _SM ＝N ₁ +(N ₂ +2×N ₃ +4×N ₄ )=8+(12+2×6+1×4)=36

The grid model shown in Figure 4 needs to save the number of non-zero elements in the upper triangular matrix of the overall stiffness matrix:

S _NZ ＝6×N ₁ +9×(N ₂ +2×N ₃ +4×N ₄ )=6×8+9×28=300

Similarly, for the grid model composed of two hexahedral units shown in Fig. 5, N ₁ =12, N ₂ =20, N ₃ =11, and N ₄ =2. When the upper triangular matrix is used for storage, the number of non-zero sub-matrices is:

S _SM ＝N ₁ +(N ₂ +2×N ₃ +4×N ₄ )=12+(20+2×11+2×4)=62

The number of non-zero elements in the triangular matrix of the overall stiffness matrix that needs to be saved in the grid model shown in Figure 5 is:

S _NZ ＝6×N ₁ +9×(N ₂ +2×N ₃ +4×N ₄ )=6×12+9×50=522

After determining the number of non-zero sub-matrixes and non-zero elements, the storage space of the upper triangular matrix of the overall stiffness matrix of the grid model can be allocated accurately, and then the storage of the upper triangular matrix of the overall stiffness matrix of the grid model can be realized.

Taking the overall stiffness matrix of the grid model composed of two hexahedral elements shown in FIG. 5 as an example, the advantages of compression storage of the overall stiffness matrix of the present invention will be described.

The data length of the J array in the overall stiffness matrix storage method involved in the present invention is the quantity of the non-zero sub-matrix of the overall stiffness matrix, that is, S _SM =62, and the data length of the K array is the number of cells 12, then the grid model The storage data volume of the overall stiffness matrix is 62+12=74. If the currently commonly used CSR (Compressed Sparse Row format) format storage method is used, the data length of the J array is the number of all non-zero elements, that is, S _NZ = 522, and the data length of the K array is the dimension of the overall stiffness equation. That is, 12×3=36, therefore, the total amount of data is 558. If only two auxiliary arrays J and K are compared, compared with the CSR method, the overall stiffness matrix storage method of the present invention saves storage space as: (558-74)/558=87%. If the floating-point array A adopts single-precision floating-point numbers to save, then the overall stiffness matrix storage method involved in the present invention is compared with the CSR method, and the saved space of the overall stiffness matrix storage is: 1-(74+552×2)/ (558+552×2)=29%. It can be seen that for the example shown in Figure 5, the present invention saves 87% of the storage space of the auxiliary arrays J and K, and saves 29% of the overall stiffness matrix storage space, and the space saving effect is obvious.

The beneficial effects of the present invention will be further described through a specific forging example below.

Accompanying drawings 6a-6d are the finite element numerical simulation results obtained by using different element division numbers at a certain moment in the forging process of a forging, in which all the elements of the mesh model are eight-node hexahedral elements, and the forging mesh in Fig. 6a The lattice model contains 1565 units and 2045 nodes. The forging mesh model in Fig. 6b contains 3284 units and 4119 nodes. The forging mesh model in Fig. 6c contains 5339 units and 6557 nodes. The forging mesh model in Fig. 6d The mesh model contains 8530 elements and 10217 nodes. In the following, the compressed storage iterative solution algorithm of the present invention is compared with the half-band storage direct solution algorithm from two aspects of memory occupation and calculation efficiency. The computer used is: CPU: Intel(R) Core(TM) 2, Duo, CPU E8400, 3.00GHz, memory: 3.25, operating system: Windows XP.

In terms of memory usage, the elements of the overall stiffness equation of the forging mesh model are stored in single-precision floating-point numbers, that is, each data occupies 4 bytes; the auxiliary array is stored in integer data, that is, each data occupies 2 bytes. Under different numbers of nodes and elements, the storage space occupied by the overall stiffness matrix of the forging mesh model stored in different ways is compared, and the results are shown in Table 1.

Table 1 Comparison of the overall stiffness matrix storage space of the forging mesh model under different storage methods

As can be seen from Table 1, the storage space occupied by the compressed storage method of the present invention is far less than that of the half-band storage method. When the number of units of the forging grid model is 1565, the storage space is saved by 93.2%. When the number of units is 8530, the storage space is saved by 97.4%. It can be seen that with the increase of the number of grid model nodes or units, the proportion of storage space saving becomes larger and larger. This is mainly because the compressed storage method of the present invention only saves the non-zero elements of the overall stiffness matrix, and the number of non-zero elements in each row of the overall stiffness matrix is only related to the number of generalized adjacent nodes of the corresponding node in the row, but not related to the overall network The number of nodes in the lattice model is irrelevant. In other words, for the compressed storage method of the present invention, the proportional relationship between the storage space and the number of nodes is only related to the increase in the number of rows of the overall stiffness matrix, and shows a linear trend. For the half-band storage method, the storage space of each row is determined by the half-bandwidth of the overall stiffness matrix. As the number of nodes increases, the half-bandwidth generally increases, and the number of rows of the overall stiffness matrix and the data that needs to be stored in each row are different. As a result, it presents a parabolic trend. Therefore, as the number of nodes increases, the storage space required by the half-band storage method increases with a quadratic function, while the storage space required by the compressed storage method of the present invention increases with a linear function, and its space-saving advantage is more significant.

In terms of calculation efficiency, here we only take the single calculation efficiency of the overall stiffness equation at a certain moment in the forging process as an example, and compare the two methods of half-band storage direct solution and compressed storage iterative solution of the present invention. The results are shown in the attached Table 2 shows.

Table 2 compares the solution efficiency of mesh models with different numbers of elements and nodes in different storage methods

It can be seen from the table that when the number of elements of the forging mesh model is 1565, the calculation efficiency is increased by 69.4%, and when the number of elements of the forging mesh model is 8530, the calculation efficiency is increased by 88.3%. It can be seen that the efficiency of the compressed storage iterative solution of the present invention is much higher than that of the half-band storage direct solution, and the greater the number of nodes, the more obvious the improvement of the solution efficiency. Therefore, the method proposed by the present invention is particularly suitable for the finite element numerical analysis of metal volume forming problems such as forging, extrusion, and rolling, which have a large number of mesh model units and nodes and are repeatedly solved iteratively.

Claims (6)

1. a stiffness matrix for finite-element analysis in metal bulk plastic forming is stored and the generation method, it is characterized in that comprising the following steps:

(1.1) at first adopt existing CAD software or three-coordinates measuring machine to obtain the 3-D geometric model of metal volume Plastic Forming Problems; Adopt existing grid generation method to generate the finite element grid model of Forming Workpiece; Next is according to the finite element grid model topology structure of Plastic Forming workpiece, determine " the broad sense adjacent node to " of its grid node, set up the corresponding relation between " broad sense adjacent node to " and workpiece finite element grid model global stiffness matrix non-zero submatrices; Described " broad sense adjacent node to " refers in the finite element grid model of bulk forming problem, and the node that any two nodes in each unit consist of pair namely comprises the node pair of node self and self formation;

(1.2) based on " broad sense adjacent node to ", store non-zero submatrices, opening up non-zero submatrices row number storage array J is the integer array of indexes, the row of the global stiffness matrix non-zero submatrices that record is corresponding with " broad sense adjacent node to " in workpiece finite element grid model number, open up another integer array of indexes K, record the position of the first non-zero submatrices of every row in array J in workpiece finite element grid rigidity of model matrix upper triangular matrix, determine position and the quantity of non-zero submatrices in workpiece finite element grid rigidity of model matrix, corresponding relation according to non-zero submatrices and nonzero element, determine position and the quantity of workpiece finite element grid rigidity of model matrix nonzero element,

(1.3) calculate the non-zero submatrices at each element stiffness matrix element place in workpiece finite element grid model global stiffness matrix in workpiece finite element grid model, utilize the data in array of indexes J and K, the position of determining unit Element of Stiffness Matrix in grid model global stiffness matrix stores array A; According to the create-rule of grid model global stiffness matrix, each element stiffness matrix element of grid model is assembled in the global stiffness matrix stores array A of grid model with sitting in the right seat; Behind all unit in the traversal grid model, finally generate the global stiffness matrix stores array of metal volume Problems in forming grid model;

(1.4) adopt the lax pre-service of existing symmetric positive definite matrix of coefficients volume gradient method altogether, global stiffness equation to grid model is that large linear systems is found the solution, and obtains the numerical solution of the speed increment field of all nodes in metal volume forming process grid model.

2. stiffness matrix for finite-element analysis in metal bulk plastic forming according to claim 1 storage and generation method, is characterized in that, in described step (1.1), determines that the method for " the broad sense adjacent node to " of grid node is:

(2.1) in the finite element grid model of bulk forming problem, the node that any two nodes in each unit consist of pair namely comprises the node pair of node self and self formation, is called " broad sense adjacent node to "; With the sequential node of adjacent node to being called " in order broad sense adjacent node to "; Otherwise, be called " unordered broad sense adjacent node to ", if in grid model, the global node of certain two node in unit is numbered (I, J), consist of two " broad sense adjacent node to " in order

With

If do not consider the order of node, consist of one " unordered broad sense adjacent node to "

Namely think

With

Be same node pair, be designated as

Or

For the convenience of using and representing, can take the ascending form of subscript to represent " unordered broad sense adjacent node to ", i.e. unified being designated as

(I≤J);

(2.2) in grid model, in the global stiffness matrix of each " in order broad sense adjacent node to " corresponding grid model one non-zero submatrices identical with institute problem analysis dimension, the position of non-zero submatrices is by the global node sequence number decision of " broad sense adjacent node to " in order;

(2.3) in grid model, in the upper triangular matrix of the global stiffness matrix of each " unordered broad sense adjacent node to " corresponding grid model one non-zero submatrices identical with institute problem analysis dimension, the position of non-zero submatrices is determined by the global node sequence number of " unordered broad sense adjacent node to ";

(2.4) the global stiffness matrix of metal volume plastic forming process finite element analysis is symmetric matrix, only stores the triangular matrix on it, non-zero submatrices in grid model in the upper triangular matrix of " unordered broad sense adjacent node to " and the global stiffness matrix of grid model is corresponding one by one, and the global node numbering that node is right and the line number of non-zero submatrices, the row correspondence, therefore by this corresponding relation, quantity and position by non-zero submatrices in the upper triangular matrix of the quantity of all in grid model " unordered broad sense adjacent nodes to " and global node number definite grid model stiffness matrix, be line number and row number, thereby further determine quantity and the position of the nonzero element of composition non-zero submatrices.

3. stiffness matrix for finite-element analysis in metal bulk plastic forming storage according to claim 1 and generation method, is characterized in that, and be in described step (1.2), further comprising the steps of:

(3.1) adopt a floating type array A by each nonzero element in row sequential storage grid model global stiffness matrix upper triangular matrix;

(3.2) adopt the row number of the sub-square of non-zero in shaping array J record and " unordered broad sense adjacent node to " corresponding global stiffness matrix upper triangular matrix;

(3.3) adopt shaping array K to record the numbering of first non-zero submatrices of every row in array J in global stiffness matrix upper triangular matrix, and determine that with this non-zero submatrices is in the line number of global stiffness matrix upper triangular matrix, because K[i] be the numbering of capable first non-zero submatrices of i in array J in global stiffness matrix upper triangular matrix, and K[i+1] be the numbering of capable first non-zero submatrices of i+1 in array J in global stiffness matrix upper triangular matrix, thus in array J numbering from K[i] to K[i+1] the sub-square of all non-zeros of-1 correspondence is all capable at i.

4. stiffness matrix for finite-element analysis in metal bulk plastic forming storage according to claim 1 and generation method, it is characterized in that, the finite element grid model after the geometric model of institute's problem analysis is divided is three-dimensional hexahedral element, tetrahedron element or its mixed cell model both; Finite element grid model after perhaps the geometric model of institute's problem analysis is divided is quadrilateral units, triangular element or its mixed cell model both of two dimension.

5. stiffness matrix for finite-element analysis in metal bulk plastic forming storage according to claim 1 and generation method, is characterized in that, in described step (1.3), for a unit in grid model, if the dimension that the nodes of unit is n, physical quantity is m, remember that the integral body of n the node in this unit is numbered I[i], i=1,2, ∧, n, this element stiffness matrix are U(m * n, m * n), the arbitrary element of this element stiffness matrix is

U _ijkl=U[(i-1)×m+k][(j-1)×m+1]

Wherein, i, j=1,2, ∧, n are the part numbering of this cell node, and k, l=1,2, ∧, m are physical quantity dimension numbering; Determine U _IjklThe process of sequence number is as follows in grid model global stiffness matrix stores array A:

(4.1) determine U _IjklIntegral node numbering: I=I[i], J=I[j], wherein I, J are U _IjklThe non-zero submatrices A at place _IJIn the line number in grid model global stiffness matrix upper triangular matrix and row number;

(4.2) determine non-zero submatrices A _IIThe quantity of non-zero submatrices: N before _I=K[I]-1;

(4.3) determine non-zero submatrices A _IIThe quantity of nonzero element before:

NN _I=N _I×m ²-(I-1)×m×(m-1)/2

Wherein, m ²Be the quantity of each non-zero submatrices nonzero element, the minimizing item of the nonzero element that m (m-1)/2 causes for the diagonal angle non-zero submatrices of every row; Non-zero submatrices A _IIBe the three-dimensional submatrix of m=3, only stored 6 nonzero elements, compare with other submatrix and deposit less 3;

(4.4) determine the capable quantity of nonzero element before of k in the capable non-zero submatrices of I:

NN _IK＝K[I+1]-K[I]×m×(k-1)

(4.5) determine that I is capable, the quantity of the non-zero submatrices before the J row.Search array J is from label K[I] to K[I+1] J[II]=the array label II of J, during I was capable, it was N that J is listed as non-zero submatrices quantity before _J=II-K[I];

(4.6) determine U _IjklNumbering in being expert in the global stiffness matrix:

NN _J1＝N _J×m-TT+1

Wherein, TT is the minimizing item of the nonzero element that causes due to the diagonal angle non-zero submatrices, TT=k (k-1)/2;

(4.7) determine U _IjklNumbering in grid model global stiffness matrix storage arrays A:

NN _IJk1=NN _I+NN _IK+NN _Jl

(4.8) according to create-rule and the U of grid model global stiffness matrix _IjklPosition Number NN in grid model global stiffness matrix storage arrays A _IJK1, with U _IjklAssemble or be added in the global stiffness matrix storage arrays A of grid model;

(4.9) repeating step (4.1) ~ (4.8), each element in all element stiffness matrixs in the traversal grid model, the global stiffness matrix one dimension that can obtain grid model is stored array A.

6. stiffness matrix for finite-element analysis in metal bulk plastic forming storage according to claim 1 and generation method, it is characterized in that: it is Solving Linear that the numerical analysis method that the analysis grid model adopts finally is summed up as storage, generation and stiffness equations thereof sparse, the symmetric population stiffness matrix, and the global stiffness equation solution parameter of the grid model of analyzing is full dose or increment.

Images