CN102184298B - Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming - Google Patents
Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming Download PDFInfo
- Publication number
- CN102184298B CN102184298B CN 201110130073 CN201110130073A CN102184298B CN 102184298 B CN102184298 B CN 102184298B CN 201110130073 CN201110130073 CN 201110130073 CN 201110130073 A CN201110130073 A CN 201110130073A CN 102184298 B CN102184298 B CN 102184298B
- Authority
- CN
- China
- Prior art keywords
- grid model
- stiffness matrix
- node
- matrix
- zero submatrices
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Complex Calculations (AREA)
Abstract
The invention relates to a method for storing and generating a stiffness matrix for finite-element analysis in metal bulk plastic forming. The method capable of remarkably saving a storage space of the stiffness matrix and improving the solving efficiency of a stiffness equation comprises the following steps of: (1.1) determining a generalized adjacent node pair of a grid node and a correspondence relationship between the generalized adjacent node pair and non-zero sub-matrixes in a total stiffness matrix of a grid model according to a finite-element grid model of a bulk plastic forming workpiece; (1.2) determining the positions and quantity of the non-zero sub-matrixes in the total stiffness matrix of the grid model based on the generalized adjacent node pair, and storing the non-zero sub-matrixes; (1.3) determining the non-zero sub-matrix of the total stiffness matrix in which each unit stiffness matrix element of the grid model is located and the position of each unit stiffness matrix element in a total stiffness matrix storage array, and assembling each unit stiffness matrix element into the total stiffness matrix storage array of the grid model; and (1.4) solving a total stiffness equation of the grid model by using a conjugate gradient method for relaxation pre-treatment of symmetrical positive definite coefficient matrixes.
Description
Technical field
The present invention relates to stiffness matrix for finite-element analysis in metal bulk plastic forming storage and generation method, be suitable for forging, the finite element numerical simulation of extruding, rolling, swaging equal-volume forming process.
Background technology
In commercial production, forging, extruding, the metal volume Plastic Forming manufacturing technology such as rolling support one of major technique of the national economic development and national defense construction, have extensive use in industries such as automobile, space flight and aviation, equipment manufacturing, weapons, the energy, shipbuildings.Plastic yield occurs with the metal material (blank) of (cold, warm, thermoforming) under room temperature or high temperature in volume plastic forming process exactly under forging equipment and the mould action that is installed on it thereof, and be configured as the process of the part of required form and performance, volume Plastic Forming material generally is assumed to be and just (glues) plastic material, therefore, large and violent plastic yield occurs in blank in mold cavity, have material nonlinearity, geometrical non-linearity and complicated boundary condition, finding the solution of forming process is very complicated.Can forming technology and Design of Dies be for the quality of drip molding and be shaped smoothly and have material impact, are also the keys of the energy-saving and cost-reducing and Accurate Shaping of impact.The key theory foundation of design forming technology and mould and grasp plastic yield behavior and the flowing law of metal material in mold cavity.Finite Element Method is as the important method of engineering problem numerical analysis, brought into play vital role at engineering field, metal volume plastic forming process numerical simulation method is the main method of research Plastic Deformation of Metal Materials behavior and flowing law, can obtain metal material flowing law and various field variable distribution results, the feasibility of checking technological design, parameter selection and die design schemes, and then realize the optimization optimization of forming technology and mould improving speed, quality and the reduction cost of development of product development.
Finite element analysis for the metal volume Plastic Forming Problems, at first by the plastic yield governing equation, set up the mathematical model of metal volume plastic yield mechanics problem, geometric model (mould and sex change material) is carried out finite element method discrete, set up the finite element stiffness equations, realize the numerical simulation of forming process, the numerical analysis methods such as finite element have particularly been brought into play vital role in fields such as metal forming manufacturing and structure analyses at engineering field.
The step that finite element numerical is found the solution the metal volume Plastic Forming Problems is: the CAD geometric model (geometric configuration of deformable material) of setting up problem, the CAD geometric model is carried out finite elements discrete, obtain the finite element grid model of geometric model, this grid model is comprised of the grid cell of limited quantity; The stiffness equations of each unit in the model grid model is assembled the stiffness equations of all unit, forms the non-linear global stiffness system of equations about all node speeds of grid model; By to the linearization process about the non-linear rigidity equation of node speed, obtain the linear rigidity system of equations about the node speed increment; Adopt the Newton-Raphson process of iteration to carry out iterative to velocity field.due to the node speed increment is found the solution, therefore, for a metal volume plastic forming process, curring time need to be decomposed into an even hundreds of time step up to a hundred (length), all need computing unit stiffness matrix and nodal force vector in each time step, be assembled into the global stiffness equation of node speed increment, apply velocity boundary conditions, the global stiffness equation (large linear systems) of solution node speed increment, obtain the speed correction, with the speed correction, velocity field is revised again, then repeat above-mentioned steps, velocity field is carried out iterative, until velocity field convergence, generally, at least the velocity field that needs 7 to 10 times iteration can obtain to restrain, according to the speed of convergence field that obtains, refresh distortion workpiece shape and die location, just can complete the numerical evaluation in a time step.And then carry out the interior numerical evaluation of next step-length, until reach desired deformation extent.
This shows, for the numerical simulation of metal volume plastic forming process, need to find the solution frequently global stiffness equation (large linear systems).For example for the 3-dimensional metal volume Plastic Forming Problems of medium complexity, generally need the simulation step number more than 200 times, each roughly needs 9 velocity fields of iteration in the simulation step, and each velocity field iteration just need to be found the solution large linear systems one time, if the finite element interstitial content of a metal volume Plastic Forming Problems is 50,000 (three-dimensional problems of medium complexity), just need to find the solution the large Linear side group that consisted of by 150,000 equations 1800 times.If analyze more challenge, its amount of calculation is huger.
In recent years, needs along with heavy constructions such as Grand Equipments, boats and ships, automobile and wind-powered electricity generations, increasing to demand large-scale, heavy duty metal volume Plastic Forming part, need the Problems in forming of analysis day by day complicated, physical dimension is large, node and the element number of finite element model are many, data volume is huge, the FEM (finite element) calculation inefficiency, required calculator memory increases severely, even when analyzing large-scale and ultra-large type problem, the problem of memory space inadequate often occurs, this has limited the application of finite element numerical method in METHOD IN METAL FORMING PROCESSES is analyzed greatly.
Metal volume Plastic Forming Problems not only, and much the finite element analysis of other engineering problem finally all can be summed up as and finds the solution large linear systems, and the solution efficiency of finite element depends primarily on the solution efficiency of storage, generation and the stiffness equations thereof of system of linear equations stiffness matrix.
metal volume Plastic Forming Problems finite element matrix is large-scale, banded sparse matrix, nonzeros generally concentrates on the diagonal line zone, for this situation with in order to improve counting yield, people have proposed the storage means of many conserve storage, for example, the half-band width storage means, the variable bandwidth storage means, sparse simultaneous equation pretreatment iteration method for solving etc., although these methods have reduced the storage space of stiffness matrix, but stored the part nonzero element also unavoidablely, wasted the part storage area, these methods all depend on node serial number simultaneously, need to be optimized the node serial number order, thereby also wasted computing time.Existing one-dimension variable bandwidth compression and storage method adopts one-dimension array to preserve the nonzero element of stiffness matrix by the order of row, adopt two auxiliary arrays record row number and the initial element of every row in the position of stiffness matrix array, saved storage space, avoided the node numbering optimization problem, but the auxiliary array that adopts and stiffness matrix storage array have same length, still have the problem that consumes storage space and waste computing time.
Therefore, during the engineering problems such as employing finite element analysis metal volume Plastic Forming, the Efficient Compression storage of finite element matrix element becomes problem in the urgent need to address with the generation method.
Summary of the invention
The objective of the invention is and stiffness equations counting yield not high problem large for the stiffness matrix data space that exists in metal volume Plastic Forming Problems finite element numerical simulation process, a kind of new stiffness matrix for finite-element analysis in metal bulk plastic forming storage and generation method are provided, it not only has remarkable saving stiffness matrix storage space and improves the characteristics of counting yield, and has an increase with metal volume Plastic Forming Problems finite element grid model node and element number, the more remarkable and higher characteristics of solution efficiency of its conserve space.
The present invention realizes by following technical scheme:
A kind of stiffness matrix for finite-element analysis in metal bulk plastic forming storage and generation method comprise the following steps:
(1.1) at first adopt the reverse equipment such as existing CAD design 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, 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 (concrete determining step sees 4.1～4.9); 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 (large linear systems) to grid model is found the solution, and obtains the numerical solution of the speed increment field of all nodes in metal volume forming process grid model.
In described step (1.1), determine 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 is to the node of node self and self formation (comprise to), 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 " unordered broad sense adjacent node a to " B
IJ, namely think
With
Be same node pair, be designated as B
IJOr B
JIFor the convenience of using and representing, can take the ascending form of subscript to represent " unordered broad sense adjacent node to ", i.e. the unified B that is designated as
IJ(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 determined by the global node sequence number of " in order broad sense adjacent node to ", for example, " in order broad sense adjacent node to " in grid model
The non-zero submatrices that in the global stiffness matrix of corresponding grid model, I is capable, J is listed as;
(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 "; " unordered broad sense adjacent node to " B in grid model for example
IJ(I≤J), the non-zero submatrices that in the upper triangular matrix of the global stiffness matrix of corresponding grid model, I is capable, J is listed as;
(2.4) the global stiffness matrix due to the finite element analysis of metal volume plastic forming process is symmetric matrix, for conserve space can only be stored its upper triangular matrix, the global stiffness matrix stores of the grid model of therefore mentioning in the present invention refers to the storage of its upper triangular matrix, because the non-zero submatrices in the upper triangular matrix of " unordered broad sense adjacent node to " in grid model 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, number can be determined the quantity of non-zero submatrices in the upper triangular matrix of grid model stiffness matrix and position (line number with row number) by the quantity of all in grid model " unordered broad sense adjacent nodes to " and global node, and can further determine to form quantity and the position of the nonzero element of non-zero submatrices.
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 non-zero submatrices 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, numbering that hence one can see that in array J is from K[i] to K[i+1] the sub-square of all non-zeros of-1 correspondence is all capable at i.
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, can adopt method of the present invention to carry out the storage of grid model global stiffness matrix.
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+l]
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; The position of non-zero submatrices as shown in Figure 1.
(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
2-(I-1)×m×(m-1)/2
Wherein, m
2Be 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; As the non-zero submatrices A in Fig. 2
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
Jl＝N
J×m-TT+l
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
IJkl＝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
IJkl, 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.
Find the solution numerical analysis method that grid model adopts and finally be summed up as storage, generation and stiffness equations thereof (being system of linear equations) sparse, the symmetric population stiffness matrix and find the solution, the parameter of the global stiffness equation solution of the grid model of analyzing is full dose or increment.
The invention has the beneficial effects as follows: propose according to the present invention " the broad sense adjacent node to " and and grid model global stiffness matrix non-zero submatrices between corresponding relation, non-zero submatrices and position thereof in save mesh model global stiffness matrix have changed the method for nonzero element in save mesh model global stiffness matrix in the past and position thereof; The required storage space of compression and storage method of the present invention is linear function with the increase of grid model number of nodes to be increased, being different from the required storage space of existing compression and storage method is quadratic function with the increase of grid model number of nodes and increases, therefore, method of the present invention has not only significantly been saved the storage space of grid model global stiffness matrix and has been improved operation efficiency, and with the increase of grid model node and element number, the ability of its conserve storage is more remarkable; Lax pre-service based on the grid model global stiffness matrix equation of storing non-zero submatrices is total to the volume gradient method, need not inverts and open up array in addition the division matrix of progressively overrelaxation iteration of symmetry deposits the division matrix, can directly carry out computing to the division nonzeros, find the solution with operation efficiency high; Therefore, being particularly suitable for metal volume is shaped and thisly need to repeatedly carries out the numerical analysis of the large complicated engineering problem of Incremental SAT and iterative.
Description of drawings
Fig. 1 is the distribution schematic diagram of non-zero submatrices in the global stiffness matrix;
Fig. 2 is the position view of nonzero element in the global stiffness matrix;
Fig. 3 is metal just (gluing) Plastic Forming grid model global stiffness matrix equation compression storage and generating algorithm process flow diagram;
The grid model of Fig. 4 for being consisted of by eight node hexahedral elements;
The grid model of Fig. 5 for being consisted of by two eight node hexahedral elements;
Fig. 6 a is the forging process finite element numerical simulation;
Fig. 6 b is the forging process finite element numerical simulation;
Fig. 6 c is the forging process finite element numerical simulation;
Fig. 6 d is the forging process finite element numerical simulation.
Embodiment
Below in conjunction with accompanying drawing and embodiment, the present invention will be further described.
Fig. 3 is metal just (gluing) Plastic Forming grid model global stiffness matrix equation compression storage and generating algorithm process flow diagram.According to shown in Figure 3, just (glue) Plastic Forming grid model global stiffness matrix equation compression of metal store and the generating algorithm flow process as follows: according to geometric model and the finite element grid model of the metal forming problem that has generated, determine " the broad sense adjacent node to " of whole grid model and calculate its quantity; Open up grid model global stiffness matrix stores array A, open up non-zero submatrices row number storage array J, open up the first non-zero submatrices of every row position storage array K; All nodes in grid model are circulated, for arbitrary node i, comprise all unit of node i in the search grid model; Search comprises all nodes that all unit of node i comprise, and determines " the broad sense adjacent node to " of node i; According to the numbering that comprises node i " broad sense adjacent node to ", generating indexes array J; According to the quantity that comprises node i " broad sense adjacent node to ", generating indexes array K; Whether the circulation that judges all nodes in grid model finishes, if also do not finish, next node is repeated above-mentioned steps, until complete to all node processing in grid model.All unit in grid model are circulated, generate the stiffness matrix U of each unit, place non-zero submatrices be expert at non-zero submatrices before and the quantity of nonzero element are calculated in the position of each element of judging unit stiffness matrix non-zero submatrices at place in grid model global stiffness matrix upper triangular matrix; The quantity of all nonzero elements before the element of the grid model global stiffness matrix upper triangular matrix that the computing unit Element of Stiffness Matrix is corresponding is expert at; The serial number that the grid model global stiffness matrix upper triangular matrix element that the computing unit Element of Stiffness Matrix is corresponding is expert at; The Position Number of grid model global stiffness matrix upper triangular matrix element in grid model global stiffness matrix stores array A that the computing unit Element of Stiffness Matrix is corresponding, according to the create-rule of grid model global stiffness matrix, with this element stiffness matrix element grid model global stiffness matrix upper triangular matrix storage array of packing into.Judge in this cell matrix whether all elements calculates complete, if not yet complete, the next element in this unit is carried out above-mentioned processing, until all elements in this unit is disposed.Then judge whether to still have the unit that does not generate stiffness matrix, if also have, carry out the cycle calculations of next unit, until the stiffness matrix of all unit is complete according to the above-mentioned steps generation, thereby obtain upper triangular matrix and the stiffness equations thereof of whole grid model global stiffness matrix.At last, adopt lax pre-service volume gradient method altogether, find the solution the global stiffness equation of whole grid model, obtain the Incremental of velocity field, then utilize this Incremental to refresh the velocity field of grid model in this moment.
In grid model, the quantity of " broad sense adjacent node to " is determined by the topological structure of grid model, any two nodes that belong to same unit can form one " broad sense adjacent node to ", therefore " broad sense adjacent node to " quantity determines, calculate exactly in whole grid model, the unduplicated node of same unit that belongs to is to quantity, for the grid model of different units type, determine that the principle of " broad sense adjacent node to " is identical, but concrete grammar is different.
Because non-zero submatrices in the global stiffness matrix upper triangular matrix of " broad sense adjacent node to " and grid model is corresponding one by one, can be determined by the quantity of " broad sense adjacent node to " quantity of non-zero submatrices in global stiffness matrix upper triangular matrix and then the quantity of definite nonzero element.Yet because the global stiffness matrix of grid model adopts the upper triangular matrix storage, so each non-zero submatrices on diagonal line also only needs the upper triangle submatrix of storage.Be the dimension of the problem of finding the solution if establish m, the nonzero element of diagonal line non-zero submatrices is m (m+1)/2, and off-diagonal non-zero submatrices data are m
2, can determine thus the quantity of the nonzero element of grid model global stiffness matrix upper triangular matrix.
The grid model that consists of take eight node hexahedral elements by limited quantity is illustrated above-mentioned computation process as example.
For the grid model that is comprised of eight node hexahedral elements, its " broad sense adjacent node to " is divided into three types, and the one, limit adjacent node pair; The 2nd, face adjacent node pair; The 3rd, body adjacent node pair.The limit adjacent node is to being determined by the quantity of adjacent node line in the unit, corresponding one of every line " broad sense adjacent node to "; The right quantity of face adjacent node determines by the quantity of elemental area, and each face has two diagonal line, corresponding one of every face diagonal " broad sense adjacent node to "; For body adjacent node pair, there are four body diagonals each unit, corresponding one of every body diagonal " broad sense adjacent node to ".If N
1Be number of nodes (being the node of node and self generation to quantity), N
2Be the quantity (being the right quantity of limit adjacent node) of adjacent node line in grid model, N
3(each face has two face adjacent nodes pair, 2N for the face number of unit in grid model
3The quantity that the presentation surface adjacent node is right), N
4(there are two 4 individual adjacent nodes pair, 4N in each unit for the quantity of unit in grid model
4Expression body adjacent node right quantity), the quantity of grid model " broad sense adjacent node to " is:
S
SM＝N
1+N
2+2×N
3+4×N
4
N wherein
1Be the right quantity of node that node self consists of, the type node to the non-zero submatrices line number of correspondence be listed as number unanimously, be positioned on diagonal line, so N
1Quantity for the diagonal line non-zero submatrices.
For three-dimensional problem (m=3), the data that need to store of each diagonal line non-zero submatrices are m (m+1)/2=6, and the data that the off-diagonal non-zero submatrices need to be stored are m
2=9.Therefore the quantity of the global stiffness matrix upper triangular matrix nonzero element of grid model is:
S
NZ＝6×N
1+9×(N
2+2×N
3+4×N
4)
The below is described further above-mentioned computation process as an example of two concrete grid models example.
For the grid model that is consisted of by a hexahedral element shown in accompanying drawing 4, N
1=8, N
2=12, N
3=6, N
4=1, when adopting the upper triangular matrix storage, its non-zero submatrices quantity is:
S
SM＝N
1+（N
2+2×N
3+4×N
4）=8+(12+2×6+1×4)=36
The global stiffness matrix upper triangular matrix nonzero element quantity that grid model shown in Figure 4 need to be preserved is:
S
NZ＝6×N
1+9×(N
2+2×N
3+4×N
4)=6×8+9×28=300
In like manner, for the grid model that is consisted of by two hexahedral elements shown in accompanying drawing 5, N
1=12, N
2=20, N
3=11, N
4=2.When adopting the upper triangular matrix storage, its non-zero submatrices quantity is:
S
SM＝N
1+（N
2+2×N
3+4×N
4）=12+(20+2×11+2×4)=62
The global stiffness matrix upper triangular matrix nonzero element quantity that grid model shown in Figure 5 need to be preserved is:
S
NZ＝6×N
1+9×(N
2+2×N
3+4×N
4)=6×12+9×50=522
Determined the quantity of non-zero submatrices and nonzero element, the storage space of the global stiffness matrix upper triangular matrix of distribution network lattice model exactly, and then realize the storage of the global stiffness matrix upper triangular matrix of grid model.
The below illustrates the compression storage advantage of global stiffness matrix of the present invention take the global stiffness matrix of the grid model that is made of two hexahedral elements shown in Figure 5 as example.
The data length of the J array in global stiffness matrix storage method involved in the present invention is the quantity of global stiffness matrix non-zero submatrices, i.e. S
SM=62, and the data length of K array is element number 12, the global stiffness matrix stores data volume of this grid model is 62+12=74.If adopt CSR(Compressed Sparse Row format commonly used at present) the form storage means, the data length of its J array is the quantity of all nonzero elements, i.e. S
NZThe data length of=522, K array is the dimension of global stiffness equation, is 12 * 3=36, and therefore, total amount of data is 558.If only compare two auxiliary array J and K, global stiffness matrix storage method involved in the present invention is compared with the CSR method, and the storage space of its saving is: (558-74)/558=87%.If floating type array A adopts single precision floating datum to preserve, global stiffness matrix storage method involved in the present invention is compared with the CSR method, and the space that the global stiffness matrix stores is saved is: 1-(74+552 * 2)/(558+552 * 2)=29%.As seen, for example shown in Figure 5, the storage space that the present invention saves auxiliary array J and K reaches 87%, saves global stiffness matrix stores space 29%, the space-saving successful.
Below by a concrete forging example, beneficial effect of the present invention is described further.
Accompanying drawing 6a-Fig. 6 d is the forging finite element result that a certain moment adopts different dividing elements quantity to obtain in forging process, wherein the unit of all grid models is all eight node hexahedral elements, forging grid model in Fig. 6 a contains 1565 unit and 2045 nodes, forging grid model in Fig. 6 b contains 3284 unit and 4119 nodes, forging grid model in Fig. 6 c contains 5339 unit and 6557 nodes, and the forging grid model in Fig. 6 d contains 8530 unit and 10217 nodes.The below compares the algorithm of compression storage iterative algorithm of the present invention with half band storage direct solution from EMS memory occupation and counting yield two aspects.The computing machine that wherein adopts is: CPU:Intel(R) Core(TM) 2, and Duo, CPU E8400,3.00GHz, internal memory: 3.25, operating system: Windows XP.
Aspect EMS memory occupation, the element of the global stiffness equation of forging grid model adopts single precision floating datum to preserve, and namely each data takies 4 bytes; Auxiliary array adopts the shaping data to preserve, and namely each data takies 2 bytes.Under different nodes and element number, the shared storage space of forging grid model global stiffness matrix of storing is in different ways compared, result is as shown in table 1.
The comparison in table 1 forging grid model global stiffness matrix stores space under different storage meanss
By as seen from Table 1, the shared storage space of compression and storage method of the present invention is far smaller than partly the storage space with storage means, when the element number of forging grid model is 1565, and conserve storage 93.2%, when the element number of forging grid model is 8530, conserve storage 97.4%.As seen, with the increase of grid model node or element number, it is increasing that storage space is saved ratio.This is main because compression and storage method of the present invention has only been preserved the global stiffness nonzeros, in the global stiffness matrix, the quantity of every delegation nonzero element is only relevant with the quantity of the broad sense adjacent node of this row corresponding node, and irrelevant with the number of nodes of integral grid model.In other words, for compression and storage method of the present invention, storage space is only relevant with the line number increase of global stiffness matrix with the proportional relation of number of nodes, is linear trends of change.And for partly being with storage means, the storage space of every delegation is determined by the half-band width of global stiffness matrix, increases with number of nodes, and half-band width is also increase thereupon generally, and the data that the line number of global stiffness matrix and every row need to be stored all increase to some extent, thereby present the para-curve variation tendency.Therefore, along with the increase of number of nodes, partly being quadratic function with the required storage space of storage means increases, and the required storage space of compression and storage method of the present invention is the linear function increase, and its space-saving advantage is more remarkable.
Aspect counting yield, herein only take the global stiffness equation single calculation efficient of finding the solution a certain moment in forging process as example, half-and-half band storage direct solution and two kinds of methods of compression of the present invention storage iterative compare, and result is as shown in subordinate list 2.
Table 2 is that the solution efficiency of grid model when different storage means of different units and number of nodes compares
As can be seen from the table, when the element number of forging grid model was 1565, counting yield improved 69.4%, when the element number of forging grid model is 8530, improved counting yield 88.3%.As seen, the efficient of compression storage iterative of the present invention is far above the efficient of half band storage direct solution, and number of nodes is more, and solution efficiency improves more obvious.Therefore, the method that proposes of the present invention is particularly suitable for the grid model unit and interstitial content is large, the finite element numerical analysis of the forging of finding the solution of iterating, extruding, the metal volume Problems in forming such as rolling.
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
2-(I-1)×m×(m-1)/2
Wherein, m
2Be 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.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 201110130073 CN102184298B (en) | 2011-05-18 | 2011-05-18 | Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN 201110130073 CN102184298B (en) | 2011-05-18 | 2011-05-18 | Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN102184298A CN102184298A (en) | 2011-09-14 |
| CN102184298B true CN102184298B (en) | 2013-06-19 |
Family
ID=44570475
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN 201110130073 Active CN102184298B (en) | 2011-05-18 | 2011-05-18 | Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN102184298B (en) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10713400B2 (en) | 2017-04-23 | 2020-07-14 | Cmlabs Simulations Inc. | System and method for executing a simulation of a constrained multi-body system |
| US10467359B2 (en) * | 2017-08-10 | 2019-11-05 | Livermore Software Technology Corp. | Special-purpose programmed computer for numerical simulation of a metal forming process having a predefined load path with corresponding mesh adjustment scheme |
| CN107908844A (en) * | 2017-11-07 | 2018-04-13 | 武汉科技大学 | Metamorphic mechanisms kinematics model auto-creating method |
| CN108724817B (en) * | 2018-05-04 | 2020-03-17 | 河南工业大学 | Intelligent manufacturing method of thermoplastic composite material metal sandwich plate product |
| CN109948279B (en) * | 2019-03-29 | 2022-12-20 | 江苏精研科技股份有限公司 | Simulation design method for shaping metal piece |
| CN111008497B (en) * | 2019-12-06 | 2022-11-11 | 集美大学 | Method for generating finite element total stiffness matrix and terminal |
| CN111965694B (en) * | 2020-08-03 | 2024-01-30 | 中国石油天然气集团有限公司 | Method and device for determining position of physical point of earthquake and earthquake observation system |
| CN112434451B (en) * | 2020-10-28 | 2024-04-02 | 西安交通大学 | Finite element analysis method based on block parallel computing |
| CN115081297B (en) * | 2022-08-23 | 2022-11-04 | 天津大学 | Overall stiffness matrix solving method in composite material plane elastic finite element analysis |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1139242A (en) * | 1996-02-02 | 1997-01-01 | 曾攀 | Numerical composite unit method for obtaining static and dynamic characteristics of engineering structure |
| CN1386201A (en) * | 2000-06-29 | 2002-12-18 | 目标储油层公司 | Method and system for solving finite element models using multi-phase physics |
| CN1609351A (en) * | 2004-11-24 | 2005-04-27 | 东南大学 | Set-up method for generalized Bingham soft soil rheological deformation analogue body |
-
2011
- 2011-05-18 CN CN 201110130073 patent/CN102184298B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1139242A (en) * | 1996-02-02 | 1997-01-01 | 曾攀 | Numerical composite unit method for obtaining static and dynamic characteristics of engineering structure |
| CN1386201A (en) * | 2000-06-29 | 2002-12-18 | 目标储油层公司 | Method and system for solving finite element models using multi-phase physics |
| CN1609351A (en) * | 2004-11-24 | 2005-04-27 | 东南大学 | Set-up method for generalized Bingham soft soil rheological deformation analogue body |
Also Published As
| Publication number | Publication date |
|---|---|
| CN102184298A (en) | 2011-09-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN102184298B (en) | Method for storing and generating stiffness matrix for finite-element analysis in metal bulk plastic forming | |
| CN110110413B (en) | Structural topology optimization method based on material field reduction progression expansion | |
| CN111709171B (en) | Isogeometric solving and heat dissipation topology generation method for heat flow strong coupling problem | |
| CN108446445B (en) | Composite material wing optimization design method based on aerodynamic reduced order model | |
| CN103324850B (en) | Twice polycondensation parallel method of finite element two-stage subregion based on multifile stream | |
| Yin et al. | On the ensemble of metamodels with multiple regional optimized weight factors | |
| CN108846212B (en) | Rigid frame pile internal force and displacement design calculation method | |
| CN111709085A (en) | Topological optimization design method for constrained damping sheet structure | |
| CN107491585B (en) | Structural topology optimization design method taking random displacement response variance as target | |
| Gao et al. | An exact block-based reanalysis method for local modifications | |
| CN104573281A (en) | Complex space curved surface thin plate forming die face designing method taking springback compensation | |
| CN105574255A (en) | Simple implementation method for predicting periodical composite material thermal conductivity coefficient in gradual and homogeneous manner | |
| CN102063525B (en) | Method for generating practical multidisciplinary design optimization model automatically | |
| Tang et al. | Lattice-skin structures design with orientation optimization | |
| CN114347029A (en) | Model order reduction method for rapid simulation of pneumatic soft robot | |
| Idczak et al. | Minimization of Poisson’s ratio in anti-tetra-chiral two-phase structure | |
| CN111680441B (en) | Gradient lattice sandwich board structure suitable for thermal working condition | |
| CN111274624B (en) | Multi-working-condition special-shaped node topology optimization design method based on RBF proxy model | |
| CN101114308A (en) | Three kinds variables range B sample strip wavelet beam unit constructing method | |
| CN112163385A (en) | Parallel non-grid method and system for solving strong heat fluid-solid coupling problem | |
| CN103309845A (en) | Partitioning solving method for linear simultaneous equations of dynamic simulation of power system | |
| CN116187074A (en) | Multi-scale topological optimization method of anisotropic periodic structure material based on isogeometry | |
| CN112966421B (en) | Calculation method for calculating buckling load factor and corresponding buckling shape of sheet structure by using p-type finite element method | |
| Shiye et al. | Topology Optimization Design of 3D Continuum Structure with Reserved Hole Based on Variable Density Method. | |
| Nigro et al. | Matrix‐free modified extended BDF applied to the discontinuous Galerkin solution of unsteady compressible viscous flows |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant |