Novel stable node integral calculation method based on hourglass control

Abstract

The invention relates to a novel stable node integral calculation method based on hourglass control, and the method comprises the steps: dividing a problem domain into a triangular or tetrahedral mesh, and dividing a mesh integral domain into a plurality of non-overlapping node integral domains based on nodes in the triangular or tetrahedral mesh; performing node domain smoothing on the node integral domain based on the triangular / tetrahedron strain component of the node to obtain a smooth domain; hourglass control between a node and an adjacent node is constructed in the smooth domain, and hourglass force between the node and the adjacent node is calculated; and superposing the hourglass force borne by the node with the original internal force calculated by node integration to obtain the internal force of a stable node integration method, and solving a discrete motion equation in explicit dynamics based on the internal force of the stable node integration method to obtain a dynamic response. According to the method, the calculation precision based on the unstructured grid is greatly improved, and the advantages that the number of integration points of node integration is small and the grid distortion resistance is high are effectively kept.

Description

Technical Field

[0001] This application relates to the technical field of high-performance numerical simulation analysis, and particularly to a new stable nodal integration calculation method based on hourglass control. Background Art

[0002] At present, although the space-time stable nodal integration algorithm constructed by introducing strain gradient into the energy equation solves the problem of numerical instability of nodal integration, multiple integration points located on the coordinate axes are artificially introduced during the explicit dynamics calculation process, resulting in an increase in the amount of calculation and the inability to reflect the efficiency advantage of nodal integration, which limits its application and popularization in dynamic large-deformation explicit dynamics calculation. Therefore, how to ensure the numerical stability of the stable nodal integration algorithm on the premise of reducing the number of artificial integration points of the stable nodal integration algorithm and maintaining the single-integration-point calculation feature of nodal integration, and how to construct a single-integration-point stabilization term are the key scientific problems that need to be solved urgently. Summary of the Invention

[0003] Based on this, it is necessary to provide a new stable nodal integration calculation method based on hourglass control, and this method includes: S1: Divide the problem domain into triangular or tetrahedral meshes, and divide the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the triangular or tetrahedral meshes; S2: Smooth the nodal integration domain based on the triangular / tetrahedral strain components of the nodes to obtain a smoothed domain; S3: Construct hourglass control between the node and its adjacent nodes in the smoothed domain, and calculate the hourglass force between the node and its adjacent nodes; S4: Superimpose the hourglass force received by the node and the original internal force calculated by nodal integration to obtain the internal force of the stable nodal integration method, and solve the discrete motion equation in explicit dynamics based on the internal force of the stable nodal integration method to obtain the dynamic response.

[0004] Preferably, in S1, dividing the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the triangular mesh includes: for the triangular mesh divided for two-dimensional problems, sequentially connect the midpoints of the surrounding edges of a node and the midpoints of the surrounding triangular elements to obtain multiple non-overlapping nodal integration domains centered on the node.

[0005] Preferably, in S1, dividing the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the tetrahedral mesh includes: for the tetrahedral mesh divided for three-dimensional problems, sequentially connect the midpoints of the surrounding edges of a node, the midpoints of the surrounding triangles, and the centroid of the surrounding tetrahedral elements to obtain multiple non-overlapping nodal integration domains centered on the node.

[0006] Preferably, in S2, the node domain smoothing of the triangle / tetrahedron strain components based on nodes includes: ; where represents the derivative of the smoothing shape function of the node integration domain of the I-th node; represents the area / volume of the I-th node; represents the derivative of the shape function of the node integration domain of the I-th node; represents the gradient; represents the outer normal of the surrounding edge / surrounding face of the I-th node; represents the surrounding edge / surrounding face of the I-th node.

[0007] Preferably, S3 includes: Calculating the theoretical relative position relationship between a node and its adjacent nodes based on the deformation gradient increment between any node in the smoothing domain and its adjacent nodes; Calculating the actual relative position relationship between a node and its adjacent nodes based on the positions of the nodes and adjacent nodes in the smoothing domain; Calculating the relative position deviation between the theoretical relative position relationship and the actual relative position relationship; Calculating the elongation of the node to the adjacent node based on the relative position deviation; Taking the derivative of the potential energy function of the virtual spring with respect to the length to obtain the virtual spring force; Calculating the internal force between the node and the adjacent node based on the elongation and the virtual spring force; Symmetrizing the internal force between the symmetric node and the adjacent node to obtain the hourglass force between the node and the adjacent node.

[0008] Preferably, the calculation of the elongation of the node to the adjacent node based on the relative position deviation includes: ; where represents the elongation of the i -th node to its j -th adjacent node; represents the relative position deviation on the side of the -th time step of the i -th node; represents the direction vector between the -th time step of the node i and j node.

[0009] Preferably, the calculation formula of the hourglass force is: ; where Represents the hourglass force between the i th node and all its adjacent nodes; Represents the hourglass force between the i th node and its j th adjacent node; Represents the number of adjacent nodes of the i th node; Represents the actual relative position relationship between the t th node and its i th adjacent node at the j th time step; Represents the L1 norm; Represents the elastic coefficient of the i th node; Represents the elongation of the i th node to its j th adjacent node; Represents the elastic coefficient of the j th node; Represents the elongation of the i th node's j th adjacent node to the i-th node; Represents the hourglass coefficient; Represents the smooth domain volume / area of node i ; Represents the smooth domain volume / area of node j ;

[0010] Preferably, the original internal force calculation formula for the node is: ; ; ; Among them, Represents the original internal force of the I-th node; Represents the smooth domain of node k ; Represents the partial derivative of the smooth shape function of node I in the j direction; Represents the internal force between node j and node i ; k Represents node k ; Represents the deformation of node k; Represents the smooth domain volume / area of node k ; Represents the partial derivative of the shape function of node I in the j direction; Represents the determinant of the deformation gradient; Represents a node j and the node k The deformation gradient therebetween; Represents a node k and the node i The stress therebetween.

[0011] Preferably, the discrete motion equation in explicit dynamics for solving internal forces based on the stabilized nodal integration method includes: Calculate the external forces on the nodes by the stabilized nodal integration method; Subtract the internal forces of the stabilized nodal integration method on the nodes from the external forces of the stabilized nodal integration method on the nodes, and use the difference as the external force vector in the discrete motion equation; The external force vector is the product of the lumped mass matrix and the acceleration vector, and the dynamic response is obtained.

[0012] Preferably, the calculation formula for the external forces on the nodes by the stabilized nodal integration method is: ; Wherein, Represents the external force on the I-th node by the stabilized nodal integration method; Represents the nodal integration domain of the I-th node; Represents the nodal smoothing domain; b represents the body force; Represents the boundary of the problem domain; Represents the surface force.

[0013] Beneficial effects: This method constructs a stabilized nodal integration method based on triangular / tetrahedral unstructured meshes, which can adapt to complex geometric boundaries and greatly improve the calculation accuracy based on unstructured meshes; at the same time, the calculation stability of nodal integration is improved based on the calculation of hourglass forces between nodes and adjacent nodes, effectively maintaining the advantages of fewer integration points and strong anti-mesh distortion ability in nodal integration, and it is a very efficient and stable explicit dynamic large deformation calculation method. Brief Description of the Drawings

[0014] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present application. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0015] Figure 1 Is the flowchart of the novel stabilized nodal integration calculation method based on hourglass control in the embodiments of the present application.

[0016] Figure 2 Is the schematic diagram of the nodal integration domain in the embodiments of the present application.

[0017] Figure 3 This is a schematic diagram of hourglass control between nodes and adjacent nodes in the embodiments of the present application. Detailed implementation manners

[0018] To make the above objects, features, and advantages of the present application more obvious and understandable, the following will describe the detailed implementation manners of the present application with reference to the accompanying drawings. Many specific details are set forth in the following description to fully understand the present application. However, the present application can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the connotation of the present application. Therefore, the present application is not limited by the specific embodiments disclosed below.

[0019] In addition, the terms "first" and "second" are only used for descriptive purposes and cannot be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include at least one of such features. In the description of the present application, "a plurality of" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0020] As Figure 1 shown, this embodiment provides a new type of stable node integration calculation method based on hourglass control, and the method includes: S1: Divide the problem domain into triangular or tetrahedral meshes, and divide the mesh integration domain into multiple non-overlapping node integration domains based on the nodes in the triangular or tetrahedral meshes.

[0021] Specifically, in this step: Dividing the mesh integration domain into multiple non-overlapping node integration domains based on the nodes in the triangular mesh includes: for the triangular mesh divided for two-dimensional problems, successively connecting the midpoints of the surrounding edges of a node and the midpoints of the surrounding triangular elements to obtain multiple non-overlapping node integration domains centered on the node.

[0022] Dividing the mesh integration domain into multiple non-overlapping node integration domains based on the nodes in the tetrahedral mesh includes: for the tetrahedral mesh divided for three-dimensional problems, successively connecting the midpoints of the surrounding edges of a node, the midpoints of the surrounding triangles, and the centroid of the surrounding tetrahedral elements to obtain multiple non-overlapping node integration domains centered on the node.

[0023] The divided node integration domains are as Figure 2 shown. The node integration domains constructed in this way can avoid regional duplication or omission during integration. It should be noted that the interpolation and integration of the finite element method are carried out within the element, and the interpolation domain and the integration domain are kept consistent, while the interpolation and integration of node integration are completed within the element and the constructed node smoothing domain respectively.

[0024] S2: Smooth the nodal integration domain based on the triangular / tetrahedral strain components of the nodes to obtain a smoothed domain.

[0025] Specifically, the calculation formula for nodal domain smoothing is: ; where represents the derivative of the smoothed shape function of the nodal integration domain of the I-th node; represents the area / volume of the I-th node; represents the derivative of the shape function of the nodal integration domain of the I-th node; represents the gradient; represents the outer normal of the surrounding edge / surface of the I-th node; represents the surrounding edge / surface of the I-th node. This method is equivalent to reconstructing the gradient based on the nodes within the nodal integration domain, which can improve the drawback that the results based on triangular / tetrahedral elements are too "rigid".

[0026] S3: Construct hourglass control between a node and its adjacent nodes in the smoothed domain and calculate the hourglass force between the node and its adjacent nodes.

[0027] When calculating the strain / deformation gradient by nodal integration, it can be simply understood as homogenizing the strain / deformation gradient of the surrounding elements of the node. This deformation method is not unique, that is, different deformation methods may result in the same strain / deformation gradient. Just a slight perturbation will cause the mesh to randomly transform between different deformation methods, which is very similar to the hourglass phenomenon caused by using reduced integration for quadrilateral / hexahedron elements. This non-uniqueness will lead to false zero-energy modes, resulting in very unstable performance of nodal integration in explicit dynamic calculations and making it impossible to be applied to actual engineering problems. Nodal integration adopts single-point integration, and relative movement is possible between a node and its adjacent nodes. In this embodiment, artificial stiffness between a node and its adjacent nodes is introduced to limit the transition between different deformation methods, thereby overcoming the hourglass phenomenon, as Figure 3 shown.

[0028] Specifically, this step includes: Based on the deformation gradient increment between any node in the smoothed domain and its adjacent nodes, calculate the theoretical relative position relationship between the node and its adjacent nodes. The calculation formula is: ; where represents the theoretical relative position relationship with the I-th node as the main one at the time step; i represents the I-th at the t time step and the iThe actual relative position relationship between a node and its j th adjacent node; represents the time step increment; represents the time step, and the incremental deformation gradient between the i th node and its j th adjacent node.

[0029] Based on the positions of the nodes within the smooth domain and the positions of their adjacent nodes, calculate the actual relative position relationship between a node and its adjacent nodes. The calculation formula is: ; where represents the actual relative position relationship between the th node and its i th adjacent node at the j time step; represents the position of the th adjacent node of the i th node at the j time step; represents the position of the th node at the i time step.

[0030] Calculate the relative position deviation between the theoretical relative position relationship and the actual relative position relationship. The calculation formula is: ; where represents the relative position deviation on the side of the th node at the i time step; represents the theoretical relative position relationship with the th node as the main node at the i time step; represents the actual relative position relationship between the th node and its i th adjacent node at the j time step.

[0031] This operation is equivalent to adding a spring between the node and its adjacent node, applying a certain stiffness to limit the deformation. To ensure that the system still maintains a linear displacement field after adding the spring, calculate the elongation of the node to its adjacent node based on the relative position deviation. The calculation formula is: ; ; where represents the elongation of the i th node to itsj The elongation of adjacent nodes; Denote the relative position deviation at the i th time step on the side of the Denote the relative position deviation at the i th time step on the side of the j th node and its adjacent node; Denote the direction vector between node i and j node at the

[0032] Derive the virtual spring force by differentiating the potential energy function of the virtual spring with respect to length. The calculation formula is: ; ; where f denotes the spring force; denotes the elastic potential energy; L denotes the length variable; E denotes the elastic coefficient; denotes the deformation of the virtual spring; denotes the initial length of the spring.

[0033] Calculate the internal force between the node and its adjacent nodes based on the elongation and the virtual spring force. The calculation formula is: ; where denotes the internal force on the side of the i th node; E denotes the elastic coefficient; denotes the t relative position relationship between the i th node and its j adjacent node at the th time step; denotes the i L1 norm; j denotes the elongation of the

[0034] Symmetrize the internal forces between the symmetric nodes and their adjacent nodes to obtain the hourglass force between the node and its adjacent nodes. The calculation formula is: ; where denotes the hourglass force between the i th node and all its adjacent nodes; denotes thei The hourglass force between a node and its j th adjacent node; Denote the number of adjacent nodes of the i th node; Denote the t th time step of the actual relative position relationship between the i th node and its j th adjacent node; Denote the L1 norm; Denote the i th node's elastic coefficient; Denote the i th node's elongation to its j th adjacent node; Denote the j th node's elastic coefficient; Denote the i th node's j th adjacent node's elongation to the i th node; Denote the hourglass coefficient; Denote the smooth domain volume / area of node i ; Denote the smooth domain volume / area of node j ;

[0035] S4: Superimpose the hourglass force on the node and the original internal force calculated by node integration to obtain the internal force of the stable node integration method. Solve the discrete motion equation in explicit dynamics based on the internal force of the stable node integration method to obtain the product of the lumped mass matrix and the acceleration vector, and obtain the dynamic response.

[0036] In this embodiment, the calculation formula for the original internal force of the node is: ; ; ; Among them, Denote the original internal force of the Ith node; Denote the smooth domain of node k ; Denote the partial derivative of the smooth shape function of node I in the j direction; Denote the internal force between node j and node i ; k Denote node k ; Denote the deformation of node k; Denote nodek Smooth domain volume / area; Denote the partial derivative of the shape function of node I with respect to j the direction; Denote the determinant of the deformation gradient; Denote node j and node k the deformation gradient between them; Denote node k and node i the stress between them.

[0037] Due to the characteristics of node integration, it can be calculated by the following formula: ; where, Denote the i th node and its j th adjacent node elongation; Denote the partial derivative; Denote the node integration domain of the I-th node; Denote the number of influencing nodes around the node; Denote node j deformation; Denote the node displacement vector.

[0038] Furthermore, the discrete motion equation in explicit dynamics for internal force solution based on the stabilized node integration method includes: Calculate the external force of the stabilized node integration method acting on the node; Subtract the external force of the stabilized node integration method acting on the node from the internal force of the stabilized node integration method, and use the difference as the external force vector in the discrete motion equation; The external force vector is the product of the lumped mass matrix and the acceleration vector, and the dynamic response is obtained.

[0039] Even further, the calculation formula for the external force of the stabilized node integration method acting on the node is: ; where, Denote the external force of the stabilized node integration method acting on the I-th node; Denote the node integration domain of the I-th node; Denote the node smooth domain; b denotes the body force; Denote the problem domain boundary; Denote the surface force.

[0040] The entire process is calculated based on the original triangular / tetrahedral mesh without modifying the existing finite element calculation framework. It can not only make full use of the advantage of fewer integration points in nodal integration to improve the calculation efficiency, but also has strong anti-mesh distortion ability, significantly improving the calculation accuracy under low-quality meshes, and is very suitable for explicit analysis involving large material deformations such as metal volume forming and contact collision.

[0041] The new stable nodal integration calculation method based on hourglass control provided in this embodiment has the following beneficial effects: 1. Traditional linear elements based on unstructured meshes such as triangles / tetrahedrons are difficult to ensure calculation accuracy, while structured meshes based on quadrilaterals / hexahedrons do not have good discretization adaptability. The calculation method provided in this embodiment constructs stable nodal integration elements based on unstructured meshes. On the premise of being able to well adapt to complex geometric boundaries, it greatly improves the calculation accuracy, and to a certain extent solves the problem of the contradiction between calculation accuracy and calculation efficiency, providing a new idea for solving practical engineering problems with complex boundaries.

[0042] 2. Construct the hourglass force based on the relative deformation between the nodes in the nodal integration domain and the adjacent nodes, and develop a stable nodal integration method with the characteristics of single integration point, reducing the number of artificial integration points in the original stable nodal integration method. It can not only maintain the high accuracy and high anti-mesh distortion performance of the stable nodal integration method, but also significantly improve the calculation efficiency, providing a new method for the rapid prediction of dynamic large deformation explicit dynamic responses.

[0043] The technical features of the above-described embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in this specification.

[0044] The above-described embodiments only represent several implementation manners of the present application, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the patent application. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several deformations and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent of the present application should be subject to the appended claims.

Claims

1. A novel stable nodal integration calculation method based on hourglass control, characterized in that, Including: S1: Divide the problem domain into triangular or tetrahedral meshes, and divide the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the triangular or tetrahedral meshes; S2: Smooth the nodal integration domain based on the triangular / tetrahedral strain components of the nodes to obtain a smoothed domain; S3: Construct hourglass control between a node and its adjacent nodes in the smoothed domain, and calculate the hourglass force between the node and its adjacent nodes; S4: Superimpose the hourglass force received by the node and the original internal force calculated by nodal integration to obtain the internal force of the stabilized nodal integration method, and solve the discrete motion equation in explicit dynamics based on the internal force of the stabilized nodal integration method to obtain the dynamic response.

2. The novel stable nodal integration calculation method based on hourglass control according to claim 1, characterized in that, In S1, dividing the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the triangular mesh includes: for the triangular mesh divided for two-dimensional problems, successively connect the midpoints of the surrounding edges of a node and the midpoints of the surrounding triangular elements to obtain multiple non-overlapping nodal integration domains centered on the node.

3. The novel stable nodal integration calculation method based on hourglass control according to claim 1, characterized in that, In S1, dividing the mesh integration domain into multiple non-overlapping nodal integration domains based on the nodes in the tetrahedral mesh includes: for the tetrahedral mesh divided for three-dimensional problems, successively connect the midpoints of the surrounding edges of a node, the midpoints of the surrounding triangles, and the centroid of the surrounding tetrahedral elements to obtain multiple non-overlapping nodal integration domains centered on the node.

4. The novel stable nodal integration calculation method based on hourglass control according to claim 1, characterized in that, In S2, the smoothing of the nodal integration domain based on the triangular / tetrahedral strain components of the nodes includes: ； Among them, represents the derivative of the smooth shape function of the nodal integration domain of the $i$-th node; represents the area / volume of the $i$-th node; represents the derivative of the shape function of the nodal integration domain of the $i$-th node; represents the gradient; represents the outer normal of the surrounding edge / surrounding face of the $i$-th node; represents the surrounding edge / surrounding face of the $i$-th node.

5. The novel stable nodal integration calculation method based on hourglass control according to claim 1, characterized in that, S3 includes: Calculate the theoretical relative position relationship between a node and its adjacent nodes based on the deformation gradient increment between any node and its adjacent nodes in the smoothed domain; Calculate the actual relative position relationship between a node and its adjacent nodes based on the positions of the node and its adjacent nodes in the smoothed domain; Calculate the relative position deviation between the theoretical relative position relationship and the actual relative position relationship; Calculate the elongation from the node to its adjacent nodes based on the relative position deviation; Take the derivative of the potential energy function of the virtual spring with respect to the length to obtain the virtual spring force; Calculate the internal force between the node and its adjacent nodes based on the elongation and the virtual spring force; Symmetrize the internal force between the symmetric node and its adjacent nodes to obtain the hourglass force between the node and its adjacent nodes.

6. The novel stable nodal integration calculation method based on hourglass control according to claim 5, characterized in that, The calculation of the elongation from the node to its adjacent nodes based on the relative position deviation includes: ； Among them, represents the elongation from the i -th node to its j -th adjacent node; represents the relative position deviation at the -th time step on the side of the i -th node; represents the direction vector between node and node i and j at the 7. The novel stable nodal integration calculation method based on hourglass control according to claim 5, characterized in that, The calculation formula for the hourglass force is: ； Among them, represents the hourglass force between the i th node and all its adjacent nodes; represents the hourglass force between the i th node and its j th adjacent node; represents the number of adjacent nodes of the i th node; represents the actual relative position relationship between the t th node and its i th adjacent node at the j th time step; represents the L1 norm; represents the elastic coefficient of the i th node; represents the elongation of the i th node to its j th adjacent node; represents the elastic coefficient of the j th node; represents the elongation of the i th adjacent node of the j th node to the th node; represents the hourglass coefficient; i represents the smooth domain volume / area of node represents the smooth domain volume / area of node j ​ 8. The novel stable nodal integration calculation method based on hourglass control according to claim 1, characterized in that, The calculation formula for the original internal force of the node is: ； ； ； Among them, represents the original internal force of the $I$-th node; represents the k smooth domain of the node; represents the partial derivative of the smooth shape function of node $I$ with respect to the j direction; represents the j internal force between node i and node k represents the k node; represents the deformation of node $k$; represents the k volume / area of the smooth domain of the node; represents the partial derivative of the shape function of node $I$ with respect to the j direction; represents the determinant of the deformation gradient; represents the j deformation gradient between node k and node represents the k stress between node i and node.

9. The novel stable node integration calculation method based on hourglass control according to claim 1, characterized in that, The solution of the discrete motion equation in explicit dynamics based on the internal force of the stabilized nodal integration method includes: Calculate the external force of the stabilized nodal integration method received by the node; Subtract the external force of the stabilized nodal integration method received by the node from the internal force of the stabilized nodal integration method, and use the difference as the external force vector in the discrete motion equation; The external force vector is the product of the lumped mass matrix and the acceleration vector, and the dynamic response is obtained.

10. The novel stable node integration calculation method based on hourglass control according to claim 1, characterized in that, The calculation formula for the external force of the stabilized nodal integration method received by the node is: ； Among them, represents the external force of the stabilized nodal integration method acting on the $I$-th node; represents the nodal integration domain of the $I$-th node; represents the nodal smoothing domain; $b$ represents the body force; represents the boundary of the problem domain; represents the surface force.

Images