Hybrid heat flux finite element method based on octree grid

Abstract

The application relates to the technical field of numerical simulation, in particular to a hybrid heat flux finite element method based on an octree grid. In a traditional temperature finite element method, a post-processing derivation method is usually used for heat flux calculation, which leads to the fact that the accuracy of the heat flux calculation result is lower than that of the temperature calculation result. According to the scheme, the heat flux is taken as an independent variable, and a heat flux function is constructed, so that higher-order solutions can be obtained in the solution of the temperature field and the heat flux field, and the accuracy of the heat flux calculation is improved. In the application, the octree grid is used for grid division, and the characteristic of the octree grid is that adaptive division can be carried out according to the complexity and requirements of a problem. This means that the grid can be more finely divided in the region of interest, so that the calculation accuracy is improved. Meanwhile, the hanging nodes of the octree grid do not have any limitation on the calculation of the heat flux, so that the calculation of the heat flux is more accurate.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology, and more specifically to the hybrid heat flux finite element method based on octree mesh. Background Technology

[0002] High-temperature environments can cause thermal deformation and thermal stress in materials, leading to performance degradation and reduced reliability. Furthermore, for high-performance particle-reinforced composite materials used in spacecraft, heat transfer is also affected by the reinforcement, causing disturbances in the material's temperature field. This can result in uncoordinated thermal deformation, thermal stress concentration, and other adverse consequences, further impacting spacecraft reliability. Therefore, calculating the temperature and overall heat flux of particle-reinforced composite materials is of great significance.

[0003] In computational mechanics, numerical simulation methods such as the finite element method (FEM) are widely used for analyzing heat conduction problems due to their simplicity and stability. Zienkiewicz first applied the FEM to steady-state thermal analysis in 1965. The temperature-based FEM uses temperature as the variable, directly calculating nodal temperatures, while heat flux is calculated using post-processing differentiation. This inevitably reduces the accuracy of heat flux calculations, resulting in a lower accuracy (one order lower) than temperature calculations for the same problem. Furthermore, for particle-reinforced composites, very fine meshes are required to obtain meaningful results for both temperature and heat flux, leading to low computational efficiency. However, in engineering, solving for the equivalent thermal conductivity of particle-reinforced composites requires simultaneously solving for both the temperature and heat flux fields. Several improved methods have been proposed to address these difficulties, but they still cannot solve problems related to composite materials.

[0004] Extending the two-dimensional hybrid heat flux finite element method to 3D is a challenging task due to differences in geometric features and solution schemes in existing technologies. This invention proposes a 3D hybrid heat flux finite element method based on a polyhedral octree. The octree mesh is particularly suitable for implementation using HF-FEM because the presence of dangling nodes does not impose any limitations on HF-FEM. HF-FEM assumes independent heat flux and temperature fields within the element domain and on the element surface, respectively. A thermal functional is derived, and the heat flux field pre-satisfies the thermal balance governing equations, thereby constructing a universal heat flux function. The heat flux parameter and temperature are variationally derived using the functional, and then the same solution equation format as traditional finite element method is derived simultaneously to obtain the nodal temperature of the heat flux parameter. When assembling the matrix in the calculation of a single rigid element, HF-FEM requires integration both within the element domain and on the element surface, which necessitates selecting the integration scheme based on the shape of the decomposed region within the element domain. To facilitate the computation of integrals in complex domains, arbitrary polyhedral elements are tetrahedralized, and Hammer integrals are used to solve for the simple domains. The results for the complex domains are then summed. This invention develops an octree-based 3DHF-FEM algorithm for steady-state heat conduction analysis of spherical inclusion composite materials.

[0005] The study of steady-state heat conduction in particle-reinforced composite materials using the hybrid finite element method is currently limited to two-dimensional problems. Traditional finite element algorithms use regularly shaped elements, which cannot well fit the structural shapes in practical problems. This invention proposes a hybrid heat flux finite element method using hybrid arbitrary polyhedral elements, and further proposes a 3D octree hybrid heat flux finite element method for composite materials containing spherical inclusions to solve steady-state heat conduction. Summary of the Invention

[0006] The purpose of this invention is to provide a hybrid heat flux finite element method based on octree meshes. Compared with the traditional temperature finite element method, the proposed method can achieve significantly higher orders in both the temperature and heat flux fields. Specifically, when solving problems with large gradients in the temperature and heat flux fields, the proposed algorithm can obtain solutions that meet accuracy requirements with a minimal number of elements, improving computational efficiency and enabling the solution of large-scale heat conduction problems that traditional methods cannot handle.

[0007] To achieve the above-mentioned technical objectives and effects, the present invention is implemented through the following technical solution: The finite element method for hybrid heat flux based on octree mesh includes the following steps: S1: Mesh generation: Construct an octree structure and divide it into different mesh cells.

[0008] S2: Mesh Adaptation: Adapt the octree mesh to the complexity and requirements of the problem to make the mesh more refined.

[0009] S3: Perform finite element analysis of three-dimensional hybrid heat flux distribution.

[0010] Furthermore, step S1 specifically includes: S1.1: Use the smallest outer cube that can completely enclose the target region as the initial model to divide the region, that is, to obtain an octree root node unit; S1.2: Divide the root element into 8 sub-elements based on the maximum element size, and then check the node type in turn to further divide the elements.

[0011] S1.3: Based on the octree unit partitioning logic, the partially mixed cube units are recursively divided into 8 equal-sized sub-units. That is, a partially mixed cube is subdivided into 8 smaller sub-units. Partitioning stops when a unit contains only mixed elements or matrix. S1.4: Further divide the cells according to the 2:1 balance principle to obtain a balanced octree mesh.

[0012] Furthermore, step S2 includes smoothing the octree mesh and optimizing the calculation of the G matrix; Furthermore, the smoothing process of the octree mesh specifically includes: S2.1: Locate the octree element on the boundary between the matrix material and the inclusion material, obtained from the inclusion material control equation. Substitute the element node coordinates into the inclusion material control equation to determine whether the nodes on the element edge are inside the inclusion material, thereby identifying the required element.

[0013] S2.2: Find the coordinates of the intersection point between the inclusion and the matrix material in the found boundary element. Iterate through the element, determining whether the edge intersects with the inclusion material (i.e., the intersection point of the straight line and the inclusion material's governing equation). Substitute the coordinates of the starting and ending points on the element edge into the inclusion governing equation: (1.1) Where x, y, z are coordinates, and the starting point coordinates are P0(x 0, y 0, ,z 0, The endpoint coordinates are P1(x1,y1,z1) and C0(x c ,y c ,z cLet be the center coordinates of the inclusion material, and a, b, and c be the three semi-principal axes of the inclusion ellipsoid. Substituting the starting and ending point coordinates into the equation, if one result is greater than 1 and the other less than 1, then there is an intersection point on that edge. When the result is exactly equal to 1, the intersection point is exactly at a node. To find the coordinates of the intersection point on an edge with an intersection point, the inclusion control equation (Equation 1.1) and the equation with the edge as a straight line (Equation 1.2) must be solved simultaneously, i.e.: (1.2) Where x, y, z are the coordinates of points on the line, x0, y0, z0 are the coordinates of the starting point, m, n, p are the line equation vectors, and t is the coordinate parameter.

[0014] The starting point coordinates are P0(x) 0, y 0, ,z 0, The endpoint coordinates are P1(x1,y1,z1), and the vectors m, n, and p are respectively: (1.3) Combining these equations, we obtain the following equation: (1.4) Solving equation (1.4) yields the intersection parameter t: (1.5) in: (1.6) (1.7) According to the generality of the solution, equation (1.4) has two solutions. Since we need to obtain the points inside the intersection, we discard the solutions where t is greater than 1. Finally, we substitute the obtained t back into equation (1.2) to obtain the coordinates of the intersection point.

[0015] S2.3: Remove duplicate point coordinates. Since faces in a unit share common edges, looping along these common edges will find duplicate intersection points, which need to be deleted.

[0016] S2.4: The calculated intersection points are used as the boundary points of the new surface. A new overall number is assigned to each boundary point, and the surface is divided into two based on this number: one for the mixed material elements and one for the matrix material elements. This forms the basis for dividing the mesh into two elements. The common tangent of these two elements is the smooth transition surface. Orderly connecting the numbering information of the transition surfaces is a key step in mesh transition processing, specifically including: First, find the coordinates of the center point of the common cross plane: (1.8) Where ct(x,y,z) represents the coordinates of the center point of the common tangent plane, (x i ,y i ,z i ) represents the coordinates of each vertex on the common facet, and h represents the number of vertices on the common facet.

[0017] Next, calculate the h vectors formed by each vertex and the center point of the common tangent plane: (1.9) in, Let x be the vector formed by the i-th vertex and the center point. i y i z i Let x be the coordinate of the i-th vertex. ct y ct z ct The coordinates of the center point.

[0018] Using the first vector as the reference direction vector, solve for h-1 cosine values ​​obtained from the first vector and the remaining vectors respectively: (1.10) in, Let cosine be the sum of the first and j-th vectors. , Let the first point and the j-th vector be... , Let be the modulus of the first point and the j-th vector.

[0019] Continue by cross-multiplying the first vector with the other vectors to obtain h-1 normal vectors: (1.11) Where j is the vector index, and j takes values ​​from 2 to h. Let the first vector be the new vector obtained by the cross product of the first vector and the j-th vector. Let be the vector of the j-th vertex; The vector pointing from the center of the inclusion to the center of the common tangent plane is defined as the reference normal vector: (1.12) in, Let the vector be the center of the interleaving point to the center of the common tangent. The normal vector result obtained from equation (1.11) is compared with the reference normal vector. Dot product: (1.13) Among them, I j The new vector is obtained by dot product of Equation 11 and the reference normal vector.

[0020] Based on the result of the dot product, when I j Greater than 0 indicates that it is related to The normal vector of the dot product is in the same direction as the reference normal vector, and vice versa. j Less than 0 indicates that... The normal vector of the dot product is in the opposite direction to the reference normal vector. This yields the ordered arrangement information on the common tangent plane, and the common tangent plane is placed into two cells.

[0021] S2.5: Remove duplicate nodes. After all boundary elements have been processed, duplicate node information may still appear between elements. This duplicate information needs to be removed. Two simple mapping relationships are used. The node number information in the mapping relationship is modified based on whether the node coordinates corresponding to the node number are duplicated. Then, the non-duplicate nodes are recounted to obtain the correct node numbers and related information. This results in a smooth octree mesh.

[0022] Furthermore, the optimization calculation of the G matrix specifically includes: The formula for calculating the G matrix is ​​as follows: (1.14) in, Represents the integration domain of the element surface, P T Represents the transpose of the heat flux polynomial coordinate matrix, n T M represents the transpose of the out-of-plane normal vector of the element, and M represents the interpolation matrix on the element boundary.

[0023] The integral formula on any polygon is: (1.15) Where S represents an arbitrary polygonal domain, P T Represents the transpose of the heat flux polynomial coordinate matrix, n T M represents the transpose of the out-of-plane normal vector of the element, and M represents the interpolation matrix on the polygon boundary.

[0024] The coordinates of the geometric center of any polygon can be obtained by interpolating the coordinates of each point. Similarly, the temperature of the geometric center of a polygon can also be obtained by interpolating the temperatures of each point of the polygon.

[0025] (1.16) in, Represents the temperature at the geometric center of the polygon. The value represents the temperature of the polygon node, and n represents the number of polygon nodes.

[0026] If we form a small triangle with two adjacent nodes and the center of any polygon, then the temperature at any point within the small triangle can be obtained by analogy from Equation 2.19. (1.17) in, For any point inside the triangle, For the temperature of node i in the polygon element, The temperature of node i+1 of the polygon element. Temperature at the center of the polygon, The interpolation weights of node i in the polygon cell. The interpolation weights for the element node i+1 on the polygon. The interpolation weights are the centers of the polygon.

[0027] Substituting equation 1.16 into equation 1.17, we get (1.18) That is, the M matrix can be optimized as follows: (1.19) It can be observed that Equation 1.18 only applies to the current calculation. and An additional term is added at the corresponding position of each node, merely changing the form of M. This eliminates the need to add points at the center of the element surface, while still ensuring the correctness of the integral of the G matrix and maintaining the original size of the G matrix. Therefore, Equation 1.15 can be rewritten as follows: (1.20) Among them, S i P represents the integration domain of the small triangles in the polygonal domain. T Represents the transpose of the heat flux polynomial coordinate matrix, n T M represents the transpose of the out-of-plane normal vector of the element, and M represents the interpolation matrix on the element boundary.

[0028] Further, step S3 includes: assigning correct boundary conditions to each mesh element, discretizing the heat flux and heat conduction equations into finite element equations, and then solving them to obtain the nodal temperature of each mesh element; obtaining the heat flux distribution by solving for the unknown parameters of the heat flux field based on the nodal temperatures; and finally integrating the heat flux field to obtain the temperature field distribution. Further, step S3 specifically includes the following sub-steps: S3.1: Functional Derivation; There is a volume in space. Limited enclosed area its boundary Including unit boundaries Heat flux boundary and temperature boundary The heat flux within the domain is from It indicates that the temperature is from This indicates that the heat flux along the outer normal direction n of the boundary is... express.

[0029] The governing equations for steady-state heat conduction within the domain are: Heat balance equation (2.1) in, For the gradient of heat flux, It is a volumetric heat source.

[0030] Fourier's law of heat conduction (2.2) make ,in, Here is the thermal conductivity matrix. Since the thermal conductivity matrix is ​​the inverse matrix, the above equation can also be expressed as: (2.3) Boundary conditions include: Forced boundary conditions (2.4) in, For temperature, For a given temperature, Given a temperature boundary.

[0031] Natural boundary conditions (2.5) in, This represents the heat flux in the direction normal to the boundary. , For a given heat flux, This is the transpose of the unit vector of the outward normal direction of the boundary. Given a temperature boundary.

[0032] By weakening Fourier's law of heat conduction and the first kind of boundary conditions using the Galerkin method, the equivalent integral form is written as follows: (2.6) in, For the transpose of the heat flux variation, For heat flux vector, This represents the temperature gradient.

[0033] Applying Gauss's formula to the first term yields... (2.7) in, The transpose of the unit vector of the outward normal direction of the boundary. It is a variation of the gradient of the transpose of the heat flux vector.

[0034] At the boundary of a given heat flux Therefore, the gradient variation of heat flux is 0, that is... (2.8) in, It is the transpose of the unit vector of the outward normal direction of the boundary.

[0035] Therefore, equation (2.6) can be written as (2.9) in, Given a temperature.

[0036] Variational construction of functionals (2.10) That is, the functional of this problem is (2.11) Abbreviated as (2.12) Consider the second type of boundary condition (2.5) and the continuity of heat flux between elements; (2.13) in, and These represent the outward normal vectors on the shared element surface of two adjacent elements. This represents the transpose of the heat flux vector on the common surface of two adjacent elements.

[0037] By applying the Lagrange multiplier method to introduce the above constraints, a modified functional is obtained.

[0038] (2.14) S3.2: Construction of the heat flux function; In the finite element method, a self-balancing heat flux function for each element is introduced. To facilitate subsequent differentiation and integration operations, the heat flux field is assumed to be a polynomial represented by an interpolation matrix and unknown interpolation coefficients. As the order of the polynomial function increases, the computational accuracy can be improved.

[0039] That is, assume the heat flux function is a complete polynomial.

[0040] in, Let be the heat flux function. To provide complete polynomial coefficients; However, in order to make the stiffness matrix invertible, we take...

[0041] in, The polynomial coefficients are reselected to make the stiffness matrix invertible.

[0042] The heat flux must satisfy equation (2.1), and based on this, we can set three heat flux components as follows: (2.15) Represented in matrix form, it is: (2.16) S3.3: Solving for the temperature interpolation function; After being divided, each element consists of multiple tetrahedrons, and therefore the surface of each element is composed of several triangles. The boundary temperature of the element can then be obtained by interpolating the nodal temperatures. For example... Figure 1 Depend on , , The area of ​​the spatial triangle formed by the three nodes is The temperature value at its node is , , One point inside the face The three vertices can each form three triangles. The area coordinates are represented as (2.17) It has the following properties (2.18) Therefore, the area coordinates of a triangle are chosen as the shape function to interpolate the element boundary temperature. Then, a point within the surface... Temperature can be expressed as (2.19) in, Let p be the temperature at the interpolation point. , , This is the interpolation function corresponding to the triangle point.

[0043] Written in matrix form as (2.20) in, Let p be the temperature at the interpolation point. Temperature interpolation function matrix , Nodal temperature matrix .

[0044] S3.4: Unit scaling; Obtain the centroid coordinates of the element, and scale the element with the centroid as the origin of the local coordinate system; divide the element into sections in the Cartesian coordinate system. The unit center is obtained by dividing a Delaunay tetrahedron into its constituent parts, calculating the center coordinates of each tetrahedron, and then indirectly calculating the unit center through the tetrahedron center.

[0045] Assume the coordinates of the four vertices of the tetrahedron are (2.23) The coordinates of the tetrahedron's center can be calculated. (2.24) If the division The coordinates of the center of each tetrahedron are: The volumes are respectively From this, the center coordinates of any three-dimensional polyhedron can be obtained. .

[0046] (2.25) Once the center of any polyhedral element is obtained, element scaling can be performed. Any point on the element... After translating the origin of the coordinate system to the center, divide by the scaling factor. Coordinates in the local coordinate system can be obtained .

[0047] (2.26) S3.5: Solve for the heat flux field; heat flux components and boundary temperature By introducing a functional in matrix form, we can obtain... (2.27) in, This is the transpose of the coefficient matrix. The coefficient matrix, This is the transpose of a polynomial matrix. It is a polynomial matrix. This is the transpose of the unit vector of the outward normal direction of the boundary. Given a heat flux.

[0048] make , , The above formula can be simplified to (2.28) Based on the stationarity condition of the functional (2.29) achievable (2.30) in, , Let be an integral matrix. The heat source load matrix is ​​obtained by integration.

[0049] After sorting, we can get (2.31) A generalized temperature control equation can be obtained. (2.32) in, For the heat conduction matrix, , The node temperature matrix, Here is the nodal thermal load matrix.

[0050] Once the node temperatures are obtained, the heat flux parameters can be calculated. Then The heat flux field can be solved by substituting back into equation (2.16).

[0051] The beneficial effects of this invention are: This invention, based on the field of numerical simulation technology, proposes a hybrid heat flux finite element method based on an octree mesh, primarily used to solve the steady-state heat conduction problem of particle-reinforced composite materials. Compared with the traditional temperature finite element method, this invention has the following specific advantages: 1. High-precision heat flux calculation: In traditional temperature finite element methods, heat flux calculations typically employ post-processing differentiation, resulting in a lower-order accuracy of the heat flux calculation compared to the temperature calculation. This invention, by treating heat flux as an independent variable and constructing a heat flux function, can obtain higher-order solutions in both the temperature and heat flux fields, thereby improving the accuracy of heat flux calculation.

[0052] 2. Adaptable Octree Mesh: This invention uses an octree mesh for mesh generation. The octree mesh is characterized by its adaptability to the complexity of the problem and specific requirements. This means that the mesh can be finely divided within the region of interest, thereby improving computational accuracy. Furthermore, the dangling nodes of the octree mesh do not impose any limitations on heat flux calculations, leading to more accurate heat flux calculations.

[0053] 3. Improved computational efficiency: For complex systems such as particle-reinforced composites, traditional methods require very fine meshes to obtain meaningful temperature and heat flux results, leading to low computational efficiency. The octree mesh in this invention can be adaptively partitioned according to the complexity of the problem, thereby improving computational efficiency while maintaining computational accuracy. Furthermore, this invention employs tetrahedral partitioning and Hammer integration methods to easily calculate integrals in complex domains, further enhancing computational efficiency.

[0054] 4. Application of Polyhedral Octree Mesh: The octree mesh in this invention is particularly suitable for the hybrid heat flux finite element method because the dangling nodes of the octree mesh do not impose any restrictions on the calculation of heat flux. Furthermore, in order to calculate the integral of the complex domain, this invention divides arbitrary polyhedral elements into tetrahedral parts and uses Hammer integration to solve the simple domain, finally summing the results to obtain the result for the complex domain.

[0055] In summary, this invention, by introducing octree meshes and the hybrid heat flux finite element method, can solve the steady-state heat conduction problem of particle-reinforced composite materials, and achieves significant improvements in the accuracy and efficiency of heat flux calculation. This is of great significance for addressing the problems of material performance degradation and reliability reduction under high-temperature environments, and improves the reliability and service life of high-performance materials such as those used in spacecraft.

[0056] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0057] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0058] Figure 1 This is a schematic diagram of the triangular surface described in an embodiment of the present invention; Figure 2 This is a schematic diagram of the Voronoi diagram and Delaunay triangular duality properties as described in the embodiments of the present invention; Figure 3 This is a schematic diagram of the Delaunay partitioning of the three-dimensional polyhedral unit according to an embodiment of the present invention; Figure 4 This is a schematic diagram showing the natural coordinates (a) of any point in the triangle and (b) of any point in the tetrahedron according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the octree mesh generation process described in an embodiment of the present invention; Figure 6 This is a schematic diagram of the sawtooth boundary (a) and the smooth boundary (b) described in an embodiment of the present invention; Figure 7 This is a schematic diagram of the interface triangle connection information according to an embodiment of the present invention; Figure 8 This is a schematic diagram of the optimized G matrix integration according to an embodiment of the present invention; Figure 9 This is a schematic diagram and result image of the homogeneous irregular unit model described in an embodiment of the present invention; Figure 10 This is a schematic diagram and result image of the homogeneous thick-walled cylindrical model described in an embodiment of the present invention; Figure 11 This is a schematic diagram and result diagram of a centrally located mixed-use example model as described in an embodiment of the present invention; Figure 12 This is a schematic diagram and result image of the corner point containing two inclusions as described in the embodiment of the present invention; Figure 13 This is a schematic diagram and result diagram of a model containing multiple mixed models as described in an embodiment of the present invention. Detailed Implementation

[0059] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] Example 1

[0061] like Figure 1-8 As shown The hybrid heat flux finite element method based on octree mesh described in this embodiment includes the following steps: S1: Mesh generation: Construct an octree structure and divide it into different mesh cells as required.

[0062] S2: Mesh Adaptation: Adapt the octree mesh to the complexity and requirements of the problem to make the mesh more refined.

[0063] S3: Perform finite element analysis of three-dimensional hybrid heat flux distribution.

[0064] In this embodiment, step S1 specifically includes: S1.1: Use the smallest outer cube that can completely enclose the target region as the initial model to divide the region, that is, to obtain an octree root node unit; S1.2: Divide the root element into 8 sub-elements based on the largest element size, such as... Figure 5 As shown in (b), the node types are then checked sequentially to further divide the units.

[0065] S1.3: Based on the octree unit partitioning logic, the partially mixed cube units are recursively divided into 8 equal-sized sub-units. That is, a partially mixed cube is subdivided into 8 smaller sub-units. Partitioning stops when a unit contains only mixed elements or matrix. S1.4: Further divide the units according to the 2:1 balance principle to obtain... Figure 5 The balanced octree mesh in (c).

[0066] In this embodiment, step S2 includes smoothing the octree mesh and optimizing the calculation of the G matrix; In this embodiment, the smoothing process of the octree mesh specifically includes: S2.1: Locate the octree elements at the boundary between the matrix material and the inclusion material. This step can be obtained from the inclusion material governing equations. Substitute the element node coordinates into the inclusion material governing equations to determine whether the nodes on the element edge are within the inclusion material, thereby identifying the required element.

[0067] S2.2: Find the coordinates of the intersection point between the inclusion and the matrix material within the found boundary elements. Iterate through the elements, determining if the edge intersects with the inclusion material, i.e., the intersection point of the line with the inclusion material's governing equation. Substitute the coordinates of the starting and ending points on the element edge into the inclusion governing equation: (1.1) in The center coordinates of the inclusion material, These are the three semi-principal axes of the mixed ellipsoid. Substituting the coordinates of the starting and ending points into this equation, if one result is greater than 1 and the other less than 1, then there is an intersection point on that edge. When the result is exactly equal to 1, the intersection point is exactly at a node. To find the coordinates of the intersection point on an edge with an intersection point, it is necessary to solve the mixed control equation (Equation 1.1) and the equation with the edge as a straight line (Equation 1.2) simultaneously, that is: (1.2) The starting point coordinates are The coordinates of the end point are ,vector They are respectively: (1.3) Combining these equations, we obtain the following equation: (1.4) Solving equation (1.4) yields the intersection parameter t: (1.5) in: (1.6) (1.7) According to the generality of the solution, equation (1.4) has two solutions. Since we need to obtain the points inside the intersection, we discard the solutions where t is greater than 1. Finally, we substitute the obtained t back into equation (1.2) to obtain the coordinates of the intersection point.

[0068] S2.3: Remove duplicate point coordinates. Since faces in a unit share common edges, looping along these common edges will find duplicate intersection points, which need to be deleted.

[0069] S2.4: The calculated intersection points are used as the boundary points of the new surface. A new overall number is assigned to each boundary point, and the surface is divided into two based on this number: one for the mixed material elements and one for the matrix material elements. This forms the basis for dividing the mesh into two elements. The common tangent of these two elements is the smooth transition surface. Orderly connecting the numbering information of the transition surface is a key step in the mesh transition processing.

[0070] First, find the coordinates of the center point of the common cross plane: (1.8) Next, calculate the h vectors formed by each vertex and the center point of the common tangent plane: (1.9) Using the first vector as the reference direction vector, solve for h-1 cosine values ​​obtained from the first vector and the remaining vectors respectively: (1.10) Continue by cross-multiplying the first vector with the other vectors to obtain h-1 normal vectors: (1.11) The vector pointing from the center of the inclusion to the center of the common tangent plane is defined as the reference normal vector: (1.12) The normal vector result obtained from equation (1.11) is compared with the reference normal vector. Dot product: (1.13) Based on the result of the dot product, when Greater than 0 indicates that it is related to The normal vector of the dot product is in the same direction as the reference normal vector, and vice versa. Less than 0 indicates that... The dot product of the normal vector is opposite to the reference normal vector. Since the ordering information of the node numbers is related to the cosine value between them, and the angles between vectors (acute and obtuse) can result in positive or negative values, a reference normal vector is introduced. This reference normal vector is then used to perform a dot product with the result of the previous cross product to determine the node order. Nodes with normal vectors in the same direction are numbered from largest to smallest using bubble sort, and nodes with normal vectors in opposite directions are numbered from smallest to largest using bubble sort as well. Then, nodes are connected in order from those with normal vectors in the same direction to those with opposite directions, and the connection information of the element faces is stored in the direction of the outward normal according to the right-hand rule. At this point, the ordered arrangement information on the common tangent plane is obtained, and the common tangent plane is placed into two elements.

[0071] S2.5: Remove duplicate nodes. After all boundary elements have been processed, duplicate node information may still appear between elements. It is necessary to remove the duplicate node information. Here, two simple mapping relationships are used. The numbering information in the mapping relationship is modified according to whether the node coordinates corresponding to the node number are duplicated. Then, the nodes that are not duplicated are recounted to obtain the correct node number and related information.

[0072] The above steps will yield a smooth octree mesh.

[0073] Furthermore, the optimization calculation of the G matrix specifically includes: The formula for calculating the G matrix is ​​as follows: (1.14) The integral formula on any polygon is: (1.15) The coordinates of the geometric center of any polygon can be obtained by interpolating the coordinates of each point. Similarly, the temperature of the geometric center of a polygon can also be obtained by interpolating the temperatures of each point of the polygon.

[0074] (1.16) like Figure 8 If we form a small triangle with the center of any polygon by taking any two adjacent nodes, the temperature at any point within the small triangle can be obtained by analogy from Equation 2.19. (1.17) Substituting equation 1.16 into equation 1.17, we get (1.18) That is, the M matrix can be optimized as follows: (1.19) It can be observed that Equation 1.18 only applies to the current calculation. and An additional term is added at the corresponding position of each node, merely changing the form of M. This eliminates the need to add points at the center of the element surface, while still ensuring the correctness of the integral of the matrix G and maintaining the original size of the G matrix. Therefore, Equation 1.15 can be rewritten as follows: (1.20) In this embodiment, step S3 includes: assigning correct boundary conditions to each grid cell, discretizing the heat flux and heat conduction equations into finite element equations, and then solving them to obtain the temperature distribution of each grid cell; obtaining the required temperature distribution based on the temperature field data through methods such as interpolation or fitting; and calculating the heat flux distribution according to the temperature distribution and the heat conduction equations.

[0075] In this embodiment, step S3 specifically includes the following sub-steps: S3.1: Functional Derivation; There is a volume in space. Limited enclosed area its boundary Including unit boundaries Heat flux boundary and temperature boundary The heat flux within the domain is from It indicates that the temperature is from Indicates the direction of the outer normal of the boundary. The heat flux on is from express.

[0076] The governing equations for steady-state heat conduction within the domain are: Heat balance equation (2.1) Fourier's law of heat conduction (2.2) make The above formula can also be expressed as (2.3) Boundary conditions include: Forced boundary conditions (2.4) Natural boundary conditions (2.5) This represents the heat flux in the direction normal to the boundary. , Let be the unit vector of the direction of the outward normal to the boundary.

[0077] By weakening Fourier's law of heat conduction and the first kind of boundary conditions using the Galerkin method, the equivalent integral form is written as follows: (2.6) Applying Gauss's formula to the first term yields... (2.7) For a given heat flux Therefore, the gradient variation of heat flux is 0, that is... (2.8) Therefore, equation (2.6) can be written as (2.9) Variational construction of functionals (2.10) That is, the functional of this problem is (2.11) Abbreviated as (2.12) Consider the second type of boundary condition (2.5) and the continuity of heat flux between elements; (2.13) By applying the Lagrange multiplier method to introduce the above constraints, a modified functional is obtained.

[0078] (2.14) S3.2: Construction of the heat flux function; In the finite element method, a self-balancing heat flux function for each element is introduced. To facilitate subsequent differentiation and integration operations, the heat flux field is assumed to be a polynomial represented by an interpolation matrix and unknown interpolation coefficients. As the order of the polynomial function increases, the computational accuracy can be improved.

[0079] That is, assume the heat flux function is a complete polynomial.

[0080] in, To provide complete polynomial coefficients; However, in order to make the stiffness matrix invertible, we take...

[0081] The heat flux must satisfy equation (2.1), and based on this, we can set three heat flux components as follows: (2.15) Represented in matrix form, it is: (2.16) S3.3: Solving for the temperature interpolation function; After being divided, each element consists of multiple tetrahedrons, and therefore the surface of each element is composed of several triangles. The boundary temperature of the element can then be obtained by interpolating the nodal temperatures. For example... Figure 1 Depend on , , The area of ​​the spatial triangle formed by the three nodes is The temperature value at its node is , , One point inside the face The three vertices can each form three triangles. The area coordinates are represented as (2.17) It has the following properties (2.18) Therefore, the area coordinates of a triangle are chosen as the shape function to interpolate the element boundary temperature. Then, a point within the surface... The temperature can be expressed as .

[0082] (2.19) Written in matrix form as (2.20) in (2.21) (2.22) S3.4: Unit scaling; To address the issue of non-invertibility of the H matrix caused by excessively large or small numerical values ​​during computation, and to avoid significant errors in the calculation results, it is necessary to scale the elements to obtain a local coordinate system for numerical calculation. Therefore, it is necessary to obtain the centroid coordinates of the elements and scale them using the centroid as the origin of the local coordinate system. However, the elements in the hybrid element method generally have complex shapes, making it impossible to directly determine their centroids. This paper proposes dividing the elements into Cartesian coordinates. The unit center is obtained by dividing a Delaunay tetrahedron into its constituent parts, calculating the center coordinates of each tetrahedron, and then indirectly calculating the unit center through the tetrahedron center.

[0083] Assume the coordinates of the four vertices of the tetrahedron are (2.23) The coordinates of the tetrahedron's center can be calculated. (2.24) If the division The coordinates of the center of each tetrahedron are: The volumes are respectively From this, the center coordinates of any three-dimensional polyhedron can be obtained. .

[0084] (2.25) Once the center of any polyhedral element is obtained, element scaling can be performed. Any point on the element... After translating the origin of the coordinate system to the center, divide by the scaling factor. Coordinates in the local coordinate system can be obtained .

[0085] (2.26) S3.5: Solve for the heat flux field; The heat flux components derived above and boundary temperature By introducing a functional in matrix form, we can obtain... (2.27) make , , The above formula can be simplified to (2.28) Based on the stationarity condition of the functional (2.29) achievable (2.30) After sorting, we can get (2.31) A generalized temperature control equation can be obtained. (2.32) in, .

[0086] Once the node temperatures are obtained, the heat flux parameters can be calculated. Then The heat flux field can be solved by substituting back into equation (2.16).

[0087] Example 2

[0088] like Figure 8-13 As shown To evaluate the performance and practicality of the derived elements, the calculation results were compared with those of traditional finite element methods. To investigate the performance of the proposed elements in solving heat conduction and thermal stress problems, the following five models were calculated: In Example 2.1, the heat conduction problem of a homogeneous hexahedral model divided by three arbitrary polyhedral elements was calculated to verify the correctness of the proposed method; in Example 2.2, the heat conduction problem of a thick-walled cylindrical model with an analytical solution was calculated, and the calculated temperature results were compared with the analytical solution to verify the calculation accuracy of HF-FEM; in Examples 2.3-2.5, steady-state heat conduction problems of a model with one central inclusion, a model with two corner inclusions, and a model with five random inclusions were calculated, respectively, to evaluate the performance of the two types of elements. The material parameters for the above five models are as follows.

[0089]

[0090] 2.1 Homogeneous Irregular Units The first example aims to verify the correctness of the proposed method. The simplest way to verify this is to calculate a uniform heat flux field model. Therefore, consider the following... Figure 9 The regular hexahedral model shown in (a) is divided into three irregular arbitrary polyhedral elements. After Delaunay tetrahedral subdivision, it is as follows: Figure 9 As shown in (b), the Abaqus model is as follows: Figure 9 As shown in (c). The model has a side length of 2 cm. The end face temperature of =0 is 30℃. =2 has an end face temperature of 50 ℃.

[0091] In Abaqus, the same model parameter size was used, the element type was set to 8-node linear heattransfer brick, and the mesh size was set to the default approximate global size, resulting in a total of 1000 elements and 1331 nodes.

[0092] The x-direction heat flux contour plots calculated by HF-FEM and Abaqus are shown below. Figure 9 (d) and (e), temperature cloud maps as follows Figure 9 (g) and (h). Numerical comparisons were performed on 50 points along the path from (0, 1, 1) to (2, 1, 1) in the fields of the two physical quantities mentioned above, as shown below. Figure 9 As shown in (f) and (i).

[0093] It can be seen that the heat flux and temperature calculated by HFFEM and Abaqus are highly consistent. Furthermore, the error between the two methods along the chosen path is negligible. This demonstrates the correctness and effectiveness of the proposed method.

[0094] 2.2 Homogeneous Thick-Walled Cylinder To verify the numerical accuracy of the HF-FEM calculation results, this example calculates the heat conduction problem of a thick-walled cylindrical model with an inner diameter of 1.3 dm, an outer diameter of 2 dm, and a height of 0.6 dm, and compares the results with analytical solutions and Abaqus calculations. The model is as follows: Figure 10 As shown in (a), the Delaunay tetrahedral partition diagram of a single element is as follows: Figure 10 As shown in (b), the Abaqus model is as follows: Figure 10 As shown in (c). The boundary conditions are an inner wall temperature of 50°C and an outer wall temperature of 30°C.

[0095] In Abaqus, the same model parameter size was used, the element type was set to 8-node linear heattransfer brick, and the mesh generation used the default approximate global size, resulting in a total of 4200 elements and 5600 nodes.

[0096] From the analytical solution of the heat conduction temperature of a thick-walled cylinder, the temperature along the radial direction is... (3.1) in, The inner wall temperature The outer wall temperature, Let the inner radius be , The outer radius is This is the distance from the selected point to the center of the circle.

[0097] The x-direction heat flux contour plots calculated by HF-FEM and Abaqus are shown below. Figure 10 (d) and (e), temperature cloud maps as follows Figure 10 (g) and (h). Fifty points in the temperature field along the path (1.3, 0.3, 0) to (2, 0.3, 0) were used for numerical comparison, as shown below. Figure 10 As shown in (f) and (i).

[0098] It can be seen that the heat flux and temperature contour maps calculated by HF-FEM are highly consistent with those calculated by Abaqus. The temperature values ​​calculated by the HF-FEM method are in good agreement with the Abaqus and analytical solutions. Moreover, the HF-FEM method only used 120 elements to obtain results that are very close to both, which shows that the method has the characteristics of high efficiency and high accuracy.

[0099] 2.3 The center contains a mixed example. This example is as follows: Figure 11 (a) shows a regular hexahedron model with a side length of 2 cm, containing a spherical inclusion with a radius of 0.5 cm at its center (1, 1, 1). The octree mesh diagram is shown below. Figure 11As shown in (b), the Abaqus model diagram is as follows: Figure 11 As shown in (c). The boundary conditions are: =2 has an end face temperature of 50 ℃. The end face temperature of =0 is 30 ℃.

[0100] In Abaqus, the same model parameter size was used, the element type was set to 10-node quadratictetrahedron, and the mesh generation used the default approximate global size, resulting in a total of 53,158 elements and 75,663 nodes.

[0101] The heat flux in the x-direction calculated by HF-FEM and Abaqus is in Cloud map of the cross section as follows Figure 11 (d) and (e), temperature cloud maps as follows Figure 11 (g) and (h). The values ​​of the two physical quantity fields are compared at 50 points along the path from (0, 1, 1) to (2, 1, 1), as shown below. Figure 11 As shown in (f) and (i).

[0102] The calculation results show that the results of the two methods are very close, but the octree mesh only uses 952 elements and 1311 nodes. This indicates that the HF-FEM and HS-FEM methods can obtain satisfactory results with fewer elements.

[0103] 2.4 Corner Point with Two Mixed Examples This example is as follows: Figure 12 (a) shows a regular hexahedral base with a side length of 2 cm. At the two corner points (2, 2, 0) and (2, 2, 2), there are spherical inclusions with a radius of 0.5 cm, each being one-eighth of a sphere. The octree mesh diagram is shown below. Figure 12 As shown in (b), the Abaqus model diagram is as follows: Figure 12 As shown in (c). The end face temperature of =2 is 50 ℃. The end face temperature is 30 ℃.

[0104] In Abaqus, using the same model parameter size, the element type was set to 10-node quadratictetrahedron, resulting in a total of 8249 elements and 12421 nodes. The octree mesh consisted of 1156 elements and 1727 nodes.

[0105] The x-direction heat flux contour plots calculated by HF-FEM and Abaqus are shown below. Figure 12 (d) and (e), temperature cloud maps as follows Figure 12(g) and (h). Temperature values ​​at 50 points along the path from (0, 2, 0) to (2, 2, 0) are extracted and compared, such as... Figure 12 (i) Extract the heat flux values ​​in the x-direction from 50 points along the path from (2, 2, 0) to (2, 2, 2) for comparison, such as... Figure 12 As shown in (f).

[0106] from Figure 12 As can be seen from (d) and (e), there is a large heat flux gradient at the hybrid junction between the matrix and the sandwich element. Thanks to the advantages of the hybrid element method, the proposed element can represent this change with only a small number of elements. This not only demonstrates the advantages of HF-FEM compared to traditional finite element method, but also reflects the improvement in computational efficiency of the proposed method.

[0107] 2.5 Contains multiple mixed models This example considers the following: Figure 13 (a) shows a regular hexahedral model with a side length of 4 cm. It contains 5 randomly located inclusions at coordinates (1, 1, 1), (1, 2, 3), (2, 3, 1), (2.5, 1, 3), and (2, 2, 2), with an inclusion radius of 0.25 cm. The octree mesh is shown below. Figure 13 As shown in (b), the Abaqus model diagram is as follows: Figure 13 As shown in (c). Boundary settings The end face temperature of =4 is 50 ℃. The end face temperature of =0 is 0 ℃.

[0108] In Abaqus, using the same model parameters and dimensions, with the element type set to 10-node quadratictetrahedron and the mesh generation using the default approximate global size, a total of 341,738 elements and 474,179 nodes were created. The octree mesh consisted of 925 elements and 1,288 nodes.

[0109] The results calculated by HF-FEM and Abaqus are as follows: The heat flux contour plot in the x-direction of the cross section is as follows Figure 13 (d) and (e), temperature cloud maps as follows Figure 13 (g) and (h). The heat flux values ​​in the x-direction at 50 points along the path from (0, 0, 0) to (4, 4, 4) are extracted and compared, as shown below. Figure 13 As shown in (f), temperature values ​​at 50 points along the path from (0, 1.25, 2.2) to (4, 1.25, 2.2) are extracted and compared. Figure 13 (i).

[0110] Although this example only has 5 inclusions, the actual number of inclusions in a calculated model can be much larger. The results show that the proposed method achieves results largely consistent with Abaqus using a significantly smaller number of elements compared to traditional methods. The temperature contour plot reflects the perturbation of the overall temperature field caused by inclusions, but the HF-FEM method, with a smaller number of elements, also captures this perturbation, demonstrating its significant advantages in efficiency and accuracy when calculating models with a large number of inclusions.

[0111] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A finite element method for hybrid heat flux based on octree meshes, characterized in that, To address the steady-state heat conduction problem in particle-reinforced composite materials, the following steps are included: S1: Mesh generation: Construct an octree structure and divide it into different mesh units; S2: Mesh Adaptation: Smoothing and G-matrix optimization calculations are performed on the octree mesh to adapt the mesh to computational requirements; S3: Perform three-dimensional hybrid heat flux finite element calculation to obtain the heat flux distribution; The three-dimensional hybrid heat flux finite element calculation refers to: treating heat flux as an independent variable, assuming independent heat flux field and temperature field in the element domain and on the element surface respectively, and solving the heat conduction problem by constructing a thermal energy functional and finding stationary values; The smoothing process for the octree mesh specifically includes: S2.1: Locate the octree element on the boundary between the matrix material and the inclusion material. Obtain it from the inclusion material control equation. Substitute the element node coordinates into the inclusion material control equation to determine whether the nodes on the element edge are inside the inclusion material, thereby determining the required element. S2.2: Find the coordinates of the intersection point between the inclusion and the matrix material in the found boundary element. Iterate through the element, determining whether the edge intersects with the inclusion material (i.e., the intersection point of the straight line and the inclusion material's governing equation). Substitute the coordinates of the starting and ending points on the element edge into the inclusion governing equation: （1.1） Where the starting point coordinates are P0(x0, y0, z0), the ending point coordinates are P1(x1, y1, z1), and C0(x c y c , z c Let be the center coordinates of the inclusion material, and a, b, and c be the three semi-principal axes of the inclusion ellipsoid. Substituting the coordinates of the starting and ending points into the equation, if one result is greater than 1 and the other less than 1, then there is an intersection point on that edge. When the result is exactly equal to 1, the intersection point is exactly at a node. To find the coordinates of the intersection point on an edge with an intersection point, it is necessary to solve the inclusion control equation and the equation with the edge as the straight line simultaneously, i.e.: （1.2） Where x, y, z are the coordinates of points on the line, x0, y0, z0 are the coordinates of the starting point, m, n, p are the line equation vectors, and t is the coordinate parameter; Vectors m, n, and p are respectively: （1.3） Where x1, y1, z1 are the coordinates of the end point of the element edge, and x0, y0, z0 are the coordinates of the start point of the element edge; Combining these equations, we obtain the following equation: （1.4） Solving equation (1.4) yields the intersection parameter t: （1.5） in: （1.6） （1.7） According to the generality of the solution, equation (1.4) has two solutions. Since we need to obtain the points inside the intersection, we discard the solutions where t is greater than 1. Finally, we substitute the obtained t back into equation (1.2) to obtain the coordinates of the intersection point. S2.3: Clear duplicate point coordinates; Since the faces in a unit have common edges, when looping along these common edges, duplicate intersection points will be found. At this time, duplicate points need to be deleted. S2.4: The calculated intersection points are used as the boundary points of the new surface. A new overall number is assigned to each boundary point, and the surface is divided according to this number, splitting the original surface into two: one for the mixed material elements and one for the matrix material elements. This forms the basis for dividing the mesh into two elements. The common tangent of these two elements is the smooth transition surface. Orderly connecting the numbering information of the transition surface is a key step in mesh transition processing, specifically including: First, find the coordinates of the center point of the common cross plane: （1.8） Where ct(x, y, z) are the coordinates of the center point of the common tangent plane, (x i y i , z i Let ) represent the coordinates of each vertex on the common face, i be the vertex index of the common face, i ranges from 1 to h, and h is the number of vertices on the common face; Next, calculate the h vectors formed by each vertex and the center point of the common tangent plane: （1.9） in, Let x be the vector formed by the i-th vertex and the center point. i y i z i Let x be the coordinate of the i-th vertex. ct y ct z ct The coordinates of the center point; Using the first vector as the reference direction vector, solve for h-1 cosine values ​​obtained from the first vector and the remaining vectors respectively: （1.10） Where j is the vector index, and j takes values ​​from 2 to h. Let cosine be the sum of the first and j-th vectors. , Let be the vector between the 1st vertex and the jth vertex. , Let be the magnitude of the first point and the j-th vector; Continue by cross-multiplying the first vector with the other vectors to obtain h-1 normal vectors: （1.11） Where j is the vector index, and j takes values ​​from 2 to h. Let the new vector be the result of the cross product of the first vector and the j-th vector. Let be the vector of the j-th vertex; The vector pointing from the center of the inclusion to the center of the common tangent plane is defined as the reference normal vector: （1.12） in, Let the vector be the center of the interleaving point to the center of the common tangent. The normal vector result obtained from equation (1.11) is compared with the reference normal vector. Dot product: （1.13） Among them, I j The new vector is obtained by dot product of the vector obtained in Equation 1.11 and the reference normal vector; Based on the dot product result, when I j Greater than 0 indicates that it is related to The normal vector of the dot product is in the same direction as the reference normal vector, and vice versa. j Less than 0 indicates that... The normal vector of the dot product is opposite to the reference normal vector; the ordered arrangement information on the common tangent is obtained, and the common tangent is placed into two units; S2.5: Remove duplicate nodes; After all boundary elements have been processed, duplicate node information will still appear between elements. It is necessary to remove duplicate node information. Two simple mapping relationships are used. The number information in the mapping relationship is modified according to whether the node coordinates corresponding to the node number are duplicated. Then, the nodes that are not duplicated are recounted to obtain the correct node number and related information; that is, a smooth octree mesh is obtained. Step S3 specifically includes the following sub-steps: S3.1: Functional Derivation; There is a volume in space. Limited enclosed area its boundary Including unit boundaries Heat flux boundary and temperature boundary The heat flux within the domain is from It indicates that the temperature is from This indicates that the heat flux along the outer normal direction n of the boundary is... express; The governing equations for steady-state heat conduction within the domain are: Heat balance equation (1.2) in, For the gradient of heat flux, As a volumetric heat source; Fourier's law of heat conduction (2.2) in, Thermal conductivity matrix The Middle Line number Thermal conductivity components of the column ; make ,in, Here is the thermal conductivity matrix. The inverse of the thermal conductivity matrix is ​​given by the above equation, which can also be expressed as: (2.3) Boundary conditions include: Forced boundary conditions (2.4) in, For temperature, For a given temperature, Given a temperature boundary; Natural boundary conditions (2.5) in, This represents the heat flux in the direction normal to the boundary. , For a given heat flux, This is the transpose of the out-of-plane normal vector of the element. Given a temperature boundary; By weakening Fourier's law of heat conduction and the first kind of boundary conditions using the Galerkin method, the equivalent integral form is written as follows: (2.6) in, For the transpose of the heat flux variation, For heat flux vector, For the temperature gradient; Applying Gauss's formula to the first term yields... (2.7) in, The variation of the gradient of the transpose of the heat flux vector; At the boundary of a given heat flux Therefore, the gradient variation of heat flux is 0, that is... (2.8) Therefore, equation (2.6) can be written as (2.9) in, Given a temperature; Variational construction of functionals (2.10) That is, the functional of this problem is (2.11) Abbreviated as (2.12) Consider the second type of boundary condition (2.5) and the continuity of heat flux between elements; (2.13) Where, n + and n - Let q represent the outward normal direction vectors of two adjacent elements on the common plane. T This is the transpose of the heat flux vector on the common surface of two adjacent units; By applying the Lagrange multiplier method to introduce the above constraints, a modified functional is obtained. (2.14) S3.2: Construction of the heat flux function; In the finite element method, a self-balancing heat flux function for each element is introduced. To facilitate subsequent differentiation and integration operations, the heat flux field is assumed to be a polynomial represented by an interpolation matrix and unknown interpolation coefficients. As the order of the polynomial function increases, the calculation accuracy is improved. That is, assume the heat flux function is a complete polynomial. in, To provide complete polynomial coefficients; However, in order to make the stiffness matrix invertible, we take... in, The polynomial coefficients were reselected to make the stiffness matrix invertible; The heat flux must satisfy equation (2.1), and based on this, let the three heat flux components be... (2.15) Represented in matrix form, it is: (2.16) S3.3: Solving for the temperature interpolation function; The divided elements are composed of multiple tetrahedrons, therefore the surface of each element is composed of several triangles, and the boundary temperature of the element is obtained by interpolation of the nodal temperatures; , , The area of ​​the spatial triangle formed by the three nodes is The temperature value at its node is , , One point inside the face Form three triangles with each of the three vertices. The area coordinates are represented as (2.17) Where n is the vertex index of the triangle, and n takes the values ​​i, j, and k respectively; L n It has the following properties: (2.18) Using the area coordinates of a triangle as the shape function, interpolation is performed on the element boundary temperature; then, a point within the surface... Temperature is expressed as (2.19) in, Let p be the temperature at the interpolation point. , , The interpolation function corresponding to the triangle point; Written in matrix form as (2.20) in, Let p be the temperature at the interpolation point. Temperature interpolation function matrix , -Nodal temperature matrix ; S3.4: Unit scaling; Obtain the centroid coordinates of the element, and scale the element with the centroid as the origin of the local coordinate system; divide the element into sections in the Cartesian coordinate system. Each Delaunay tetrahedron is divided into several parts. The center coordinates of the center of each tetrahedron are then calculated, and the unit center is indirectly calculated through the center of the tetrahedron. Assume the coordinates of the four vertices of the tetrahedron are (2.21) The coordinates of the tetrahedron's center were calculated. (2.22) If the division The coordinates of the center of each tetrahedron are: The volumes are respectively Thus, the center coordinates of any three-dimensional polyhedron are obtained. ; (2.23) After obtaining the center of any polyhedral element, the element is scaled; any point on the element... After translating the origin of the coordinate system to the center, divide by the scaling factor. Obtain coordinates in the local coordinate system ; (2.24) S3.5: Solve for the heat flux field; heat flux components and boundary temperature Introducing the functional into the matrix form, we obtain (2.25) in, This is the transpose of the coefficient matrix. The coefficient matrix, This is the transpose of a polynomial matrix. It is a polynomial matrix. This is the transpose of the outer normal vector of the element boundary. Given a heat flux; make , , The above formula simplifies to (2.26) Based on the stationarity condition of functionals (2.27) have to (2.28) in, , Let be the integral matrix. The heat source load matrix is ​​obtained by integration; After sorting (2.29) The generalized temperature control equation is obtained. (2.29) in, For the heat conduction matrix, , The node temperature matrix, Here is the nodal thermal load matrix; After obtaining the node temperature, calculate the heat flux parameters. Then Substituting back into equation (2.16) solves for the heat flux field.

2. The hybrid heat flux finite element method based on octree mesh as described in claim 1, characterized in that: Step S1 specifically includes: S1.1: Use the smallest outer cube that completely surrounds the target region as the initial model to divide the region, that is, to obtain an octree root node unit; S1.2: Divide the root element into 8 sub-elements based on the maximum element size, and then check the node type in turn to further divide the elements; S1.3: Based on the partitioning logic of the octree unit, the cubic unit containing some inclusions is recursively divided into 8 sub-units of equal size; that is, a cubic unit containing some inclusions is subdivided into 8 smaller sub-units; if a unit contains only inclusions or matrix, the partitioning stops. S1.4: Further divide the cells according to the 2:1 balance principle to obtain a balanced octree mesh.

3. The hybrid heat flux finite element method based on octree mesh as described in claim 1, characterized in that: The optimization calculation of the G matrix specifically includes: The formula for calculating the G matrix is ​​as follows: （1.14） in, Represents the integration domain of the element surface, P T Represents the transpose of the heat flux polynomial coordinate matrix, n T M represents the transpose of the out-of-plane normal vector of the element, and M represents the interpolation matrix on the element boundary; The integral formula on any polygon is: （1.15） Where S represents an arbitrary polygonal domain, and M represents the interpolation matrix on the polygonal boundary; The coordinates of the geometric center point of any polygon are obtained by interpolating the coordinates of each point. Similarly, the temperature of the geometric center of a polygon is obtained by interpolating the temperatures of each point of the polygon. （1.16） in, Represents the temperature at the geometric center of the polygon. The temperature of the polygon node is represented by n, and the number of polygon nodes is represented by n. If we form a small triangle with two adjacent nodes of any polygon and its center, then the temperature at any point within the small triangle can be obtained by interpolating the values ​​of the three points: （1.17） in, For any point inside the triangle, For the temperature of element node i on the polygon, The temperature of node i+1 of the polygon element. Temperature at the center of the polygon, The interpolation weights of node i in the polygon cell, The interpolation weights for the element node i+1 on the polygon. Interpolation weights for the polygon center; Substituting equation 1.16 into equation 1.17, we get （1.18） That is, the optimized M matrix is ​​written as （1.19） Equation 1.15 rewritten as （1.20） Among them, S i P represents the integration field of the small triangle in the polygonal domain. T Represents the transpose of the heat flux polynomial coordinate matrix, n T M represents the transpose of the out-of-plane normal vector of the element, and M represents the interpolation matrix on the element boundary. i This represents the boundary interpolation matrix corresponding to the integration domain of the i-th small triangle; Step S3 includes: assigning correct boundary conditions to each grid cell, discretizing the heat flux and heat conduction equations into finite element equations, and then solving them to obtain the nodal temperature of each grid cell; solving for the unknown parameters of the heat flux field based on the nodal temperature to obtain the heat flux distribution; and finally integrating the heat flux field to obtain the temperature field distribution.

Images