Finite element calculation method of nearly incompressible material and computer equipment

Abstract

The invention discloses a finite element calculation method of a nearly incompressible material and computer equipment, and the finite element calculation method comprises the following steps: constructing a geometric model, endowing parameters of the nearly incompressible material to the geometric model, endowing boundary conditions and loads to the geometric model, and obtaining a finite element parameter model of the nearly incompressible material; performing mesh generation on the finite element parameter model according to a pre-selected unit type to obtain a finite element model; establishing a unit stiffness matrix of the units; assembling the unit stiffness matrix to form a global stiffness matrix; discretizing the load to generate a global node force vector, and calculating to obtain a node displacement vector; and performing post-processing by utilizing the node displacement vector to obtain displacement, stress and strain of the finite element model. According to the method, the volume self-locking characteristic of the unit is eliminated by correcting the part, related to the volume, in the shape function derivative matrix B, so that the volume self-locking problem of the nearly incompressible material in the finite element calculation process can be well solved.

Description

Technical Field

[0001] The present invention relates to the technical field of engineering digital simulation calculation, and more particularly to a finite element calculation method and computer equipment for nearly incompressible materials. Background Art

[0002] Nearly incompressible materials (Poisson's ratio close to 0.5) are materials whose volume changes very little when subjected to stress. In practical engineering and scientific research, it is often necessary to study the properties of nearly incompressible materials, such as in bionics research, soft matter engineering, and the design of rubber parts. The model of a nearly incompressible material under linear elastic material is:

[0003] σ=2με+λε v I,

[0004] In the above formula, σ is the stress tensor, ε is the strain tensor, I is the unit tensor, and the shear modulus is Lame coefficient Volumetric strain ε v =tr(ε)=ε xx +ε yy +ε zz , E is the elastic modulus. Therefore, in numerical simulations, since the Poisson's ratio v of nearly incompressible materials approaches 0.5, μ is always constant, and λ tends to infinity, the volume strain must satisfy ε v ≈0, the standard displacement element will force ε v = 0 is assigned to each Gaussian integration point, which causes the system volume to self-lock (overconstraint), which has a significant impact on engineering design and material performance analysis and requires special treatment. Summary of the Invention

[0005] The purpose of the present invention is to provide a finite element calculation method and computer equipment for nearly incompressible materials. The present invention eliminates the volume self-locking characteristics of the unit by modifying the volume-related part of the shape function derivative matrix B, thereby better solving the volume self-locking problem that occurs in the finite element calculation process of nearly incompressible materials.

[0006] To achieve the above object, the technical solution of the present invention is as follows:

[0007] A finite element calculation method for nearly incompressible materials comprises the following steps:

[0008] S1: constructing a geometric model, assigning parameters of the nearly incompressible material to the geometric model, assigning boundary conditions and loads to the geometric model, and obtaining a finite element parameter model of the nearly incompressible material;

[0009] S2: Meshing the finite element parameter model according to the pre-selected unit type to obtain a finite element model;

[0010] S3: Establish the element stiffness matrix of the element:

[0011]

[0012] The element stiffness matrix k e is a function of the shape function derivative matrix B and the elasticity matrix D, and the integration domain is the unit, wherein the shape function derivative matrix B of each unit is a matrix composed of m 6×3 sub-matrices.

[0013] B=[B1 B2 … B m ]

[0014] Where m is the number of nodes of the unit type, and the submatrix B of each node is i (i=1~m) is:

[0015]

[0016] Among them B j (j=1~9) is the element shape function N i (i=1~m) for the global coordinate x s (s=1~3), where

[0017]

[0018] B5=B1+B4,

[0019] B7=B2+B6, B9=B3+B8

[0020] The integration domain Ω e is the unit, N i Determined by the unit type;

[0021] S4: Assemble the unit stiffness matrix to form a global stiffness matrix K; discretize the load to generate a global node force vector r, and calculate the node displacement vector u=K -1 r;

[0022] S5: performing post-processing using the node displacement vector to obtain the displacement, stress and strain of the finite element model.

[0023] The present invention also discloses a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the above method.

[0024] The implementation of the present invention will have the following beneficial effects:

[0025] The present invention adopts the average expansion method to correct the volume-related part of the shape function derivative matrix B, decouples the volume deformation from the shear deformation, and eliminates the self-locking problem.

[0026] The present invention adopts the Gaussian integral method to construct a first-order triangular prism unit and a first-order pyramid unit that eliminate the volume self-locking characteristic, thereby improving the calculation accuracy of the deformation of nearly incompressible materials. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0028] in:

[0029] Figure 1 It is a flow chart of a finite element calculation method for nearly incompressible materials provided by an embodiment of the present invention.

[0030] Figure 2 This is a schematic diagram of the 8×4×1 (×2) grid division of a thick-walled cylinder 1 / 4 structure provided by an embodiment of the present invention (each hexagonal unit is further divided into 2 triangular prism units).

[0031] Figure 3 This is a schematic diagram of the 8×4×1 (×6) grid division of a thick-walled cylinder 1 / 4 structure provided by an embodiment of the present invention (each hexagonal unit is further divided into 6 pyramid units).

[0032] Figure 4 It is a structural diagram of a computer device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0033] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0034] refer to Figure 1 The present invention discloses a finite element calculation method for nearly incompressible materials, comprising the following steps:

[0035] S1: Construct a geometric model, assign parameters of the nearly incompressible material to the geometric model, assign boundary conditions and loads to the geometric model, and obtain a finite element parameter model of the nearly incompressible material.

[0036] Nearly incompressible materials are special materials with a Poisson's ratio close to 0.5. This paper develops a finite element calculation method for linear elastic near-incompressible materials. The parameters of linear elastic near-incompressible materials can be characterized solely by the elastic modulus E and Poisson's ratio v, or solely by the shear modulus μ and the Lame coefficient λ. The two characterization methods are consistent and have the following equivalent relationship:

[0037]

[0038] S2: According to the pre-selected unit type, the finite element parameter model is meshed to obtain a finite element model.

[0039] Finite element methods offer a wide variety of element types. Based on the problem to be solved in this paper, the most suitable element type is a three-dimensional element. Within the three-dimensional element, the number of element nodes can be selected. Furthermore, the order of the nodal interpolation function can be categorized to select first-order or higher-order elements.

[0040] In the first embodiment of the present invention, a first-order pyramid unit is selected, and in the second embodiment, a first-order triangular prism unit is selected.

[0041] S3: Establish the element stiffness matrix of the element:

[0042]

[0043] Element stiffness matrix k e It is a function of the shape function derivative matrix B and the elasticity matrix D. The integration domain is the unit, where the shape function derivative matrix B of each unit is a matrix composed of m 6×3 sub-matrices.

[0044] B=[B1 B2 … B m ]

[0045] Where m is the number of nodes of the unit type, and the submatrix B of each node i (i=1~m)

[0046]

[0047] Among them B j (j=1~9) is the element shape function N i (i=1~m) for the global coordinate x s (s=1~3), where

[0048]

[0049] B5=B1+B4,

[0050] B7=B2+B6, B9=B3+B8

[0051] The integration domain Ω e is a unit, N i Determined by the element type.

[0052] Preferably, the element stiffness matrix k is calculated from the shape function derivative matrix B and the elastic matrix D e In the process, numerical integration methods are used for calculation, and the numerical integration methods involved include but are not limited to Gaussian integration, Simpson integration or midpoint integration method.

[0053] Gaussian integral is the most commonly used integration method. Due to its high precision and efficiency, it is suitable for various units and problems. Simpson integral and midpoint integral are mostly used in simple scenarios or teaching verification.

[0054] Preferably, Gaussian integral is used for calculation, and the Gaussian integral expression of the element stiffness matrix is:

[0055]

[0056] Where B is the shape function derivative matrix, D is the elasticity matrix, n is the number of Gaussian integration points, w i is the weight of the i-th Gaussian integration point, B i is the shape function derivative matrix at the i-th Gaussian integration point, B i T Is the transposed matrix, J i is the Jacobian matrix of the i-th Gaussian integration point,

[0057]

[0058] |J i |It's J i The determinant value of , ξ1, ξ2, ξ3 are the three coordinate axes in the unit's parameter coordinate system, and x1, x2, x3 are the three coordinate axes in the unit's global coordinate system.

[0059] In finite element analysis, the isoparametric coordinate system (Isoparametric Coordinate System) is a method that uses coordinate transformation to map elements with complex geometric shapes (such as curved boundary elements) to a simple coordinate system (such as a regular square or cube) for analysis. This method uses shape functions to achieve the mapping between the physical coordinate system and the natural coordinate system, allowing integration and differentiation operations on complex geometries to be converted into standard operations in a simple coordinate system, thereby simplifying the calculation process.

[0060] Under the isoparametric element technology, the parametric coordinate system is a standardized local coordinate system. The present invention adopts a three-dimensional cubic coordinate system with a coordinate range of [-1, 1]. The unit shape is regular, which is convenient for mathematical operations (such as integration and differentiation).

[0061] The transformation between the parametric coordinate system and the global coordinate system is carried out through isoparametric element mapping, that is, the shape function method. Under this method, the differential operator in the global coordinate system needs to be converted into the operator in the parametric coordinate system through the Jacobian matrix, that is,

[0062]

[0063] For nearly incompressible materials, the Poisson's ratio approaches 0.5, which causes some elements in the elastic matrix D to approach 0.

[0064]

[0065] in When υ→0.5, it approaches 0 and volume self-locking occurs. In the present invention, the shape function derivative matrix B is specially processed and the part related to the Lame coefficient λ is corrected. The derivation process is as follows:

[0066] 1) Decompose the elastic matrix and record the part containing only μ as The part containing only λ is recorded as Right now,

[0067]

[0068] 2) The unit stiffness matrix is ​​also decomposed into two parts

[0069]

[0070] 3) Observation From the form of When , only the upper half of the B matrix (rows 1 to 3) is involved in the calculation, that is, the volume deformation does not cause shear deformation.

[0071] 4) In the finite element method, strain is related to displacement through the B matrix, which can be decomposed into deflection (B dev ) and volume (B dil )Two parts:

[0072]

[0073]

[0074] In the above formula, I is the unit tensor, It is the tensor cross product operator symbol.

[0075] 5) In the present invention, the average expansion method is used to i Take the average in the unit domain Ω for correction

[0076]

[0077] 6) The modified shape function derivative matrix B is denoted as

[0078]

[0079] For ease of writing,

[0080]

[0081] B5=B1+B4,

[0082] B7=B2+B6, B9=B3+B8

[0083] In the first embodiment of the present invention, the unit type is a first-order triangular prism unit with a total of 6 nodes. The node shape function of the node is

[0084]

[0085]

[0086] Where ξ1, ξ2, ξ3 are the three coordinate axes in the unit's parameter coordinate system.

[0087] In Example 1, the first-order triangular prism unit uses 6 Gaussian integration points with a weight of 1 / 6. In the parameter coordinate system, the coordinates of the Gaussian integration points are

[0088]

[0089] In the second embodiment of the present invention, the unit type is a first-order pyramid unit, with a total of 5 nodes, and the node shape function of the node is

[0090]

[0091] Where ξ1, ξ2, ξ3 are the three coordinate axes in the unit's parameter coordinate system.

[0092] In the second embodiment of the present invention, the first-order pyramid unit uses 8 Gaussian integration points, all with a weight of 1. In the parameter coordinate system, the coordinates of the Gaussian integration points are

[0093]

[0094] S4: Assemble the unit stiffness matrix to form the global stiffness matrix K; discretize the load to generate the global node force vector r, and calculate the node displacement vector u = K -1 r.

[0095] S5: Post-process the node displacement vectors to obtain the displacement, stress, and strain of the finite element model.

[0096] Next, the beneficial effects of the present invention are described using examples.

[0097] Consider a circular ring of linear elastic incompressible material with a height of 1 mm, an inner diameter of 1 mm, an outer diameter of 4 mm, an elastic modulus of 30 MPa, a Poisson's ratio of 0.5, a fixed displacement in the height direction, and a pressure value of 1 MPa on the inner wall.

[0098] The calculation formula of the theoretical solution for the deformation of thick-walled cylindrical structures is:

[0099]

[0100] where u r is the radial displacement of a point on the cylinder; v is the Poisson's ratio of the material; P is the internal pressure on the cylinder; R o is the outer radius of the cylinder; R i is the inner radius of the cylinder; E is the Young's modulus of the material; r is the distance from the point to the center of the cylinder.

[0101] According to the formula, the maximum theoretical value of the radial displacement of the axisymmetric thick-walled cylinder under internal pressure is 0.053340 mm, and its direction is away from the center of the cylinder.

[0102] When performing finite element calculations, considering that the thick-walled cylinder has a symmetrical structure, 1 / 4 of its structure is used for analysis. Utilizing the symmetric properties, the entire model is fixed in the Z direction (the height direction of the cylinder), the side parallel to the XOZ plane is fixed in the Y direction, and the side parallel to the YOZ plane is fixed in the X direction.

[0103] refer to Figure 2In Example 1, hexahedral unit division is first performed, and then each hexahedral mesh is divided into two triangular prism units, where the inner and outer walls are divided into triangles. There are three numbers of hexahedral unit divisions: 8×4×1 / 16×8×2 / 32×16×2 (number of units in the circumferential direction×radial direction×thickness direction). The calculation results are shown in Table 1 below:

[0104] Table 1: Comparison of results between ANSYS and this method

[0105]

[0106]

[0107] refer to Figure 3 In Example 2, hexahedral unit division is first performed, and then each hexahedral mesh is divided into six pyramid units (each face of the hexahedron is a quadrilateral face of a pyramid unit). There are three numbers of hexahedral unit divisions: 8×4×1 / 16×8×2 / 32×16×2 (number of units in the circumferential direction×radial direction×thickness direction). The calculation results are shown in Table 2 below:

[0108] Table 2: Comparison of results between ANSYS and this method

[0109]

[0110] According to the comparison of the results, it can be seen that for the compression problem of a nearly incompressible cylinder, the accuracy and convergence of the method disclosed in the present invention are significantly better than those of ANSYS, and are closer to the theoretical solution.

[0111] Figure 4 FIG1 shows an internal structure diagram of a computer device in an embodiment. The computer device can be a terminal or a server. Figure 4 As shown, the computer device includes a processor, a memory, and a network interface connected via a system bus. The memory includes a non-volatile storage medium and an internal memory. The non-volatile storage medium of the computer device stores an operating system and may also store a computer program. When the computer program is executed by the processor, the processor can implement the finite element calculation method of the present invention. The internal memory may also store a computer program. When the computer program is executed by the processor, the processor can implement the finite element calculation method of the present invention. It will be understood by those skilled in the art that Figure 4 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than shown in the figure, or combine certain components, or have a different component arrangement.

[0112] In one embodiment, the present invention provides a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the following steps:

[0113] S1: constructing a geometric model, assigning parameters of the nearly incompressible material to the geometric model, assigning boundary conditions and loads to the geometric model, and obtaining a finite element parameter model of the nearly incompressible material;

[0114] S2: Meshing the finite element parameter model according to the pre-selected unit type to obtain a finite element model;

[0115] S3: Establish the element stiffness matrix of the element:

[0116]

[0117] The element stiffness matrix k e is a function of the shape function derivative matrix B and the elasticity matrix D, and the integration domain is the unit, wherein the shape function derivative matrix B of each unit is a matrix composed of m 6×3 sub-matrices.

[0118] B=[B1 B2 … B m ]

[0119] Where m is the number of nodes of the unit type, and the submatrix B of each node is i (i=1~m) is:

[0120]

[0121] Among them B j (j=1~9) is the element shape function N i (i=1~m) for the global coordinate x s (s=1~3), where

[0122]

[0123] B5=B1+B4,

[0124] B7=B2+B6, B9=B3+B8

[0125] The integration domain Ω e is the unit, N i Determined by the unit type;

[0126] S4: Assemble the unit stiffness matrix to form a global stiffness matrix K; discretize the load to generate a global node force vector r, and calculate the node displacement vector u=K -1 r;

[0127] S5: performing post-processing using the node displacement vector to obtain the displacement, stress and strain of the finite element model.

[0128] The above-described embodiments merely represent several implementation methods of the present invention. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent application. It should be noted that a person skilled in the art would be able to make various modifications and improvements without departing from the spirit of the present invention, all of which fall within the scope of protection of the present invention. Therefore, the scope of protection of the patent for this invention shall be determined by the appended claims.

Claims

1. A finite element calculation method for nearly incompressible materials, characterized in that: The following steps are involved: S1: constructing a geometric model, assigning parameters of the nearly incompressible material to the geometric model, assigning boundary conditions and loads to the geometric model, and obtaining a finite element parameter model of the nearly incompressible material; S2: Meshing the finite element parameter model according to the pre-selected unit type to obtain a finite element model; S3: Establish the element stiffness matrix of the element: The element stiffness matrix k e is a function of the shape function derivative matrix B and the elasticity matrix D, and the integration domain is the unit, wherein the shape function derivative matrix B of each unit is a matrix composed of m 6×3 sub-matrices, B=[B1 B2 … B m ] Where m is the number of nodes of the unit type, and the submatrix B of each node is i (i=1~m) is: Among them B j (j=1~9) is the element shape function N i (i=1~m) for the global coordinate x s (s=1~3), where The integration domain Ω e is the unit, N i Determined by the unit type; S4: Assemble the element stiffness matrix to form a global stiffness matrix K; discretize the load to generate a global node force vector r, and calculate the node displacement vector u=K -1 r; S5: performing post-processing using the node displacement vector to obtain the displacement, stress and strain of the finite element model.

2. The finite element calculation method according to claim 1, characterized in that: In step S3, the element stiffness matrix k is calculated from the shape function derivative matrix B and the elastic matrix D. e In the process, numerical integration method is used for calculation.

3. The finite element calculation method according to claim 2, characterized in that: The numerical integration method is a Gaussian integration method, a Simpson integration method or a midpoint integration method.

4. The finite element calculation method according to claim 3, characterized in that: In step S3, when the Gaussian integral method is used for calculation, the Gaussian integral expression of the element stiffness matrix is: Among them, B is the shape function derivative matrix, D is the elasticity matrix, n is the number of Gaussian integration points, w i is the weight of the i-th Gaussian integration point, B i is the shape function derivative matrix at the i-th Gaussian integration point, B i T It's B i The transposed matrix, J i is the Jacobian matrix of the i-th Gaussian integration point, |J i |It's J i The determinant value of , ξ1, ξ2, ξ3 are the three coordinate axes in the parameter coordinate system of the unit, and x1, x2, x3 are the three coordinate axes of the global coordinate system of the unit.

5. The finite element calculation method according to claim 4, characterized in that: The element type is a first-order triangular prism element with a total of 6 nodes. The node shape function of the node is: Wherein ξ1, ξ2, ξ3 are the three coordinate axes in the parameter coordinate system of the unit.

6. The finite element calculation method according to claim 5, characterized in that: The first-order triangular prism element uses 6 Gaussian integration points with a weight of 1 / 6. In the parametric coordinate system, the coordinates of the Gaussian integration points are:

7. The finite element calculation method according to claim 4, characterized in that: The unit type is a first-order pyramid unit with a total of 5 nodes. The node shape function of the node is: Wherein ξ1, ξ2, ξ3 are the three coordinate axes in the parameter coordinate system of the unit.

8. The finite element calculation method according to claim 7, characterized in that: Eight Gaussian integration points are used, and the weights are all 1. In the parameter coordinate system, the coordinates of the Gaussian integration points are 9. The finite element calculation method according to claim 1, characterized in that: The parameters include shear modulus and Poisson's ratio; Alternatively, the parameters include shear modulus and Lame coefficient.

10. A computer device comprising a memory and a processor, characterized in that: The memory stores a computer program, and when the computer program is executed by the processor, the processor is caused to perform the steps of the method according to any one of claims 1 to 9.

Images