A method for solving mechanical problems based on hybrid finite element method

Abstract

The present application discloses a method for solving mechanical problems based on a hybrid finite element method. It includes: analyzing each unit of a force-bearing body, constructing a coefficient matrix and a load vector; processing the boundary conditions on the boundary units, adjusting the load vector and obtaining a set of constraint equations; analyzing each unit interface of a Lagrangian unit and an intrinsic hybrid finite element unit, obtaining a local coefficient matrix of each unit interface, and assembling the local coefficient matrix of each unit interface into a coefficient matrix; processing the load vector and the coefficient matrix based on the set of constraint equations to obtain a set of linear algebraic equations and solving them to obtain stress and displacement results. The present application combines the Lagrangian unit method and the intrinsic hybrid finite element method to solve the linear elastic mechanics problem of a force-bearing body, and solves the problem of low accuracy of the results obtained by using the traditional finite element method to solve the linear elastic mechanics problem in the related technology while ensuring the solution efficiency.

Description

Technical Field

[0001] The present application relates to the technical field of structural mechanics, and in particular to a method for solving mechanical problems based on a hybrid finite element method. Background Art

[0002] Linear elasticity mainly studies the stress, strain and displacement of linear elastic bodies under the action of specific external factors, so as to solve the problems of stiffness, strength and stability in the design of products or engineering structures. It has always been the focus of research for technicians in the fields of vehicles, civil engineering, aerospace, etc. However, the classical solutions to linear elasticity problems have limitations. When solving complex geometric domains and complex boundaries, analytical solutions are often not obtained, so various numerical solutions have emerged.

[0003] As one of the most widely used numerical solution methods, the finite element method uses the variational principle as the basic idea for solving linear elasticity problems, that is, converting the problem of solving partial differential equations into the problem of solving functional energy minimization. The solution is different from the previous method of finding approximate functions over the entire area. The finite element method divides the area into a finite number of units, analyzes each unit, and then assembles them.

[0004] At present, the relevant technologies generally use intrinsic mixed finite elements or linear elements to solve linear elastic mechanics problems. Among them, although the intrinsic mixed finite element has a high solution accuracy, it has a large computational cost and low efficiency, resulting in a small practical application range; while the linear element has a low computational cost, but due to its low order, the error of the solution result is large. Therefore, there is an urgent need for a method that can take into account both high efficiency and high accuracy in solving linear elastic mechanics problems.

[0005] To address the above-mentioned problems, no effective solution has been proposed yet. Summary of the invention

[0006] The embodiment of the present application provides a method for solving mechanical problems based on a hybrid finite element method, so as to at least solve the technical problem that the results obtained by the related art using the traditional finite element method to solve linear elastic mechanics problems are of low accuracy.

[0007] According to one aspect of an embodiment of the present application, a method for solving mechanical problems based on a hybrid finite element method is provided, comprising: step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide a force-bearing body, wherein the unit data at least includes: unit shape, node coordinates, number of units, unit type, unit dimension, and the unit type includes: intrinsic hybrid finite element type, Lagrangian element type, material parameters include: Young's modulus, Poisson's ratio, density, and boundary conditions include: surface force, displacement; step S2, looping each unit, obtaining multiple integration points on the current unit, and combining the unit data, material parameters and body forces obtained in step S1 to determine the value of the body force at each integration point and the value corresponding to the unit type The values, gradient values ​​or divergence values ​​of the finite element basis functions at each integration point include: the values ​​of the intrinsic mixed finite element basis functions at each integration point, the divergence values ​​of the intrinsic mixed finite element basis functions at each integration point, the values ​​of the discontinuous Lagrangian basis functions at each integration point, the values ​​of the continuous Lagrangian basis functions at each integration point and the gradient values ​​of the continuous Lagrangian basis functions at each integration point; the overall degrees of freedom of the finite element basis functions of each unit are determined according to the unit type, and the coefficient matrix A and the load vector b are constructed in combination with the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis functions corresponding to the unit type at each integration point in the unit; step S4, for each boundary unit Step S5: loop through each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit to obtain the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface; Step S6: calculate the values ​​of the k1-order Lagrangian unit basis function and the k2-order intrinsic hybrid finite element basis function on the current unit interface; Step S7: combine the Lagrangian unit and the intrinsic hybrid finite element unit to obtain the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface. Steps S5 to S6 determine the local coefficient matrix C of the current unit interface; Step S8, based on the overall degrees of freedom of the finite element basis function of the unit obtained in step S3 and the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the current unit interface determined in step S6; Step S9, based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8, assemble the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A; Step S10, using the constraint equation group Gx=s obtained in step S4, process the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain the linear algebraic equation group Ax=b;Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results. ;

[0008] Optionally, the above step S2 includes the following steps: traverse each unit and perform the following steps: step S21, for each current unit, obtain the unit volume, unit order k, multiple integration points and the integration weight of each integration point of the current unit; step S22, using the physical force obtained in step S1, determine the value of the physical force at each integration point on the current unit; step S23, determine the unit type of the current unit; step S24, when the unit type is a Lagrangian element type, use step S21 to calculate the value of the k-order Lagrangian element basis function of the current unit at each integration point and the k-order Lagrangian element basis function The gradient value of the function at each integration point; Step S25, when the unit type is an intrinsic hybrid finite element type, obtain the tensor basis of the three-dimensional symmetric matrix space of the current unit, and combine the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point calculated in step S21 to obtain the value of the intrinsic hybrid finite element basis function of the current unit at each integration point and the divergence value of the intrinsic hybrid finite element basis function at each integration point, and at the same time use step S21 to calculate the value of the k-1-order discontinuous Lagrangian basis function of the current unit at each integration point.

[0009] Optionally, the above-mentioned step S25 includes the following steps: step S251, using step S21 to calculate the value of the k-order Lagrangian basis function at each integration point, the gradient value of the k-order Lagrangian basis function at each integration point, and the value of the k-1-order discontinuous Lagrangian basis function at each integration point; step S252, determining the type of each interpolation point on the current unit, wherein the type of interpolation point includes: unit vertex, unit internal point, point on the unit edge and point on the unit surface; step S253, determining the tensor basis of the three-dimensional symmetric matrix space according to the type of interpolation point; step S254, combining the tensor basis of the three-dimensional symmetric matrix space corresponding to each type of interpolation point with the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point obtained in step S251, to obtain the value of the intrinsic mixed finite element basis function at each integration point and the divergence value of the intrinsic mixed finite element basis function at each integration point.

[0010] Optionally, the above step S253 includes the following steps: Step 1: When the type of the interpolation point is a unit vertex and a unit internal point, obtain the standard basis of the third-order symmetric tensor space; Step 2: When the type of the interpolation point is a point on the unit edge, obtain a unit tangent vector and two unit normal vectors for each of the six edges of the current unit, and use the unit tangent vector and two unit normal vectors of each edge to determine the edge basis of the unit edge; Step 3: When the type of the interpolation point is a point on the unit face, obtain two unit tangent vectors and one unit normal vector for each of the four faces of the current unit, and use the two unit tangent vectors and one unit normal vector of each face to determine the face basis of the unit face.

[0011] Optionally, the above step S3 includes the following steps: traverse each unit and execute the following steps: step S31, determine the unit type of the current unit; step S32, when the unit type is a Lagrangian element type, construct the coefficient matrix A and the load vector b according to the following steps, including: step S321, sort the interpolation points corresponding to the k-order Lagrangian element basis function according to the order of continuous degrees of freedom to obtain the overall degrees of freedom of the k-order Lagrangian element basis function of the current unit; step S322, use the gradient value of the k-order Lagrangian element basis function obtained in step S24 at each integration point to calculate the unit stiffness matrix , and combined with the overall degrees of freedom of the k-order Lagrangian basis function obtained in step S321, the calculated unit stiffness matrix is ​​assembled into the coefficient matrix A; step S323, using the k-order Lagrangian basis function obtained in step S24 and the value of the body force at each integration point to calculate the unit load vector, and combined with the overall degrees of freedom of the k-order Lagrangian basis function obtained in step S321, the calculated unit load vector is assembled into the load vector b; step S33, when the unit type is an intrinsic mixed finite element type, the coefficient matrix A and the load vector b are constructed according to the following steps, including: step S331, determine when The overall degrees of freedom of the intrinsic mixed finite element basis function of the previous unit; step S332, using the values ​​of the intrinsic mixed finite element basis function obtained in step S25 at each integration point to calculate the unit coefficient matrix, and combining step S331 to construct the coefficient matrix A; step S333, sorting the interpolation points corresponding to the k-1 order discontinuous Lagrangian basis function obtained in step S25 according to the order of discontinuous degrees of freedom, to obtain the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function of the current unit; step S334, using the values ​​of the k-1 order discontinuous Lagrangian basis function obtained in step S25 at each integration point, the intrinsic mixed The unit divergence operator matrix is ​​calculated by combining the divergence values ​​of the finite element basis function at each integration point, and the unit divergence operator matrix is ​​assembled into the coefficient matrix A in combination with the overall degrees of freedom of the intrinsic mixed finite element basis function obtained in step S331 and the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function obtained in step S333; step S335, the unit load vector is calculated using the k-1 order discontinuous Lagrangian basis function obtained in step S25 and the value of the body force at each integration point, and the unit load vector is assembled into the load vector b in combination with the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function obtained in step S333.

[0012] Optionally, the above step S331 includes the following steps: Step 1: determine the type of interpolation point corresponding to the intrinsic hybrid finite element basis function; Step 2: when the type of the interpolation point is a unit vertex and a unit internal point, determine the overall degree of freedom of the intrinsic hybrid finite element basis function according to the order of continuous degree of freedom; Step 3: when the type of the interpolation point is a point on the unit edge or a point on the unit surface, determine whether the intrinsic hybrid finite element basis function is composed of pure tangential vectors; Step 4: when the intrinsic hybrid finite element basis function is composed of pure tangential vectors, determine the overall degree of freedom of the intrinsic hybrid finite element basis function according to the order of discontinuous degree of freedom; Step 5: when the intrinsic hybrid finite element basis function is not composed of pure tangential vectors, determine the overall degree of freedom of the intrinsic hybrid finite element basis function according to the order of continuous degree of freedom.

[0013] Optionally, the above step S4 includes the following steps: traverse each boundary condition of each boundary unit and execute the following steps: step S41, determine the unit type of the boundary unit to which the current boundary condition belongs; step S42, when the unit type is an intrinsic hybrid finite element type, process the current boundary condition according to the following steps, including: step S421, when the current boundary condition is displacement, obtain the unit external normal vector of the current boundary, call the intrinsic hybrid finite element basis function, and load the integral result of the intrinsic hybrid finite element basis function, the unit external normal vector and the displacement into the load vector b; step S422, when the current boundary condition is displacement, When it is a surface force, the least squares method is called to process the surface force, and the processing result is added to the constraint equation group Gx=s; step S43, when the unit type is a Lagrangian element type, the current boundary condition is processed according to the following steps, including: step S431, when the current boundary condition is a surface force, the k-order Lagrangian element basis function of the Lagrangian element is called, and the integral result of the k-order Lagrangian element basis function and the surface force is loaded into the load vector b; step S432, when the current boundary condition is a displacement, the least squares method is called to process the displacement, and the processing result is added to the constraint equation group Gx=s.

[0014] Optionally, the above step S7 includes the following steps: traverse each unit interface and execute the following steps: step S71, calculate the product of the unit normal vector obtained in step S5 and the value of the k2-order intrinsic mixed finite element basis function obtained in step S6 at each integration point in the current unit interface; step S72, calculate the product result and the value of the k1-order Lagrangian basis function obtained in step S6 at each integration point in the current unit interface, and use the obtained product result as the local coefficient matrix C of the current interface.

[0015] According to another aspect of an embodiment of the present application, a computer program product is also provided, which includes a stored computer program, wherein when the computer program is executed by a processor, the method for solving mechanical problems based on the hybrid finite element method is implemented.

[0016] According to another aspect of an embodiment of the present application, an electronic device is also provided, which includes: a memory and a processor, wherein a computer program is stored in the memory, and the processor is configured to execute the above-mentioned method for solving mechanical problems based on the hybrid finite element method through the computer program.

[0017] Compared with the prior art, the present invention exchanges boundary information on the unit interface between the Lagrangian element type unit and the intrinsic hybrid finite element type unit, so that the entire algorithm combines the flexibility of the Lagrangian element and the high precision of the intrinsic hybrid finite element, thereby solving the technical problem of low precision of the results obtained by the related technology using the traditional finite element method to solve linear elastic mechanics problems while ensuring the solution efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] The drawings described herein are used to provide a further understanding of the present application and constitute a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0019] Figure 1 It is a hardware structure block diagram of a computer terminal (or mobile device) for implementing a method for solving mechanical problems based on a hybrid finite element method according to an embodiment of the present application;

[0020] Figure 2 It is a flow chart of an optional method for solving mechanical problems based on the hybrid finite element method according to an embodiment of the present application;

[0021] Figure 3 is a flow chart of an optional step S2 according to an embodiment of the present application;

[0022] Figure 4 is a flowchart of an optional step S25 according to an embodiment of the present application;

[0023] Figure 5 is a flowchart of an optional step S253 according to an embodiment of the present application;

[0024] Figure 6 is a schematic diagram of an optional tetrahedral unit according to an embodiment of the present application;

[0025] Figure 7 is a schematic diagram of an optional tangent normal vector of a tetrahedron unit according to an embodiment of the present application;

[0026] Figure 8 is a flowchart of an optional step S3 according to an embodiment of the present application;

[0027] Fig. 9 is a flowchart of an optional step S331 according to an embodiment of the present application;

[0028] Fig.10 is a flow chart of an optional step S4 according to an embodiment of the present application;

[0029] Fig.11a It is a schematic diagram of the structure of an L-shaped structure;

[0030] Fig.11b It is a schematic diagram of a finite element model of an L-shaped structure;

[0031] Fig.12 It is a schematic diagram of mesh division of an L-type structure model No. 1 analyzed by Lagrange element method;

[0032] Fig.13 It is a schematic diagram of the mesh distribution of the L-type structure model No. 2 analyzed by Lagrange element method;

[0033] Fig.14 It is a schematic diagram of mesh division of the L-type structure model No. 3 using Lagrange element method;

[0034] Fig.15 It is a schematic diagram of mesh division of the L-type structure model No. 4 using Lagrange element method;

[0035] Fig.16 It is a schematic diagram of mesh division of the L-type structure model No. 5 analyzed by Lagrange element method;

[0036] Fig.17 It is a schematic diagram of mesh division of the L-type structure model No. 6 using Lagrange element method;

[0037] Fig.18 It is a schematic diagram of mesh division of the L-type structure model No. 7 using Lagrange element method;

[0038] Fig.19 It is a schematic diagram of stress distribution analysis of an L-type structural model using a hybrid finite element method. DETAILED DESCRIPTION

[0039] In order to enable those skilled in the art to better understand the solution of the present application, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of the present application.

[0040] It should be noted that the terms "first", "second", etc. in the specification and claims of the present application and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that the data used in this way can be interchangeable where appropriate, so that the embodiments of the present application described herein can be implemented in an order other than those illustrated or described herein. In addition, the terms "including" and "having" and any of their variations are intended to cover non-exclusive inclusions, for example, a process, method, system, product or device comprising a series of steps or units is not necessarily limited to those steps or units clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0041] In addition, the relevant information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data for display, data for analysis, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties. For example, an interface is set up between this system and the relevant users or organizations. Before obtaining relevant information, it is necessary to send an acquisition request to the aforementioned user or organization through the interface, and obtain relevant information after receiving the consent information fed back by the aforementioned user or organization.

[0042] In order to better understand the embodiments of the present application, the technical terms involved in the embodiments of the present application are explained as follows:

[0043] Interpolation points: a concept used in interpolation theory, which refers to a set of known data points that are used to estimate or construct a function through interpolation methods so that the function takes known values ​​at these points. In simple terms, interpolation points are a set of data points used to generate an approximate function.

[0044] Example 1

[0045] According to an embodiment of the present application, an embodiment of a method for solving mechanical problems based on a hybrid finite element method is provided. It should be noted that the steps shown in the flowchart of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and although a logical order is shown in the flowchart, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0046] The method embodiments provided in the embodiments of the present application can be executed in a mobile terminal, a computer terminal or a similar computing device. Figure 1 The hardware structure block diagram of a computer terminal (or mobile device) for implementing a method for solving mechanical problems based on a hybrid finite element method is shown. Figure 1 As shown, the computer terminal 10 (or mobile device 10) may include one or more (102a, 102b, ..., 102n are used to illustrate) processors 102 (the processor 102 may include but is not limited to a processing device such as a microprocessor MCU or a programmable logic device FPGA), a memory 104 for storing data, and a transmission device 106 for communication functions. In addition, it may also include: a display, an input / output interface (I / O interface), a universal serial bus (USB) port (which may be included as one of the ports of the BUS bus), a network interface, a power supply and / or a camera. It can be understood by those skilled in the art that Figure 1 The structure shown is only for illustration and does not limit the structure of the above electronic device. Figure 1 More or fewer components as shown, or with Figure 1 Different configurations shown.

[0047] It should be noted that the one or more processors 102 and / or other data processing circuits described above may generally be referred to herein as "data processing circuits". The data processing circuits may be embodied in whole or in part as software, hardware, firmware, or any other combination thereof. In addition, the data processing circuit may be a single independent processing module, or may be incorporated in whole or in part into any of the other components in the computer terminal 10 (or mobile device). As involved in the embodiments of the present application, the data processing circuit acts as a processor control (e.g., selection of a variable resistor terminal path connected to an interface).

[0048] The memory 104 can be used to store software programs and modules of application software, such as the program instructions / data storage device corresponding to the method for solving mechanical problems based on the hybrid finite element method in the embodiment of the present application. The processor 102 executes various functional applications and data processing by running the software programs and modules stored in the memory 104, that is, the method for solving mechanical problems based on the hybrid finite element method of the above-mentioned application program is realized. The memory 104 may include a high-speed random access memory, and may also include a non-volatile memory, such as one or more magnetic storage devices, flash memory, or other non-volatile solid-state memory. In some examples, the memory 104 may further include a memory remotely arranged relative to the processor 102, and these remote memories may be connected to the computer terminal 10 via a network. Examples of the above-mentioned network include, but are not limited to, the Internet, an intranet, a local area network, a mobile communication network, and a combination thereof.

[0049] The transmission device 106 is used to receive or send data via a network. The specific example of the above network may include a wireless network provided by a communication provider of the computer terminal 10. In one example, the transmission device 106 includes a network adapter (Network Interface Controller, NIC), which can be connected to other network devices through a base station so as to communicate with the Internet. In one example, the transmission device 106 can be a radio frequency (RF) module, which is used to communicate with the Internet wirelessly.

[0050] The display may be, for example, a touch screen liquid crystal display (LCD) that enables a user to interact with a user interface of the computer terminal 10 (or mobile device).

[0051] In the above operating environment, Figure 2 is a flow chart of an optional method for solving mechanical problems based on the hybrid finite element method according to an embodiment of the present application, such as Figure 2 As shown, the method at least includes steps S1-S11, wherein:

[0052] Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body.

[0053] In the technical solution provided in the above step S1, the above-mentioned force-bearing body refers to a force-bearing body with a curved surface. The above-mentioned unit data at least includes: unit shape, node coordinates, unit number, unit type, and the unit type includes: intrinsic mixed finite element type and Lagrangian finite element type. The above-mentioned material parameters include: Young's modulus, Poisson's ratio, density. The above-mentioned boundary conditions include: surface force and displacement.

[0054] Specifically, when simulating and analyzing structural mechanics, the type of mesh, mesh density, and mesh shape used to divide the force-bearing body play a vital role in the accuracy of the simulation results. In order to simplify the structure into a simulation result that meets the accuracy requirements and can be solved quickly, a mesh of a suitable shape can be selected according to the actual situation to divide the force-bearing area of ​​the force-bearing body. In the embodiment of the present application, a tetrahedral mesh (i.e., a unit shape) is preferably used to divide the force-bearing area of ​​the force-bearing body, and the number of units and unit dimensions obtained by the specific division can be determined according to the shape of the divided mesh and the size and dimension of the actual force-bearing area. In addition, the boundary conditions, material coefficients, and body forces of the units can be determined according to the constraints, material coefficients, and body forces of the actual force-bearing body in the force-bearing area.

[0055] Step S2, looping each unit, obtaining multiple integration points on the current unit, and combining the unit data, material parameters and body force obtained in step S1 to determine the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point.

[0056] In the technical solution provided in the above step S2, the above-mentioned determined values ​​include: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point.

[0057] Optionally, Figure 3 is a flow chart of an optional step S2 according to an embodiment of the present application, such as Figure 3 Traverse each unit of the force-bearing body and execute the following steps S21-S25, including:

[0058] Step S21, obtaining the unit order k, multiple integration points and the integration weight of each integration point of the current unit;

[0059] Step S22, using the physical strength obtained in step S1 and step S21 to determine the value of the physical strength at each integral point on the current unit;

[0060] Step S23, determining the unit type of the current unit, wherein when the unit type is a Lagrangian element type, executing step S24; when the unit type is an intrinsic mixed finite element type, executing step S25;

[0061] Step S24, using step S21, calculate the value of the k-order Lagrangian basis function of the current unit at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point.

[0062] Step S25 obtains the tensor basis of the three-dimensional symmetric matrix space of the current unit, and combines the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point calculated by step S21 to obtain the value of the intrinsic mixed finite element basis function at each integration point of the current unit and the divergence value of the intrinsic mixed finite element basis function at each integration point. At the same time, step S21 is used to calculate the value of the k-1-order discontinuous Lagrangian basis function of the current unit at each integration point.

[0063] It should be noted that there is no specific restriction on the execution order of the above steps S24 and S25.

[0064] Further, Figure 4 is a flow chart of an optional step S25 according to an embodiment of the present application, such as Figure 4 The above step S25 can be implemented according to the following steps S251-S254, including:

[0065] Step S251, using step S21 to calculate the value of the k-th order Lagrangian basis function at each integration point, the gradient value of the k-th order Lagrangian basis function at each integration point, and the value of the k-1-th order discontinuous Lagrangian basis function at each integration point.

[0066] Step S252, determining the type of each interpolation point on the current unit.

[0067] The types of interpolation points include: unit vertices, unit internal points, points on unit edges, and points on unit faces.

[0068] Step S253, determining the tensor basis of the three-dimensional symmetric matrix space according to the type of the interpolation point.

[0069] In the technical solution provided in the above step S253, the above third-order symmetric matrix space can be understood as a vector space composed of all third-order symmetric matrices. Therefore, the basis of the third-order symmetric matrix space refers to a set composed of several linearly independent third-order symmetric matrices, and these third-order symmetric matrices can represent all third-order symmetric matrices in the space, and any third-order symmetric matrix can be expressed as six specific symmetric matrices.

[0070] Furthermore, Figure 5 is a flow chart of an optional step S253 according to an embodiment of the present application, such as Figure 5 This step S253 can be subdivided into the following steps, including:

[0071] Step 1: Loop through the Lagrange interpolation points on the current cell.

[0072] Step 2: Determine the type of Lagrange interpolation point.

[0073] Among them, when the type of Lagrangian interpolation point is a unit vertex and a unit internal point, the third step is executed; when the type of Lagrangian interpolation point is a point on the unit edge, the fourth step is executed; when the type of Lagrangian interpolation point is a point on the unit surface, the fifth step is executed.

[0074] Step 3: Obtain the standard basis of the third-order symmetric tensor space.

[0075] Specifically, the standard basis is the tensor basis of the Lagrange interpolation point. For example, the standard basis can be a linear combination of the following matrices: , , , , , .

[0076] Step 4: Get a unit tangent vector and two unit normal vectors for each of the six edges of the current unit, and use the unit tangent vector and two unit normal vectors of each edge to determine the edge basis of the unit edge.

[0077] Specifically, Figure 6 The element edges of the tetrahedral element are shown For example, for this unit edge The process of determining the edge basis is described as follows.

[0078] First, define the unit edge A unit tangent vector on the current tetrahedral unit, and record the unit edge shared by the current tetrahedral unit The unit external normal vectors of the two faces of the unit edge The two unit normal vectors of Figure 7 As shown, their expressions are:

[0079]

[0080]

[0081]

[0082] Among them, the above Represents the unit edge The unit normal vector of Represents the unit edge The unit tangent vector of .

[0083] Then, through the unit tangent vector We can define a symmetric matrix of the first kind with a rank of 1, denoted as , which can be expressed as:

[0084]

[0085] At the same time, two unit normal vectors A class of symmetric matrices formed by Zhang and composed of unit normal vectors and the unit tangent vector A type of symmetric matrix that is jointly spanned is called a second-type symmetric matrix, and belongs to the third-order symmetric matrix space. The second-type symmetric matrix can be recorded as , which can be expressed as:

[0086]

[0087]

[0088]

[0089]

[0090]

[0091] Therefore, we can use the first kind of symmetric matrix and the symmetric matrix of the second kind Together they form a unit edge The edge base.

[0092] Step 5: Get two unit tangent vectors and one unit normal vector for each of the four faces of the current unit, and use the two unit tangent vectors and one unit normal vector of each face to determine the face basis of the unit face.

[0093] Specifically, Figure 6 The element faces of the tetrahedral element are shown For example, for this unit surface The process of determining the surface base is described as follows.

[0094] First, define the element surface The unit external normal vector on the unit surface is taken as a unit tangent vector, and the unit surface The first unit tangent vector of the first edge and the second unit tangent vector of the second edge are used as the unit face Two unit tangent vectors, such as Figure 7 As shown, their expressions are:

[0095] Among them, the above represents the unit normal vector, and Respectively represent unit tangent vectors. Among them, if

[0096] If not 0, .

[0097] Then, through the two unit tangent vectors We can define three rank-1 symmetric matrices of the third kind and denote them as , which can be expressed as:

[0098]

[0099]

[0100]

[0101] At the same time, the unit normal vector A class of symmetric matrices formed by Zhang and composed of unit normal vectors and two unit tangent vectors The symmetric matrix of the fourth kind is a symmetric matrix of the fourth kind, and the fourth kind of symmetric matrix belongs to the third-order symmetric matrix space. The fourth kind of symmetric matrix can be recorded as , whose expressions can be written as:

[0102]

[0103]

[0104]

[0105] Therefore, the third kind of symmetric matrix can be and the fourth kind of symmetric matrix Together they form a unit surface Surface base.

[0106] Step 6: Determine whether to loop to the last Lagrange interpolation point of the current unit. If yes, execute step 8, otherwise continue to execute step 2:

[0107] Step 8: Output the tensor basis of the third-order symmetric matrix space.

[0108] Step S254, combines the tensor basis of the three-dimensional symmetric matrix space corresponding to each type of interpolation point with the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point obtained in step S251, to obtain the value of the intrinsic mixed finite element basis function at each integration point and the divergence (div) value of the intrinsic mixed finite element basis function at each integration point.

[0109] In the technical solution provided in the above step S254, the value of the above intrinsic mixed finite element basis function at the integration point can be divided into the following cases according to the type of Lagrange interpolation point:

[0110] When the type of the Lagrangian interpolation point is a unit vertex, the value of the intrinsic mixed finite element basis function at the integration point is obtained according to the product of the standard basis determined in step S421 and the value of the k-order Lagrangian basis function corresponding to the current vertex at the integration point;

[0111] When the type of the Lagrangian interpolation point is an internal point of the unit, the value of the intrinsic mixed finite element basis function at the integration point is obtained according to the product of the standard basis determined in step S421 and the value of the k-order Lagrangian basis function corresponding to the current internal point of the unit at the integration point;

[0112] When the type of the Lagrangian interpolation point is a unit edge point, the value of the intrinsic mixed finite element basis function at the integration point is obtained by multiplying the edge basis of the current unit determined in step S422 by the value of the k-order Lagrangian basis function corresponding to the current unit edge point at the integration point;

[0113] When the type of the Lagrangian interpolation point is a point on the element surface, the value of the intrinsic mixed finite element basis function at the integration point is obtained by multiplying the surface basis of the current element determined in step S423 and the value of the k-order Lagrangian basis function corresponding to the point on the current element surface at the integration point.

[0114] In other words, if The three-dimensional force-bearing area of ​​the force-bearing body is respectively The unit vertex set, the internal point set on the unit edge, the internal point set on the unit face, and the internal point set of the unit for all tetrahedral units obtained in , the value of the intrinsic mixed finite element basis function at each integration point can be determined according to the following rules:

[0115] For a given , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0116]

[0117] For a given , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0118]

[0119] For a given , which is an internal point on the element face F, and F is a tetrahedral element and tetrahedral elements

[0120] The common surface of , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0121]

[0122] For a given , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0123]

[0124] For a given , which is an internal point on the element edge E, and the element E is a tetrahedral element , , …, The common edge of , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0125]

[0126] For a given , then the value of the intrinsic mixed finite element basis function at each integration point can be written as:

[0127]

[0128] Among them, the above is the value of the Lagrange basis function at the integration point, represents the k-th order Lagrangian finite element.

[0129] It should be noted that the above intrinsic mixed finite element space , are continuous; and To be discontinuous.

[0130] Step S3, determine the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and construct the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit.

[0131] In the technical solution provided in the above step S3, each unit is traversed, and the overall degrees of freedom of the finite element basis function of the unit are determined according to the unit type (i.e., Lagrangian element type, intrinsic mixed finite element type) of the unit, and the coefficient matrix A and the load vector b are constructed by combining the values ​​of the body force on the unit at each integration point in the current unit and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the current unit.

[0132] Optionally, Figure 8 is a flow chart of an optional step S3 according to an embodiment of the present application, such as Figure 8 Traverse each unit and execute the following steps S31-S33, including:

[0133] Step S31, determining the unit type of the current unit, wherein when the unit type of the current unit is a Lagrangian element type, executing step S32; when the unit type of the current unit is an intrinsic mixed finite element type, executing step S33;

[0134] Step S32, constructing the coefficient matrix A and the load vector b according to the following steps, including:

[0135] Step S321, sorting the interpolation points corresponding to the k-order Lagrangian basis function in the order of continuous degrees of freedom to obtain the overall degrees of freedom of the k-order Lagrangian basis function of the current unit;

[0136] Step S322, using the gradient value of the k-order Lagrangian basis function obtained in step S24 at each integration point to calculate the unit stiffness matrix, and combining the overall degree of freedom of the k-order Lagrangian basis function obtained in step S321, the calculated unit stiffness matrix is ​​loaded into the coefficient matrix A;

[0137] Step S323, calculate the unit load vector using the k-order Lagrangian basis function obtained in step S24 and the value of the body force at each integration point, and assemble the calculated unit load vector into the load vector b in combination with the overall degrees of freedom of the k-order Lagrangian basis function obtained in step S321.

[0138] Step S33, constructing the coefficient matrix A and the load vector b according to the following steps, including:

[0139] Step S331, determining the overall degrees of freedom of the intrinsic mixed finite element basis function of the current unit;

[0140] Step S332, using the values ​​of the intrinsic mixed finite element basis functions at each integration point obtained in step S25 to calculate the unit coefficient matrix, and combining step S331 to construct the coefficient matrix A;

[0141] Step S333, sorting the interpolation points corresponding to the k-1 order discontinuous Lagrangian basis function obtained in step S25 according to the order of discontinuous degrees of freedom, to obtain the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function of the current unit;

[0142] Step S334, using the values ​​of the k-1 order discontinuous Lagrangian basis functions at each integration point obtained in step S25 and the divergence values ​​of the intrinsic mixed finite element basis functions at each integration point to calculate the unit divergence operator matrix, and combining the overall degrees of freedom of the intrinsic mixed finite element basis functions obtained in step S331 and the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis functions obtained in step S333 to assemble the unit divergence operator matrix into the coefficient matrix A;

[0143] Step S335, calculate the unit load vector using the values ​​of the k-1 order discontinuous Lagrangian basis function at each integration point and the value of the body force at each integration point obtained in step S25, and assemble the unit load vector into the load vector b in combination with the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function obtained in step S333.

[0144] Further, Fig. 9 is a flow chart of an optional step S331 according to an embodiment of the present application, such as Fig. 9 As shown. It includes the following steps:

[0145] Step 1: Determine the type of interpolation point corresponding to the intrinsic mixed finite element basis function, where when the type of Lagrangian interpolation point is a point on the element edge or a point on the element surface, execute the second step; when the type of Lagrangian interpolation point is a unit vertex or an internal point of the unit, execute the third step;

[0146] Step 2: Determine the overall degrees of freedom of the intrinsic mixed finite element basis function according to the order of continuous degrees of freedom;

[0147] Step 3: determine whether the intrinsic hybrid finite element basis function is composed of pure tangent vectors, wherein, if the intrinsic hybrid finite element basis function is composed of pure tangent vectors, execute the fourth step, and if the intrinsic hybrid finite element basis function is not composed of pure tangent vectors, execute the fifth step;

[0148] Step 4: Determine the overall degrees of freedom of the intrinsic hybrid finite element basis function according to the discontinuous degrees of freedom sequence;

[0149] The fifth step is to determine the overall degrees of freedom of the intrinsic mixed finite element basis function according to the order of continuous degrees of freedom;

[0150] Step 6: Determine whether the loop reaches the last intrinsic mixed finite element basis function. If so, execute step 7; otherwise, execute step 1.

[0151] Step 7: Integrate the continuous degree of freedom sequence and the discontinuous degree of freedom sequence to obtain the overall degrees of freedom of the intrinsic mixed finite element basis function.

[0152] Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s.

[0153] In the technical solution provided in the above step S4, the above boundary unit is a unit for processing the edge position of the physical domain in multiple units obtained by dividing the force-bearing body. Therefore, multiple boundary units are traversed, and each boundary condition on the boundary unit is processed according to the type of boundary condition on the boundary unit (including surface force and displacement) to adjust the load vector b determined in the above step S3, and at the same time obtain the constraint equation group Gx=s, wherein G in the constraint equation group Gx=s represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector.

[0154] Optionally, Fig.10 is a flow chart of an optional step S4 according to an embodiment of the present application, such as Fig.10 The boundary conditions of each unit are traversed and the following steps S41-S43 are performed, including:

[0155] Step S41, determining the unit type of the boundary unit to which the current boundary condition belongs, wherein when the unit type is an intrinsic mixed finite element type, executing step S42, and when the unit type is a Lagrangian element type, executing step S43;

[0156] Step S42, processing the current boundary condition according to the following steps, including:

[0157] Step S421, when the current boundary condition is displacement, obtain the unit external normal vector of the current boundary, call the intrinsic hybrid finite element basis function, and load the integral result of the intrinsic hybrid finite element basis function, the unit external normal vector and the displacement into the load vector b;

[0158] Step S422, when the current boundary condition is surface force, the least square method is called to process the surface force, and the processing result is added to the constraint equation group Gx=s.

[0159] Step S43, processing the current boundary condition according to the following steps, including:

[0160] Step S431, when the current boundary condition is the surface force, calling the k-order Lagrangian basis function of the Lagrangian element, and loading the integral result of the k-order Lagrangian basis function and the surface force into the load vector b;

[0161] Step S432, when the current boundary condition is displacement, the least square method is called to process the displacement, and the processing result is added to the constraint equation group Gx=s.

[0162] Step S5, looping over the unit interface of the Lagrangian sub-type unit and the intrinsic hybrid finite element type unit, obtaining the order k1 of the Lagrangian sub-type unit, the order k2 of the intrinsic hybrid finite element type unit, and the unit normal vector from the intrinsic hybrid finite element type unit to the Lagrangian sub-type unit on the current unit interface.

[0163] In the technical solution provided in the above step S5, the purpose of using k1 and k2 is to distinguish the orders of the Lagrangian element type unit and the intrinsic hybrid finite element type unit, and their specific values ​​are still the same as the order k of the Lagrangian element basis function and the order k of the intrinsic hybrid finite element basis function. It should be noted that in actual application scenarios, the values ​​of k1 and k2 can be the same or different.

[0164] Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit.

[0165] In the technical solution provided in the above step S6, the value of the k1-order Lagrangian basis function at the current unit interface is equal to the value of the k1-order Lagrangian basis function at each integration point on the current unit interface, and the value of the k2-order intrinsic mixed finite element basis function at the current unit interface is equal to the value of the k2-order intrinsic mixed finite element basis function at each integration point on the current unit interface.

[0166] Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface.

[0167] As an optional implementation, the local coefficient matrix C of the current unit interface can be obtained according to the following steps S71-S72, including:

[0168] Step S71, calculating the product of the unit normal vector obtained in step S5 and the value of the k2-order intrinsic mixed finite element basis function obtained in step S6 at each integration point in the current unit interface;

[0169] Step S72, calculate the product of the calculated product result and the value of the k1-order Lagrangian basis function obtained in step S6 at each integration point in the current unit interface, and use the obtained product result as the local coefficient matrix C of the current interface.

[0170] Step S8, based on the overall degrees of freedom of the finite element basis function of the unit obtained in step S3 and the first degrees of freedom of the k1-order Lagrangian basis function and the second degrees of freedom of the k2-order intrinsic mixed finite element basis function of the current unit interface in step S6.

[0171] In the technical solution provided by the above step S8, the above step S3 is used to determine the overall degrees of freedom of the finite element basis functions of the Lagrangian sub-type unit and the intrinsic hybrid finite element type unit to which the unit interface belongs, wherein the first degree of freedom of the k1-order Lagrangian sub-basis function of the unit interface is determined based on the overall degrees of freedom of the finite element basis functions of the Lagrangian sub-type unit and the value of the k1-order Lagrangian sub-basis function at the current unit interface; the second degree of freedom of the k2-order intrinsic hybrid finite element basis function of the current unit interface is determined based on the overall degrees of freedom of the finite element basis functions of the intrinsic hybrid finite element type unit and the value of the k2-order intrinsic hybrid finite element basis function at the unit interface.

[0172] Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8.

[0173] In the technical solution provided in the above step S9, the local coefficient matrix C of the current unit interface is assembled into the coefficient matrix A according to the first degree of freedom and the second degree of freedom of the current unit interface.

[0174] It should be noted that, for each unit interface of the Lagrangian element type unit and the intrinsic hybrid finite element type unit, the above steps S5-S9 are repeated to assemble the local coefficient matrix C of each unit interface into the coefficient matrix A to obtain the final coefficient matrix A.

[0175] Step S10, using the constraint equation group Gx=s obtained in step S4, the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 are processed to obtain a linear algebraic equation group Ax=b.

[0176] As an optional implementation, the linear algebraic equation system Ax=b may be obtained according to the following method, including:

[0177] Step S10, simplify the stress constraint equation group Gx=s to obtain the equivalent constraint ,in, represents the constraint variable, represents the free variable, m and n represent the number of bound variables and free variables respectively, and It should be noted that in the above step S101, Gaussian elimination method can be used to simplify the stress constraint equation group Gx=s;

[0178] Step S102, adjust the order of x, and obtain the adjusted x as:

[0179]

[0180] Step S103, block processing is performed on the coefficient matrix A obtained in steps S3-S9 and the load vector b obtained in steps S3-S4, and the coefficient matrix A after block processing is obtained as follows:

[0181]

[0182] And the load vector b after block processing is:

[0183]

[0184] Step S104, rearrange the coefficient matrix A after block processing to obtain the adjusted coefficient matrix A:

[0185]

[0186] The load vector b after block processing is rearranged to obtain the adjusted load vector b:

[0187]

[0188] Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

[0189] Specifically, the linear algebraic equation Ax=b may be solved by a direct solution method or an iterative solution method, and no specific limitation is imposed on the solution method.

[0190] It should be noted that the technical solution provided in the embodiments of the present application can be widely used in the fields of vehicles, civil engineering, aerospace, etc., and its specific application scenarios include stiffness analysis of chassis frames, deformation and damage analysis of tunnels, and lightweight analysis of rocket brackets.

[0191] Based on the scheme defined by the above steps S1 to S11, it can be known that, in the embodiment, the intrinsic hybrid finite element method and the Lagrangian finite element method are combined so that a high-precision solution can still be obtained at the edge position of the region, thereby solving the technical problem that the related technology uses the traditional finite element method to solve the linear elastic mechanics problem and obtains the result with low precision.

[0192] The following will perform strength analysis on structures where there are geometric discontinuities and compare the results with those of commercial software (Abaqus or Standard) to demonstrate that the method for solving mechanical problems based on the hybrid finite element method provided in the embodiment of the present application can greatly improve the numerical stability in the iterative process, increase convergence, and improve computational efficiency.

[0193] Fig.11a is a schematic diagram of an optional L-shaped structure according to an embodiment of the present application, such as Fig.11aThe model has a geometric discontinuity at the corner with radius R=0.025m. In the finite element analysis based on displacement solution, the convergence of stress of such a geometric body with a sharp change in stress cross section is strongly correlated with the density of the grid. Fig.11b is a schematic diagram of a finite element model of an optional L-shaped structure according to an embodiment of the present application, such as Fig.11b As shown, it is not difficult to see that the maximum stress will appear at the chamfer. If the traditional displacement solution finite element analysis method is used, the mesh density at the chamfer needs to be extremely fine to obtain a converged solution.

[0194] Specifically, during the verification process, the finite element unit type used was the plane stress unit CPS (Capcom PlaySystem) 3 (also known as the 3-node triangular plane stress unit), and boundary conditions and body forces were set on each finite element unit, where the type of body force was concentrated force.

[0195] First, under the same boundary conditions and body forces, the Lagrangian finite element method is used to analyze the unit size at the chamfer, the total degree of freedom of the model, the degree of freedom at the chamfer, the maximum stress value and the stress change rate of the first-order finite element model (number 1-6) and the second-order finite element model (number 7) with different densities, and the calculation results are shown in Table 1 below.

[0196] Table 1

[0197]

[0198] At the same time, get Figure 12-18 Combined with the above Table 1 and Figure 12-18 It is not difficult to see that when the mesh size at the chamfer is fine to m, the stress change rate of the L-shaped structure has been reduced to 0.43%, so the stress result under this size can be used as the exact solution of the model. That is, when the L-shaped structure is fixedly constrained on the right end face of the model and a concentrated force of -100N in the y direction is applied to the upper end face, the maximum stress value is Pa.

[0199] In addition, for Fig.11a The L-shaped model shown in the figure uses the third-order intrinsic mixed finite element method to analyze the corners of the model (i.e., the chamfers), while the rest of the model is analyzed using the Lagrangian method, and the mesh size at the chamfers is controlled within m, obtained by commercial software analysis, the maximum stress value obtained by comprehensive analysis using this intrinsic hybrid finite element and Lagrangian element is Pa, and the stress cloud diagram is shown in Fig.19 Compared with the highest calculation accuracy in Table 1 above For example, the stress error value is 1.086%, and when the intrinsic mixed finite element is used to analyze the chamfer of the model, the model degree of freedom at the chamfer is 32369, which is higher than the highest calculation accuracy in Table 1. For the corresponding chamfered corners with 429474 degrees of freedom, the reduction is 92.46%. In other words, when solving the same engineering structure model, especially for the stress concentration problem caused by the large difference in stress section changes, the intrinsic mixed finite element can obtain an accurate solution to the structural model with a small number of degrees of freedom, which greatly saves computing resources and improves simulation accuracy.

[0200] Example 2

[0201] According to an embodiment of the present application, a non-volatile storage medium is also provided, in which a program is stored, wherein when the program is running, the device where the non-volatile storage medium is located is controlled to execute the method for solving mechanical problems based on the hybrid finite element method in Example 1.

[0202] Optionally, the device where the non-volatile storage medium is located implements the following steps by running the program:

[0203] Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body, wherein the unit data at least includes: unit shape, node coordinates, unit number, unit type, unit dimension, and the unit type includes: intrinsic mixed finite element type and Lagrangian element type, the material parameters include: Young's modulus, Poisson's ratio, density, and the boundary conditions include: surface force and displacement;

[0204] Step S2, looping each unit, obtaining multiple integration points on the current unit, and determining the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point in combination with the unit data, material parameters and body force obtained in step S1, including: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point;

[0205] Step S3, determining the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and constructing the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit;

[0206] Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s, where G represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector;

[0207] Step S5, looping each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit, obtaining the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit, and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface;

[0208] Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit;

[0209] Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface;

[0210] Step S8, determining the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the interface of the current unit based on the overall degree of freedom of the finite element basis function of the unit obtained in step S3 and step S6;

[0211] Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8;

[0212] Step S10, using the constraint equation group Gx=s obtained in step S4, processing the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain a linear algebraic equation group Ax=b;

[0213] Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

[0214] According to an embodiment of the present application, a computer program product is also provided, which includes a stored computer program, wherein when the computer program is executed by a processor, the method for solving mechanical problems based on the hybrid finite element method in Example 1 is implemented.

[0215] Optionally, the computer program executes the following steps:

[0216] Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body, wherein the unit data at least includes: unit shape, node coordinates, unit number, unit type, unit dimension, and the unit type includes: intrinsic mixed finite element type and Lagrangian element type, the material parameters include: Young's modulus, Poisson's ratio, density, and the boundary conditions include: surface force and displacement;

[0217] Step S2, looping each unit, obtaining multiple integration points on the current unit, and determining the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point in combination with the unit data, material parameters and body force obtained in step S1, including: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point;

[0218] Step S3, determining the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and constructing the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit;

[0219] Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s, where G represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector;

[0220] Step S5, looping each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit, obtaining the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit, and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface;

[0221] Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit;

[0222] Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface;

[0223] Step S8, determining the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the interface of the current unit based on the overall degree of freedom of the finite element basis function of the unit obtained in step S3 and step S6;

[0224] Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the unit interface obtained in step S8;

[0225] Step S10, using the constraint equation group Gx=s obtained in step S4, processing the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain a linear algebraic equation group Ax=b;

[0226] Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

[0227] According to an embodiment of the present application, a processor is also provided, which is used to run a program, wherein when the program is running, the method for solving mechanical problems based on the hybrid finite element method in Example 1 is executed.

[0228] Optionally, the following steps are performed when the program is running:

[0229] Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body, wherein the unit data at least includes: unit shape, node coordinates, unit number, unit type, unit dimension, and the unit type includes: intrinsic mixed finite element type and Lagrangian element type, the material parameters include: Young's modulus, Poisson's ratio, density, and the boundary conditions include: surface force and displacement;

[0230] Step S2, looping each unit, obtaining multiple integration points on the current unit, and determining the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point in combination with the unit data, material parameters and body force obtained in step S1, including: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point;

[0231] Step S3, determining the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and constructing the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit;

[0232] Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s, where G represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector;

[0233] Step S5, looping each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit, obtaining the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit, and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface;

[0234] Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit;

[0235] Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface;

[0236] Step S8, determining the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the interface of the current unit based on the overall degree of freedom of the finite element basis function of the unit obtained in step S3 and step S6;

[0237] Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8;

[0238] Step S10, using the constraint equation group Gx=s obtained in step S4, processing the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain a linear algebraic equation group Ax=b;

[0239] Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

[0240] According to an embodiment of the present application, an electronic device is also provided, wherein the electronic device includes one or more processors; a memory for storing one or more programs, so that when the one or more programs are executed by the one or more processors, the one or more processors are implemented to run the programs, wherein the program is configured to execute the method for solving mechanical problems based on the hybrid finite element method in the above-mentioned embodiment 1 when running.

[0241] Optionally, the processor is configured to implement the following steps by executing a computer program:

[0242] Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body, wherein the unit data at least includes: unit shape, node coordinates, unit number, unit type, unit dimension, and the unit type includes: intrinsic mixed finite element type and Lagrangian element type, the material parameters include: Young's modulus, Poisson's ratio, density, and the boundary conditions include: surface force and displacement;

[0243] Step S2, looping each unit, obtaining multiple integration points on the current unit, and determining the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point in combination with the unit data, material parameters and body force obtained in step S1, including: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point;

[0244] Step S3, determining the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and constructing the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit;

[0245] Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s, where G represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector;

[0246] Step S5, looping each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit, obtaining the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit, and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface;

[0247] Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit;

[0248] Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface;

[0249] Step S8, determining the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the interface of the current unit based on the overall degree of freedom of the finite element basis function of the unit obtained in step S3 and step S6;

[0250] Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8;

[0251] Step S10, using the constraint equation group Gx=s obtained in step S4, processing the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain a linear algebraic equation group Ax=b;

[0252] Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

[0253] The serial numbers of the above-mentioned embodiments of the present application are for description only and do not represent the advantages or disadvantages of the embodiments.

[0254] In the above embodiments of the present application, the description of each embodiment has its own emphasis. For parts that are not described in detail in a certain embodiment, please refer to the relevant description of other embodiments.

[0255] In the several embodiments provided in this application, it should be understood that the disclosed technical content can be implemented in other ways. Among them, the device embodiments described above are only schematic. For example, the division of the units can be a logical function division. There may be other division methods in actual implementation. For example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of units or modules, which can be electrical or other forms.

[0256] The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed on multiple units. Some or all of the units may be selected according to actual needs to achieve the purpose of the present embodiment.

[0257] In addition, each functional unit in each embodiment of the present application may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of software functional units.

[0258] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application, or the part that contributes to the relevant technology or all or part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium, including a number of instructions to enable a computer device (which can be a personal computer, server or network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage medium includes: U disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), mobile hard disk, disk or optical disk, etc. Various media that can store program codes.

[0259] The above is only a preferred implementation of the present application. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present application. These improvements and modifications should also be regarded as the scope of protection of the present application.

Claims

1. A method for solving mechanical problems based on a hybrid finite element method, characterized in that: The method comprises the following solving steps: Step S1, obtaining unit data, material parameters, boundary conditions and body forces of multiple units that divide the load-bearing body, wherein the unit data at least includes: unit shape, node coordinates, unit number, unit type, unit dimension, and the unit type includes: intrinsic mixed finite element type and Lagrangian element type, the material parameters include: Young's modulus, Poisson's ratio, density, and the boundary conditions include: surface force and displacement; Step S2, looping each unit, obtaining multiple integration points on the current unit, and determining the value of the body force at each integration point and the value, gradient value or divergence value of the finite element basis function corresponding to the unit type at each integration point in combination with the unit data, material parameters and body force obtained in step S1, including: the value of the intrinsic mixed finite element basis function at each integration point, the divergence value of the intrinsic mixed finite element basis function at each integration point, the value of the discontinuous Lagrangian basis function at each integration point, the value of the continuous Lagrangian basis function at each integration point and the gradient value of the continuous Lagrangian basis function at each integration point; Step S3, determining the overall degrees of freedom of the finite element basis function of each unit according to the unit type, and constructing the coefficient matrix A and the load vector b by combining the values ​​of the body force at each integration point in the unit obtained in step S2 and the values, gradient values ​​or divergence values ​​of the finite element basis function corresponding to the unit type at each integration point in the unit; Step S4, looping for each boundary unit, processing each boundary condition according to the type of boundary condition on the boundary unit, adjusting the load vector b, and obtaining the constraint equation group Gx=s, where G represents the constraint matrix, x represents the vector to be solved, and s represents the constraint vector; Step S5, looping each unit interface of the Lagrangian unit and the intrinsic hybrid finite element unit, obtaining the order k1 of the Lagrangian unit, the order k2 of the intrinsic hybrid finite element unit, and the unit normal vector from the intrinsic hybrid finite element unit to the Lagrangian unit on the current unit interface; Step S6, calculating the values ​​of the k1-order Lagrangian basis function and the k2-order intrinsic mixed finite element basis function at the interface of the current unit; Step S7, combining steps S5 to S6 to determine the local coefficient matrix C of the current unit interface, includes: Traverse each unit interface and perform the following steps: Step S71, calculating the product of the unit normal vector obtained in step S5 and the value of the k2-order intrinsic mixed finite element basis function obtained in step S6 at each integration point in the current unit interface; Step S72, calculating the product of the calculated product result and the value of the k1-order Lagrangian basis function at each integration point in the current unit interface obtained in step S6, and using the obtained product result as the local coefficient matrix C of the current interface; Step S8, determining the first degree of freedom of the k1-order Lagrangian basis function and the second degree of freedom of the k2-order intrinsic mixed finite element basis function of the interface of the current unit based on the overall degree of freedom of the finite element basis function of the unit obtained in step S3 and step S6; Step S9, assembling the local coefficient matrix C of the current unit interface obtained in step S7 into the coefficient matrix A based on the first degree of freedom and the second degree of freedom of the current unit interface obtained in step S8; Step S10, using the constraint equation group Gx=s obtained in step S4, processing the coefficient matrix A obtained in step S9 and the load vector b obtained in step S4 to obtain a linear algebraic equation group Ax=b; Step S11, solving the linear algebraic equation system Ax=b to obtain stress results and displacement results.

2. The method according to claim 1, characterized in that: The step S2 comprises the following steps: Traverse each cell and perform the following steps: Step S21, for each current unit, obtaining the unit volume, unit order k, multiple integration points and integration weight of each integration point of the current unit; Step S22, using the physical strength obtained in step S1, determining the value of the physical strength at each integral point on the current unit; Step S23, determining the unit type of the current unit; Step S24, when the unit type is a Lagrangian element type, using step S21 to calculate the value of the k-order Lagrangian element basis function of the current unit at each integration point and the gradient value of the k-order Lagrangian element basis function at each integration point; Step S25, when the unit type is an intrinsic hybrid finite element type, obtain the tensor basis of the three-dimensional symmetric matrix space of the current unit, and combine the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point calculated in step S21 to obtain the value of the intrinsic hybrid finite element basis function of the current unit at each integration point and the divergence value of the intrinsic hybrid finite element basis function at each integration point, and at the same time use step S21 to calculate the value of the k-1-order discontinuous Lagrangian basis function of the current unit at each integration point.

3. The method according to claim 2, characterized in that The step S25 comprises the following steps: Step S251, using step S21 to calculate the value of the k-order Lagrangian basis function at each integration point, the gradient value of the k-order Lagrangian basis function at each integration point, and the value of the k-1-order discontinuous Lagrangian basis function at each integration point; Step S252, determining the type of each interpolation point on the current cell, wherein the types of interpolation points include: cell vertices, cell internal points, points on cell edges, and points on cell faces; Step S253, determining a tensor basis of the three-dimensional symmetric matrix space according to the type of the interpolation point; Step S254, combines the tensor basis of the three-dimensional symmetric matrix space corresponding to each type of interpolation point with the value of the k-order Lagrangian basis function at each integration point and the gradient value of the k-order Lagrangian basis function at each integration point obtained in step S251, to obtain the value of the intrinsic mixed finite element basis function at each integration point and the divergence value of the intrinsic mixed finite element basis function at each integration point.

4. The method according to claim 3, characterized in that: The step S253 comprises the following steps: Step 1: When the interpolation point type is a unit vertex or a unit internal point, obtain the standard basis of the third-order symmetric tensor space; Step 2: When the interpolation point type is a point on the unit edge, obtain a unit tangent vector and two unit normal vectors for each of the six edges of the current unit, and use the unit tangent vector and two unit normal vectors of each edge to determine the edge basis of the unit edge; Step 3: When the interpolation point type is a point on a unit face, obtain two unit tangent vectors and one unit normal vector for each of the four faces of the current unit, and use the two unit tangent vectors and one unit normal vector of each face to determine the surface basis of the unit face.

5. The method according to claim 1, characterized in that The step S3 comprises the following steps: Traverse each cell and perform the following steps: Step S31, determining the unit type of the current unit; Step S32, when the element type is a Lagrangian element type, construct the coefficient matrix A and the load vector b according to the following steps, including: Step S321, sorting the interpolation points corresponding to the k-order Lagrangian basis function in the order of continuous degrees of freedom to obtain the overall degrees of freedom of the k-order Lagrangian basis function of the current unit; Step S322, using the gradient value of the k-order Lagrangian basis function obtained in step S24 at each integration point to calculate the unit stiffness matrix, and combining the overall degree of freedom of the k-order Lagrangian basis function obtained in step S321, the calculated unit stiffness matrix is ​​loaded into the coefficient matrix A; Step S323, calculating the unit load vector using the k-order Lagrangian basis function obtained in step S24 and the value of the body force at each integration point, and assembling the calculated unit load vector into the load vector b in combination with the overall degree of freedom of the k-order Lagrangian basis function obtained in step S321; Step S33, when the element type is an intrinsic mixed finite element type, construct the coefficient matrix A and the load vector b according to the following steps, including: Step S331, determining the overall degrees of freedom of the intrinsic mixed finite element basis function of the current unit; Step S332, using the values ​​of the intrinsic mixed finite element basis functions at each integration point obtained in step S25 to calculate the unit coefficient matrix, and combining step S331 to construct the coefficient matrix A; Step S333, sorting the interpolation points corresponding to the k-1 order discontinuous Lagrangian basis function obtained in step S25 according to the order of discontinuous degrees of freedom, to obtain the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function of the current unit; Step S334, using the values ​​of the k-1 order discontinuous Lagrangian basis functions at each integration point obtained in step S25 and the divergence values ​​of the intrinsic mixed finite element basis functions at each integration point to calculate the unit divergence operator matrix, and combining the overall degrees of freedom of the intrinsic mixed finite element basis functions obtained in step S331 and the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis functions obtained in step S333 to assemble the unit divergence operator matrix into the coefficient matrix A; Step S335, calculate the unit load vector using the k-1 order discontinuous Lagrangian basis function obtained in step S25 and the value of the body force at each integration point, and assemble the unit load vector into the load vector b in combination with the overall degrees of freedom of the k-1 order discontinuous Lagrangian basis function obtained in step S333.

6. The method according to claim 5, characterized in that The step S331 includes the following steps: Step 1: Determine the type of interpolation points corresponding to the intrinsic mixed finite element basis function; Step 2: When the interpolation point type is the unit vertex and the unit internal point, the overall degree of freedom of the intrinsic mixed finite element basis function is determined according to the order of continuous degree of freedom; Step 3: When the interpolation point type is a point on the element edge or a point on the element surface, determine whether the intrinsic mixed finite element basis function is composed of pure tangent vectors; Step 4: When the intrinsic hybrid finite element basis function is composed of pure tangent vectors, determine the overall degrees of freedom of the intrinsic hybrid finite element basis function according to the discontinuous degree of freedom sequence; Step 5: When the intrinsic hybrid finite element basis function is not composed of pure tangent vectors, determine the overall degrees of freedom of the intrinsic hybrid finite element basis function according to the order of continuous degrees of freedom.

7. The method according to claim 1, characterized in that The step S4 comprises the following steps: Traverse each boundary condition of each boundary cell and perform the following steps: Step S41, determining the unit type of the boundary unit to which the current boundary condition belongs; Step S42, when the element type is an intrinsic mixed finite element type, the current boundary condition is processed according to the following steps, including: Step S421, when the current boundary condition is displacement, obtain the unit external normal vector of the current boundary, call the intrinsic hybrid finite element basis function, and load the integral result of the intrinsic hybrid finite element basis function, the unit external normal vector and the displacement into the load vector b; Step S422, when the current boundary condition is surface force, the least square method is used to process the surface force, and the processing result is added to the constraint equation group Gx=s; Step S43, when the cell type is a Lagrangian element type, the current boundary condition is processed according to the following steps, including: Step S431, when the current boundary condition is the surface force, calling the k-order Lagrangian basis function of the Lagrangian element, and loading the integral result of the k-order Lagrangian basis function and the surface force into the load vector b; Step S432, when the current boundary condition is displacement, the least square method is called to process the displacement, and the processing result is added to the constraint equation group Gx=s.

8. A computer program product, characterized in that include: A computer program, wherein when the computer program is executed by a processor, the method for solving mechanical problems based on the hybrid finite element method as described in any one of claims 1 to 7 is implemented.

9. An electronic device, characterized in that: include: A memory and a processor, wherein the processor is used to run a program stored in the memory, wherein when the program is run, the method for solving mechanical problems based on the hybrid finite element method as described in any one of claims 1 to 7 is executed.

Images