A spline finite element method for contact stress calculation

Abstract

The present invention discloses a spline finite element method for contact stress calculation, comprising the following steps: establishing a geometric model of multiple deformable body contacts, obtaining spline units, inputting material parameters, applying displacement and surface force boundary conditions and specifying a main contact surface and a slave contact surface; on the slave contact surface, using a control point projection method to obtain the projection points on the deformable body corresponding to the control points of the spline unit; calculating a transformation matrix based on the projection points to obtain an interpolation spline basis function with interpolation characteristics associated with the projection points on the slave contact surface; calculating a system stiffness matrix using a hybrid form of conventional spline basis functions and interpolation spline basis functions, and applying contact constraints to solve and obtain contact stress. Under the condition of ensuring the geometric accuracy of the model, the present invention can effectively solve and calculate contact stress compared with the existing point-to-surface method, and avoids the heavy amount of calculation of the existing face-to-face contact analysis model in the spline finite element method.

Description

Technical Field

[0001] The invention relates to the field of computer-aided engineering, and in particular to a spline finite element method for contact stress calculation. Background Art

[0002] Finite element numerical simulation based on computational contact mechanics is one of the effective means to evaluate the contact stress of mechanical structures for mechanical structure design and optimization.

[0003] However, the standard finite element method requires the use of triangular / quadrilateral-based piecewise grids to discretize the structure, which results in geometric discretization errors, making the discretized grid model insufficient to describe the original geometric contact surface of the structure, making it difficult to ensure the accuracy and reliability of the contact stress analysis results. For this reason, it is necessary to refine the grid to improve the convergence and accuracy of the calculation, but mesh refinement leads to an increase in the scale of calculation, which further limits the engineering application of the finite element method in contact analysis. Summary of the invention

[0004] The technical problem to be solved by the present invention is to provide a spline finite element method for contact stress calculation, which can effectively solve and calculate the contact stress under the condition of ensuring the geometric accuracy of the model, compared with the existing point-to-surface method, and avoid the heavy calculation amount of the existing face-to-face contact analysis model in the spline finite element method.

[0005] In order to solve the above technical problems, the present invention provides a spline finite element method for contact stress calculation, comprising the following steps:

[0006] S1: Establish a geometric model of multiple deformable body contacts, obtain spline elements, input material parameters, apply displacement and traction boundary conditions, and specify the master contact surface and slave contact surface;

[0007] S2: On the contact surface, the control point projection method is used to obtain the projection points of the control points of the spline unit corresponding to the deformable body;

[0008] S3: Calculate the transformation matrix based on the projection points to obtain an interpolation spline basis function with interpolation characteristics associated with the projection points on the contact surface;

[0009] S4: The system stiffness matrix is ​​calculated using a hybrid form of conventional spline basis functions and interpolation spline basis functions, and contact constraints are imposed to solve for the contact stress.

[0010] Furthermore, in step S1, the deformable body in the geometric model is represented by non-uniform rational B-spline, ie, NURBS.

[0011] Further, in step S2, for the specified slave contact surface, its boundary information in the calculation area is obtained, including node vectors and control points, wherein the projection point is located on the boundary;

[0012] Furthermore, for the I-th projection point on the boundary, its parameter coordinates are: ξI+i belongs to the elements of the node vector, p is the order of the basis function, and then the physical coordinates of the projection point can be obtained by combining the conventional spline basis function.

[0013] Furthermore, in step S3, the specific steps of calculating the transformation matrix based on the projection points and obtaining the interpolation spline basis function are as follows:

[0014] S31: Calculate the transformation matrix according to the parameter coordinates of the projection point. The calculation formula is:

[0015]

[0016] Among them: R I，J =R J (ξ′ I ), (I = 1, 2, ..., m and J = 1, 2, ..., m, m is the number of projection points on the contact surface) is the conventional spline basis function J at the projection point ξ′ I The value at

[0017] S32: Calculate the interpolation spline basis function at the projection point using the transformation matrix obtained in S31 is the inverse of the transposed transformation matrix, R A is a conventional spline basis function, ξ is an unknown quantity in the general sense of the parameter coordinates of any projection point in the computational domain.

[0018] Furthermore, in step S4, the specific steps of calculating the system stiffness matrix are:

[0019] S41: Traverse all projection points from the contact surface, find a point on the main contact surface with the shortest distance from each projection point to the main contact surface, and calculate the contact gap function value g n ,

[0020]

[0021] in, is the external normal vector of the nearest point on the main contact surface from the contact surface projection point, is the node sequence of the contact element, the superscript (2) belongs to the main contact surface, (x I ,y I ) is the physical coordinate of the projection point from the contact surface, A is the shape function matrix,

[0022] The elements marked with (2) above belong to the shape function at the point on the main contact surface that is shortest from the projection point;

[0023] S42: Penalty function method is used to impose constraints on the projection points with negative contact gap function values, and the contact virtual work is calculated;

[0024] S43: Calculate the linearized form of the contact virtual work and obtain the contact stiffness matrix;

[0025] S44: Calculate the tangent stiffness matrix of the deformable body, and superimpose it with the contact stiffness matrix to obtain the final system stiffness matrix;

[0026] Among them, the tangent stiffness matrix is ​​obtained using the strain matrix B:

[0027]

[0028] in, is the interpolation spline basis function, A=1,2,…,n el , n el is the number of control points of the unit, x′=x or y, and the rest are conventional spline basis functions;

[0029] S45: After obtaining the final stiffness matrix of the system from S44, a nonlinear set of equations for the contact problem is formed and solved iteratively using the Newton-Raphson method.

[0030] Furthermore, in S42, the contact virtual work is calculated as:

[0031]

[0032] where ∈ n is the penalty factor.

[0033] Furthermore, in S43, the linearized form of the contact virtual work is:

[0034]

[0035] in, is the metric tensor, k is the curvature of the principal contact surface at the point closest to the projection point.

[0036] Furthermore, the tangent stiffness matrix of the deformed body is calculated using the linear elastic small deformation or finite elastic deformation assumption.

[0037] Beneficial effects of the present invention:

[0038] (1) Different from the existing spline finite element method for point-to-surface contact, the present invention interpolates at the projection point of the contact surface, and the contact force at the projection point has direct physical meaning;

[0039] (2) Using spline functions that describe geometric data to approximate physical quantities in finite element analysis can avoid geometric errors caused by the discretization of triangular / quadrilateral meshes, improve the accuracy and reliability of calculation results, and thus promote the development of spline analysis technology.

[0040] (3) Compared with the existing spline finite element method for point-to-surface contact, the present invention can improve the accuracy of contact force and the convergence of the Newton-Raphson method for solving nonlinear contact problems, while eliminating the computational complexity of the existing spline finite element method for face-to-face contact;

[0041] (4) The geometric accuracy of the model is guaranteed, making it easy to implement adaptive analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a schematic flow chart of the spline finite element method for contact stress analysis of the present invention;

[0043] Figure 2 is a schematic diagram of a contact model of two blocks in the first embodiment of the present invention;

[0044] Figure 3 is a schematic diagram of the change of contact stress with penalty factor obtained by the conventional point-to-surface spline contact model compared with the first embodiment;

[0045] Figure 4 This is a schematic diagram of the change of the contact stress with the penalty factor obtained in Example 1 of the present invention;

[0046] Figure 5 This is a schematic diagram of a slider-substrate sliding contact model according to a second embodiment of the present invention;

[0047] Figure 6 It is a schematic diagram comparing the total contact pressure in the vertical direction of the slider-substrate sliding contact model in the sliding stage of Example 2. DETAILED DESCRIPTION

[0048] The present invention is further described below in conjunction with the accompanying drawings and specific embodiments so that those skilled in the art can better understand the present invention and implement it, but the embodiments are not intended to limit the present invention.

[0049] Embodiment 1

[0050] Reference Figure 1 As shown, an embodiment of the spline finite element method for contact stress calculation of the present invention includes the following steps:

[0051] First, a geometric model of multiple deformable bodies in contact is established, spline units are obtained, material parameters are input, displacement and surface force boundary conditions are applied, and the main contact surface and the secondary contact surface are specified. The object of this embodiment is two cubes with a side length of L = 1m in contact and pressure, and the top pressure of the upper cube is P = 10N / m 2 , the bottom of the block below and the left sides of the two blocks are all normal constraints. Non-uniform rational B-splines are used to describe the deformable body, and the spline unit distribution is as follows Figure 2 As shown, the upper block is the area where the contact surface is located, and its number of spline units is greater than that of the lower block. Material parameters: Young's modulus E = 1000N / m 2 With Poisson's ratio v = 0.3, plane strain condition.

[0052] On the slave contact surface, the control point projection method is used to calculate the projection point corresponding to the control point on the slave contact surface. Specifically, for the specified slave contact surface, the boundary information in the calculation area is obtained, including the node vector and the control point. For this example, the node vector of the slave contact surface is [0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1]. For the Ith projection point on the boundary, its parameter coordinates are: ξ I+i The elements of the node vector, p is the order of the basis function. For this example, the parameter coordinates of all projection points are [0, 0.125, 0.375, 0.625, 0.875, 1]; then combined with the conventional spline basis function, the physical coordinates of the projection points can be obtained for contact detection. The projection points are as follows: Figure 2 Circle shown.

[0053] Then the transformation matrix is ​​calculated according to the parameter coordinates of the projection point, and the calculation formula is:

[0054]

[0055] Among them: R I，J =R J (ξ′ I ), (I = 1, 2, ..., m and J = 1, 2, ..., m, m is the number of projection points on the contact surface) is the conventional spline basis function J at the projection point ξ′ I The value at the projection point is calculated by the obtained transformation matrix. Where m is the number of projection points from the contact surface, is the inverse of the transposed transformation matrix, R A is a conventional spline basis function, ξ is an unknown quantity in the general sense of the parameter coordinates of any projection point in the computational domain. That is, an interpolation spline basis function with interpolation characteristics associated with the projection point from the contact surface is obtained;

[0056] Finally, the system stiffness matrix is ​​calculated by using a hybrid form of conventional spline basis functions and interpolation spline basis functions and contact constraints are imposed to obtain the contact stress. For this embodiment, the system tangent stiffness matrix is ​​calculated using linear elastic small deformation;

[0057] The specific steps of calculating the stiffness matrix of the above system are as follows: first traverse all the projection points on the contact surface, find a point on the main contact surface with the shortest distance from each projection point to the main contact surface, and calculate the contact gap function value g n ,

[0058] in, is the external normal vector of the nearest point on the main contact surface from the contact surface projection point, is the node sequence of the contact element, the superscript (2) belongs to the main contact surface, (x I ,y I ) is the physical coordinate of the projection point from the contact surface, A is the shape function matrix,

[0059] The elements marked with (2) above belong to the shape function at the point on the main contact surface that is shortest from the projection point;

[0060] Then, the penalty function method is used to impose constraints on the projection points with negative contact gap function values, and the contact virtual work is calculated. The contact stiffness matrix is ​​obtained by using the linearized form of calculating the contact virtual work.

[0061] The contact virtual work is calculated as:

[0062]

[0063] where ∈ n is the penalty factor.

[0064] Then the tangent stiffness matrix of the deformed body is calculated and superimposed with the contact stiffness matrix to obtain the final system stiffness matrix;

[0065] The tangent stiffness matrix is ​​obtained using the strain matrix B:

[0066]

[0067] in, is the interpolation spline basis function, A=1,2,…,n el , x′=x or y, and the rest are conventional spline basis functions.

[0068] Finally, after obtaining the final stiffness matrix of the system, a nonlinear set of equations for the contact problem is formed and solved iteratively using the Newton-Raphson method;

[0069] The linearized form of contact virtual work is:

[0070]

[0071] in, is the metric tensor, κ is the curvature of the principal contact surface at the point closest to the projection point;

[0072] The contact stresses are calculated by postprocessing the solution of the contact nonlinear equations.

[0073] Reference Figure 3 and Figure 4 As shown in the figure, the conventional point-to-surface spline contact method and the method proposed in the present invention are respectively shown. For this embodiment, as the penalty factor ∈ n Take 10 respectively 3 , 10 4 with 10 5 When the contact stress distribution on the contact surface is , it can be seen that as the penalty factor increases, the contact stress of the conventional method cannot converge to an exact solution, while the method proposed in the present invention can quickly obtain accurate contact stress as the penalty factor increases.

[0074] Embodiment 2:

[0075] An embodiment of a spline finite element method for contact stress analysis is also provided, comprising the following steps:

[0076] First, establish a geometric model of multiple deformable body contacts, obtain spline units, input material parameters, apply displacement and surface force boundary conditions, and specify the main contact surface and the secondary contact surface. The object of this embodiment is radius R = 1m, Young's modulus E = 50N / m 2 The module with a semicircular cross section is compressed and slides on a base with a size of 6m×2m, such as Figure 5 As shown. Matrix Young's modulus E = 1N / m 2 The bottom of the base is a fixed constraint, and the two sides are normal constraints. The surface of the base that contacts the module is the main contact surface, and the top of the module is first subjected to downward displacement. Then keep the horizontal displacement to the right. The lower arc surface is the secondary contact surface. NURBS splines are used to describe the deformable body, and plane strain assumptions are used.

[0077] For the contact surface, the control point projection method is used to calculate the projection point corresponding to the control point on the contact surface, which is consistent with the above embodiment and will not be described in detail;

[0078] Then, a transformation matrix is ​​calculated based on the projection points, thereby obtaining an interpolation spline basis function with interpolation characteristics associated with the projection points from the contact surface, which is consistent with the above embodiment and will not be described in detail;

[0079] Then, the system stiffness matrix is ​​calculated by using a hybrid form of conventional spline basis functions and interpolation spline basis functions, and contact constraints are imposed, and the contact stress is obtained by solving. In this embodiment, the matrix tangent stiffness matrix is ​​calculated by finite elastic deformation assumption, the material is a compressible Neo-Hooke constitutive relation, and the penalty factor ∈ n Take 10 2 .

[0080] Table 1 shows the number of incremental steps required for solving three different grid models using the Newton-Raphson method in Example 2:

[0081]

[0082]

[0083] Table 1 shows the number of incremental steps required for solving three mesh models with different numbers of units using the Newton-Raphson method. It can be seen that the number of incremental steps required by the method of the present invention is less than that of the conventional point-to-surface method, whether in the vertical downward pressure stage or the horizontal sliding stage, indicating that the method of the present invention can promote the convergence of contact nonlinear calculations. Figure 6 The comparison of the contact stress obtained by the conventional point-to-surface contact method and the method proposed in the present invention during the sliding process is shown, and it can be seen that the results obtained by the method of the present invention are consistent with the conventional method. Based on the above conclusions, it can be seen that compared with the prior art, the present invention interpolates from the projection point of the contact surface, and the contact force at the projection point has a direct physical meaning; the method of the present invention can improve the accuracy of the contact force and improve the convergence of the Newton-Raphson method for solving nonlinear contact problems. At the same time, it does not require the computational complexity of the existing spline finite element method for face-to-face contact, which ensures the geometric accuracy of the model and makes it easy to implement adaptive analysis.

[0084] The above embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Any equivalent substitution or change made by a person skilled in the art based on the present invention is within the protection scope of the present invention.

Claims

1. A spline finite element method for contact stress calculation, characterized in that: The steps include: S1: Establish a geometric model of multiple deformable body contacts, obtain spline elements, input material parameters, apply displacement and traction boundary conditions, and specify the master contact surface and slave contact surface; S2: On the contact surface, the control point projection method is used to obtain the projection points of the control points of the spline unit corresponding to the deformable body; S3: Calculate a transformation matrix based on the projection points to obtain an interpolation spline basis function with interpolation characteristics associated with the projection points on the contact surface; S4: The system stiffness matrix is ​​calculated using a hybrid form of conventional spline basis functions and interpolation spline basis functions, and contact constraints are imposed to obtain the contact stress; In step S2, for the specified slave contact surface, the boundary information in the calculation area is obtained, including the node vector and the control point, wherein the projection point is located on the boundary; For the I-th projection point on the boundary, its parameter coordinates are: ξ I+i The elements belonging to the node vector, p is the order of the basis function, and then the physical coordinates of the projection point are obtained by combining the conventional spline basis function; In step S3, the specific steps of calculating the transformation matrix based on the projection points and obtaining the interpolation spline basis function are as follows: S31: Calculate the transformation matrix according to the parameter coordinates of the projection point. The calculation formula is: in; R I，J =R J (ξ′ I ), (I = 1, 2, ..., m and J = 1, 2, ..., m, m is the number of projection points on the contact surface) is the conventional spline basis function J at the projection point ξ′ I The value at S32: Calculate the interpolation spline basis function at the projection point using the transformation matrix obtained in S31 is the inverse of the transposed transformation matrix, R A is the conventional spline basis function, ξ is the parameter coordinate of any projection point in the computational domain; In step S4, the specific steps of calculating the system stiffness matrix are: S41: Traverse all projection points from the contact surface, find a point on the main contact surface with the shortest distance from each projection point to the main contact surface, and calculate the contact gap function value g n , in, is the external normal vector of the nearest point on the main contact surface from the contact surface projection point, is the node sequence of the contact element, the superscript (2) belongs to the main contact surface, (x I ,y I ) is the physical coordinate of the projection point from the contact surface, A is the shape function matrix, The elements marked with (2) above belong to the shape function at the point on the main contact surface that is shortest from the projection point; S42: Penalty function method is used to impose constraints on the projection points with negative contact gap function values, and the contact virtual work is calculated; S43: Calculate the linearized form of the contact virtual work and obtain the contact stiffness matrix; S44: Calculate the tangent stiffness matrix of the deformable body, and superimpose it with the contact stiffness matrix to obtain the final system stiffness matrix; Among them, the tangent stiffness matrix is ​​obtained using the strain matrix B: in, is the interpolation spline basis function, A=1,2,…,n el , n el is the number of control points of the unit, x′=x or y, and the rest are conventional spline basis functions; S45: After obtaining the final stiffness matrix of the system from S44, a nonlinear set of equations for the contact problem is formed and solved iteratively using the Newton-Raphson method.

2. The spline finite element method for contact stress calculation according to claim 1, characterized in that: In step S1, the deformable body in the geometric model is represented by non-uniform rational B-splines, ie, NURBS.

3. The spline finite element method for contact stress calculation according to claim 1, characterized in that: In S42, the contact virtual work is calculated as: where ∈ n is the penalty factor.

4. The spline finite element method for contact stress calculation according to claim 3, characterized in that: In S43, the linearized form of contact virtual work is: in, is the metric tensor, κ is the curvature of the principal contact surface at the point closest to the projection point.

5. The spline finite element method for contact stress calculation according to claim 1, characterized in that: The tangent stiffness matrix of the deformed body is calculated using the linear elastic small deformation or finite elastic deformation assumption.

Images