Heat energy storage heat recovery rate evaluation method based on finite element analysis

Abstract

The invention relates to the technical field of computer aided design, in particular to a heat energy storage heat recovery rate evaluation method based on finite element analysis, which comprises the following steps: importing a CAD geometric model of a heat energy storage system, analyzing the topological relation, shape attribute and component naming information of the geometric model, and acquiring original data of geometric characteristics. According to the method, a functional area system is constructed based on original data of geometric features, classification operation of a fluid domain, an energy storage material domain and a heat transfer channel is completed in a structural analysis stage, and an initial binding relation between structural attributes and functional attributes is formed. In the material parameter and grid division processing process, directional configuration of material thermophysical property parameters is achieved through structural function mapping, meanwhile, heterogeneous grid parameters and partition control precision and density are set for a fluid boundary layer and a phase change area, and the resolution capacity for details of the heat transfer process is enhanced.

Description

Technical Field

[0001] The present invention relates to the technical field of computer-aided design, and particularly to a method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis. Background Art

[0002] Computer-Aided Design (CAD for short) is a comprehensive technology that uses a computer system to perform product structure modeling, function simulation, performance analysis, and generation of manufacturing drawings.

[0003] In the prior art, material properties are often set with a unified parameter to cover all regions, ignoring the differences in local heat transfer mechanisms, which limits the modeling accuracy of heat-sensitive regions such as boundary layers and phase change regions. The application of boundary conditions is simplified based on the model shape, and the joint determination process of the structural thermal physical properties and direction characteristics is not introduced, affecting the physical consistency of the thermal simulation results. Therefore, improvements are needed. Summary of the Invention

[0004] The object of the present invention is to solve the drawbacks existing in the prior art, and to propose a method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis.

[0005] To achieve the above object, the present invention adopts the following technical solutions. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis includes the following steps:

[0006] Import the CAD geometric model of the thermal energy storage system, analyze the topological relationship, shape attributes, and component naming information of the geometric model, and obtain the original geometric feature data; based on the original geometric feature data, classify and identify the fluid domain, solid domain, energy storage material domain, heating channel, cooling channel, and insulation layer functional regions, and establish a TES geometric feature set;

[0007] Based on the TES geometric feature set, retrieve the material properties from the predefined library and assign them to the corresponding functional regions, set the mesh division parameters for the fluid boundary layer and the phase change region to obtain the material and mesh configuration parameters, and based on the material and mesh configuration parameters and the channel features identified by the TES geometric feature set, set the inlet and outlet boundary condition types to generate an FEA preprocessing configuration file;

[0008] Read the finite element analysis result file targeted by the FEA preprocessing configuration file, extract the temperature distribution data and the heat flux density vector data to obtain the original FEA result data; based on the original FEA result data, map the temperature distribution data and the heat flux density vector data to the surface mesh nodes of the original CAD model to obtain the surface heat flux density gradient field;

[0009] Based on the TES geometric feature set, determine the energy input surface, energy output surface, and the interface of the energy storage material domain. Combine the mapped finite element analysis results, calculate the evaluation value of the heat recovery rate of thermal energy storage, and superimpose and display the surface heat flux density gradient field on the 3D CAD model through color mapping to obtain the evaluation result of the heat recovery rate of thermal energy storage.

[0010] Preferably, the steps for obtaining the original geometric feature data are as follows:

[0011] Import the CAD geometric model of the thermal energy storage system, call the topological relationship information, component shape parameter information, and component naming structure information, generate the geometric boundary layer structure based on the mapping of the geometric body connection node structure, and divide the functional component combinations according to the naming structure to obtain the topological structure and component naming boundary association information;

[0012] According to the topological structure and component naming boundary association information, identify the structural connections of adjacent geometric bodies on each boundary layer, perform geometric contour fitting of node and edge combinations, and perform triangulation calculation on the connection area in the order of functional combinations to generate the boundary feature dot matrix and the set of triangular elements in the functional area;

[0013] According to the boundary feature dot matrix and the set of triangular elements in the functional area, sequentially extract the geometric dimension parameters, centroid coordinate positions, side length distributions, and curvature change information of each triangular element to generate the original geometric feature data.

[0014] Preferably, the steps for obtaining the TES geometric feature set are as follows:

[0015] According to the original geometric feature data, extract the node coordinate sequence, boundary contour curvature value, normal direction distribution, boundary node density, and regional surface area of all closed boundary regions. Group the boundary regions into an independent candidate region set based on the structural connection relationship to obtain the preliminary regional structure division data;

[0016] According to the preliminary regional structure division data, calculate the morphological distribution dispersion index of each candidate region. The calculation formula is:

[0017]

[0018] where D z is the morphological distribution dispersion index, A i is the surface area of the i-th candidate region, L i is the boundary node density of the i-th candidate region, is the average normal vector of the i-th candidate region, is the boundary center vector of the i-th candidate region, θ i is the maximum boundary normal angle of the i-th candidate region, and k is the number of candidate regions;

[0019] According to the morphological distribution dispersion index, a set of candidate regions with different structural attributes is screened, and combined with the normal vector direction distribution and node density change trend of each candidate region, functional region classification and labeling are carried out to generate a TES geometric feature set.

[0020] Preferably, the steps for obtaining the material and mesh configuration parameters are as follows:

[0021] According to the TES geometric feature set, the predefined material identifiers corresponding to each functional region label are matched, and the thermal conductivity, specific heat capacity, density, heat conduction structure number, and surface smoothness level are extracted to generate a regional material property parameter set;

[0022] According to the regional material property parameter set, the local grid resolution control value is calculated, and the calculation formula is:

[0023]

[0024] where R m is the local grid resolution control value, k t is the regional thermal conductivity, μ is the surface smoothness level, c p is the specific heat capacity, ρ is the density, γ is the heat conduction structure number, σ s is the maximum surface area of the region, λ s is the boundary gradient change ratio;

[0025] According to the local grid resolution control value, the local grid resolution control value is assigned to the corresponding functional region, and the grid accuracy of the fluid boundary layer and the grid shape factor of the phase change region are respectively labeled to form the material and mesh configuration parameters.

[0026] Preferably, the steps for obtaining the FEA preprocessing configuration file are as follows:

[0027] According to the channel features identified by the material and mesh configuration parameters and the TES geometric feature set, the geometric boundaries corresponding to the heat flow inlet and heat flow outlet in each functional region are located, and based on the boundary node normal vector direction, node distribution density, and geometric region type, the identification and label mapping of the inlet unit and the outlet unit are carried out to generate the fluid channel boundary node attribute information;

[0028] According to the fluid channel boundary node attribute information, the surface area, normal vector distribution, and grid cell type corresponding to each inlet unit and outlet unit are extracted, and combined with the material thermophysical properties and grid cell distribution form in the material and mesh configuration parameters, the boundary condition types of the inlet and outlet in each functional region are set, including the dominant mode and action dimension of the heat flux boundary or the velocity boundary, to generate the boundary condition setting information;

[0029] According to the boundary condition setting information, call the meshing details and material thermodynamics properties of each unit area in the material and mesh configuration parameters, construct the node, element, boundary term structure and constraint type syntax content of the finite element preprocessing input field, and generate the FEA preprocessing configuration file.

[0030] Preferably, the steps for obtaining the original FEA result data are as follows:

[0031] Read the finite element analysis result file targeted by the FEA preprocessing configuration file, retrieve the result index label corresponding to the analysis task, locate the temperature distribution result path and heat flux density vector result path in the output file through the mapping field, and generate a set of finite element result target paths;

[0032] According to the set of finite element result target paths, sequentially load the temperature field and heat flux density field information within each time step or iteration step, extract all temperature values and heat flux vector values according to the node numbers, and uniformly organize the extraction results into a structured matrix format to generate a temperature distribution matrix and a heat flux density vector matrix;

[0033] According to the temperature distribution matrix and the heat flux density vector matrix, compare the mesh structure and node numbering rules set in the FEA preprocessing configuration file, perform matrix node order verification and boundary matching detection, process invalid numerical units and perform node structure binding to generate the original FEA result data.

[0034] Preferably, the steps for obtaining the surface heat flux density gradient field are as follows:

[0035] According to the original FEA result data, extract the surface mesh node index and geometric coordinate information of the original CAD model, and make the node coordinates correspond one by one with the temperature distribution data and the heat flux density vector data through the index to generate a set of combined attributes of the mesh node thermophysical and geometric structures;

[0036] According to the set of combined attributes of the mesh node thermophysical and geometric structures, calculate the heat flux density gradient intensity value of each node. The calculation formula is:

[0037]

[0038] where G x is the heat flux density gradient intensity value, is the heat flux density vector of the i-th node, is the temperature gradient vector of the i-th node, is the node position vector, is the node normal vector, AG i is the mesh area where the node is located, φ i is the angle between the node heat flux vector and the normal vector;

[0039] According to the heat flux density gradient intensity value, value interpolation reconstruction is performed on each grid node in the order of the numbers, and in combination with the geometric patch topology connection structure, a node heat flux gradient isoline distribution map is constructed. Field data filling is performed at the spatial positions mapped according to the grid node coordinates on the surface of the original CAD model to generate a surface heat flux density gradient field.

[0040] Preferably, the steps for obtaining the evaluation result of the heat recovery rate of the thermal energy storage are as follows:

[0041] According to the TES geometric feature set, the energy input surface, the energy output surface and the interface of the energy storage material domain are identified. The boundary edges are marked through node indexing and the grid patch topology structure and the facing directions are bound. At the same time, the heat flux density vector and the node temperature field data in the original FEA result data are called to form a target surface thermophysical mapping structure;

[0042] According to the target surface thermophysical mapping structure, the evaluation value of the heat recovery rate of the thermal energy storage is calculated, and the calculation formula is:

[0043]

[0044] where Ψ is the evaluation value of the heat recovery rate of the thermal energy storage, is the surface heat flux density vector, is the unit vector of the surface normal direction, τ s is the node temperature change rate, ξ s is the average temperature gradient of the interface of the energy storage material domain, S in is the energy input surface, S out is the energy output surface, and t0 and t1 are the start and end boundaries of the preset time window;

[0045] According to the evaluation value of the heat recovery rate of the thermal energy storage, in combination with the node heat flux direction distribution map in the surface heat flux density gradient field, the gradient direction and the heat recovery rate correlation layer are superimposed on the grid node coordinate system of the 3D CAD model through color mapping to generate the evaluation result of the heat recovery rate of the thermal energy storage.

[0046] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0047] Based on the functional region system constructed from the original geometric feature data, the present invention completes the classification operations of the fluid domain, energy storage material domain, and heat transfer channels during the structural analysis stage, forming an initial binding relationship between the structural attributes and functional attributes. During the material parameter and mesh generation process, the directional configuration of the material's thermal physical properties is achieved through structural-functional mapping. At the same time, heterogeneous mesh parameters are set for the fluid boundary layer and phase change regions, and the accuracy and density are controlled by region to enhance the resolution ability of the details of the heat transfer process. In the boundary condition setting process, the material thermal parameters and geometric direction features are nested, and the constraint type and action region are selected according to the boundary direction and local heat conduction behavior to improve the response accuracy of the boundary loading method. The node-level mapping operation of the heat flux density and temperature distribution data is performed based on the grid node structure for one-to-one mapping, providing a basis for subsequent spatial gradient calculations. The integration calculation process takes the node heat flux direction and the energy flow change trend within the time window as inputs, dynamically obtains the total input and output energy values, and thus derives the heat recovery rate value. This value is visually mapped on the structural surface, forming a spatial linkage between the numerical results and the heat flow path, realizing the state presentation on the 3D model, and enhancing the model backtracking and heat efficiency change recognition capabilities. In the overall process, the improvement of the structural identification granularity, data organization method, and energy evaluation logic constructs a simulation closed-loop of structure-parameter-time domain integration, enabling the heat response results to have the accurate expression ability for structural feedback. Brief Description of the Drawings

[0048] Figure 1 It is a schematic diagram of the steps of the present invention. Detailed Embodiment

[0049] In order to make the objectives, technical solutions, and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0050] Please refer to Figure 1 , the present invention provides a technical solution, a method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis, including the following steps:

[0051] Import the CAD geometric model of the thermal energy storage system, analyze the topological relationship, shape attributes, and component naming information of the geometric model, and obtain the original geometric feature data; based on the original geometric feature data, classify and identify the functional regions of the fluid domain, solid domain, energy storage material domain, heating channel, cooling channel, and insulation layer, and establish a TES geometric feature set;

[0052] Based on the TES geometric feature set, retrieve the material properties from a predefined library and assign them to the corresponding functional regions, set the meshing parameters for the fluid boundary layer and the phase change region, obtain the material and mesh configuration parameters, and based on the material and mesh configuration parameters and the channel features identified by the TES geometric feature set, set the types of inlet and outlet boundary conditions to generate the FEA preprocessing configuration file;

[0053] Read the finite element analysis result file targeted by the FEA preprocessing configuration file, extract the temperature distribution data and the heat flux density vector data to obtain the original FEA result data; based on the original FEA result data, map the temperature distribution data and the heat flux density vector data to the surface mesh nodes of the original CAD model to obtain the surface heat flux density gradient field;

[0054] Based on the TES geometric feature set, determine the energy input surface, the energy output surface, and the interface of the energy storage material domain, combine the mapped finite element analysis results, calculate the evaluation value of the thermal recovery rate of the thermal energy storage, and superimpose and display the surface heat flux density gradient field on the 3D CAD model by color mapping to obtain the evaluation result of the thermal recovery rate of the thermal energy storage.

[0055] The steps for obtaining the original geometric feature data are as follows:

[0056] Import the CAD geometric model of the thermal energy storage system, call the topological relationship information, the component shape parameter information, and the component naming structure information, generate the geometric boundary layer structure based on the mapping of the geometric body connection node structure, and divide the functional component combinations according to the naming structure to obtain the topological structure and the boundary association information of the component naming;

[0057] According to the topological structure and the boundary association information of the component naming, identify the structural connections of the adjacent geometric bodies on each boundary layer, perform the geometric contour fitting of the node and edge combinations, and perform the triangulation calculation on the connection regions in the order of the functional combinations to generate the boundary feature dot matrix and the functional region triangular element set;

[0058] According to the boundary feature dot matrix and the functional region triangular element set, sequentially extract the geometric dimension parameters, the centroid coordinate positions, the side length distributions, and the curvature change information of each triangular element to generate the original geometric feature data.

[0059] Specifically, based on the imported CAD geometric model of the thermal energy storage system, combined with the topological relationship information, component shape parameter information and component naming structure information obtained in the previous stage, these data are respectively compared in coordinate dimension and attribute dimension, and the association matching process is performed according to the label content of the specific item. For example, for the shape parameter marked as "fluid transport component", it is necessary to first confirm its node coordinate distribution in the CAD model. If it is found that the distribution is consistent with the component category description in the naming information, the component is deemed to belong to the fluid transport category, and its name information is bidirectionally bound to the shape parameter. Then, the information of all components is processed in a loop. During this period, if it is detected that the data is incomplete or the matching degree with the naming structure is less than 95%, it is necessary to check for possible tag conflicts. According to experience, if the matching degree is less than 70%, it is set The situation requires additional confirmation. These matching thresholds are obtained by referring to the data collection experience of the same type of system in the past and after weighted calculation. The specific weighted calculation method is to take the average of the similarities between the component names and component shape parameters in multiple historical projects and deduct a 5% safety factor, thereby forming the currently used 70% and 95% matching thresholds. Subsequently, after all legal bindings are completed, a hierarchical analysis is performed based on the coordinate topology of the node structure connected by the geometric body, and the boundary layer structure corresponding to the component name is generated one by one. During this period, it is necessary to arrange them in order according to the degree of coordinate proximity, and through simple adjacent node judgment, for example, when the node number difference is within five and the coordinate difference is within 0.1, it is regarded as the same level. Finally, all functional components are classified and combined according to the naming structure to obtain the topological structure and component naming boundary association information.

[0060] According to the topological structure and component naming boundary association information obtained previously, when identifying the structural connections of adjacent geometric bodies on each boundary layer, first construct a corresponding list of nodes and edges. The nodes are assigned unique indices and their XYZ values in the spatial coordinate system are recorded. The edges represent their connection relationships according to the indices of two adjacent nodes. Then, in accordance with the set minimum adjacent determination criterion, for example, a direct connection relationship is considered to exist only when the spatial distance between two nodes is less than 0.05 and both belong to the same functional type in the naming information. This distance threshold of 0.05 is comprehensively estimated based on material processing tolerances and actual equipment installation errors. Then, iterate through each node to compare and confirm all legal edge combinations. For each combination, call the established geometric profile fitting process, limit the fitting target to the form of third-order polynomial approximation, and gradually update the coefficients of the third-order polynomial to meet the condition that the error is less than 0.001. These errors and approximation conditions are obtained by performing linear regression on multiple actual assembly error data. For example, when the cumulative error exceeds 0.005, the polynomial coefficients are reset. Subsequently, perform triangulation calculations on each fitted connection area one by one in the order of functional combinations. A vertical dividing line is inserted at the midpoint of each edge to ensure uniform distribution of the unit side lengths after triangulation. Then, mark the functional category to which each basic unit generated by the triangulation belongs and associate it with the indices of each node, finally forming a boundary feature lattice and a set of functional region triangular elements.

[0061] Based on the above boundary feature lattice and set of functional region triangular elements, read the spatial coordinates of the three vertices of each triangular element one by one, and first calculate the side lengths of the triangle by means of analytic geometry methods. Subsequently, use Heron's formula to obtain the area of the triangle. Heron's formula can be expressed as where a, b, and c are the three side lengths respectively, and p is (a + b + c) / 2. All the letters in these formulas are bound to fixed variable names in the program. Then, input the area and the coordinate vectors of the three sides into the centroid calculation formula, obtain the centroid position of each unit by taking the average value of the coordinates and record it. Then, perform a second-order difference operation on the edges of the unit to obtain the curvature change information. If the curvature change is higher than the empirically measured threshold of 0.1, it indicates that there is a large geometric arc here. This 0.1 is obtained by using the quantile calibration method from the statistics of a large number of typical triangulation samples. When performing quantile calibration, the maximum curvature value is sorted from all samples and the value at the 90th percentile is taken and 0.02 is added to get the result. Finally, store and summarize the geometric dimension parameters, centroid coordinate positions, side length distributions, and curvature change information of all triangular elements in a fixed order to form a set of integrated data for subsequent calls, generating the original geometric feature data.

[0062] The steps to obtain the TES geometric feature set are as follows:

[0063] Extract the node coordinate sequences, boundary contour curvature values, normal direction distributions, boundary node densities, and regional surface areas of all closed boundary regions based on the original geometric feature data. Group the boundary regions into an independent set of candidate regions based on the structural connection relationships to obtain the preliminary data for regional structure division.

[0064] Calculate the morphological distribution dispersion index for each candidate region based on the preliminary data for regional structure division. The calculation formula is:

[0065]

[0066] where D z is the morphological distribution dispersion index, A i is the surface area of the i-th candidate region, L i is the boundary node density of the i-th candidate region, is the average normal vector of the i-th candidate region, is the boundary center vector of the i-th candidate region, θ i is the maximum boundary normal angle of the i-th candidate region, and k is the number of candidate regions;

[0067] Screen the set of candidate regions with different structural attributes according to the morphological distribution dispersion index, and combine the normal vector direction distribution and node density change trend of each candidate region to perform functional region classification and annotation, generating the TES geometric feature set.

[0068] Specifically, based on the original geometric feature data obtained previously, the node coordinate sequence of each closed boundary area is first read, and the corresponding curvature data, normal direction distribution and node density information are extracted point by point from the three-dimensional space according to the recorded coordinate values. Then, these node data are regrouped according to the division of the boundary to which they belong. In the process, a comparison range needs to be established for the number of nodes in each boundary. For example, the number of nodes is compared with the pre-established valid range of 10 to 200. When the number of nodes is within this range, it is regarded as a valid boundary, otherwise the boundary is recorded for further inspection. Subsequently, the contour curvature value of each valid boundary area is analyzed one by one, and each curvature value is matched with the preset curvature range of 0.05 to 0.3. This range is the quantile range obtained after measuring multiple devices of the same type, and the average value of 0.1 is used as the obvious curvature. Detection benchmark, when the curvature is detected to be greater than 0.3, it is determined to be a boundary area that needs to be segmented and identified, and then the node density information of the corresponding boundary is read, and the density value is compared with the standard range of 3 to 50. This is based on the experience of the node distribution of the same scale system in the past. At the same time, the surface area of each area is recorded and the regional identification index is established. Then, based on the structural connection relationship, the boundaries containing higher curvature or a certain number of node densities are combined into an independent candidate area set. For those areas with too low node density or small curvature, they are divided into ordinary area sets. If it is found in the process that the surface area of some areas is less than 0.1 square meters, it is considered that they may correspond to small transition segments when grouping. By marking such areas as micro-component areas and retaining the corresponding indexes, the preliminary regional structure division data is obtained after confirming that there is no duplication or missing data in all independent candidate areas.

[0069] formula: The formula is useful because it takes into account the surface area A i and the boundary node density L i The coupling relationship between them is obtained, and the dot product of the average normal vector and the boundary center vector and the maximum boundary normal angle θ are introduced in the denominator. i , thus more comprehensively describing the spatial geometric attributes and distribution of the region when measuring the morphological distribution characteristics of the candidate region;

[0070] A i The steps for obtaining the parameters are as follows: first, using the coordinate information of the candidate region boundary nodes obtained above, the surface of the candidate region is divided into several tiny triangular units based on the polygon splitting method, and then the area of each triangular unit is calculated using the Heron formula. The surface area of the candidate region is obtained by adding up the areas of all triangles, and this value is defined as A. i, the surface area is usually verified by measuring the contour and thickness characteristics of the candidate area. For the purpose of quantification, the measured area of a certain candidate area in the actual system is 2.5 square meters. At this time, A i = 2.5;

[0071] L i The steps to obtain the L parameter are as follows: First, record the number of all boundary nodes on the candidate area, define the ratio of this number to the perimeter of the area as the preliminary density value, and then correct it according to the corresponding relationship between the number of intermediate nodes and the area to obtain the final density value L i , for the sake of clear illustration, it can be assumed that the perimeter is 5 meters and the total number of nodes is 30. Then the preliminary density value is 30 / 5 = 6. If there are an additional 5 functional interface nodes in the central part of the candidate area and these interface nodes are located inside the area rather than at the boundary, then they need to be added to the overall density calculation with the area of the area as the distribution factor. Assuming that this part of the increment is equivalent to 2 boundary nodes, the total number of nodes is 32, and finally L i = 32 / 5 = 6.4;

[0072] The steps to obtain the parameter are as follows: By extracting the normal vectors of all tiny triangular surface elements in the candidate area, sum the three components respectively according to the same coordinate reference to obtain the average normal vector If the normal vector of each triangular surface is represented by (n x , n y , n z ), after adding all the normal vectors of the triangular surfaces, we can get (15, 20, 25). Then calculate the modulus of this vector. The modulus Then, after normalization, we get which is correspondingly

[0073] The steps to obtain the parameter are as follows: First, calculate the geometric center coordinates of the candidate area. Specifically, calculate the average values of the x, y, and z coordinates of all boundary nodes respectively to obtain the center point of the area Then convert this coordinate point into vector form For example, if the average value of the x coordinates of all boundary nodes of a certain candidate area is 2.2 meters, the average value of the y coordinates is 3.1 meters, and the average value of the z coordinates is -1.5 meters, then

[0074] θ i The steps to obtain the θ parameter are as follows: Calculate the included angle between the normal vector of each triangular surface element in the candidate area and one by one, and select the largest included angle as the maximum boundary normal angle θ i , for example, if the statistical result shows that the normal vectors of most triangular surfaces in this area are The included angle is less than 5 degrees, but the included angle in a certain local area reaches 12 degrees, then take θ i = 12°;

[0075] The steps to obtain the k parameter are as follows: count the total number of regions that can satisfy the complete attribute record among all candidate regions in the current stage. If 40 candidate regions are screened at this time, then k = 40;

[0076] Calculation process:

[0077] When substituting the above parameters, if a certain candidate region A i = 2.5, L i = 6.4, θ i = 12°, then the denominator part is |15.3| + 12 3 = 15.3 + 1728 = 1743.3, and then calculate the numerator Among them So inside the numerator Take the natural logarithm again to get ln(0.00907) ≈ -4.705. If it is assumed that the sum of the values obtained for all 40 candidate regions is 156.35 after adding, then the summation term Divide it by k = 40 to get approximately 3.90875, and finally perform the square root operation to obtain D z ≈ 1.9769;

[0078] This result indicates that if the D z value is large, it means that the morphological distribution of the overall candidate regions is relatively discrete. If the D z value is relatively small, it means that the candidate regions are more consistent in terms of geometric shape and distribution, providing a quantitative basis for screening candidate regions with different structural attributes in the follow-up.

[0079] According to the morphological distribution dispersion index D obtained previously z , after summarizing the index values corresponding to each candidate region, compare and analyze them with the normal vector direction distribution and node density change trend. In the process, first record the average normal vector of each region as a numerical vector and split it into components according to the three-dimensional coordinate axes. If any coordinate component exceeds the value range of plus or minus 10, then mark that the normal vector offset of the region is large and record the node density value at the same index. Then compare these node density values with the effective range of 5 to 100. If it exceeds this range, then divide the region into the category of larger node density, otherwise classify it into the ordinary density category. Then observe the D of each region zWhether the value exceeds 2.0. This 2.0 is obtained by adding 0.2 to the median value from the previous morphological test statistics of similar systems. When it is found that the D z values in some regions are higher than 2.0 and the normal vector offset is also high, they are classified as special functional areas, and the areas with node density greater than the pre-recorded threshold are marked as fluid guiding areas or heat dissipation enhancement areas. Then, further confirmation is carried out according to the specific component usage. After all the allocations are completed, corresponding functional labels are assigned to each area in the record. These label information, together with the D z values and normal vector data, constitute a set of regional property tables, and finally, they are associated with the previously identified geometric positioning information to obtain the TES geometric feature set.

[0080] The steps for obtaining material and mesh configuration parameters are as follows:

[0081] According to the TES geometric feature set, match the predefined material identifiers corresponding to each functional area label, and extract the thermal conductivity, specific heat capacity, density, heat conduction structure number, and surface smoothness level to generate a set of regional material property parameters;

[0082] According to the set of regional material property parameters, calculate the local mesh resolution control value. The calculation formula is:

[0083]

[0084] where, R m is the local mesh resolution control value, k t is the regional thermal conductivity, μ is the surface smoothness level, c p is the specific heat capacity, ρ is the density, γ is the heat conduction structure number, σ s is the maximum surface area of the region, λ s is the boundary gradient change ratio;

[0085] According to the local mesh resolution control value, assign the local mesh resolution control value to the corresponding functional area, and respectively mark the grid accuracy of the fluid boundary layer and the grid shape factor of the phase change region to form the material and mesh configuration parameters.

[0086] Specifically, according to the TES geometric feature set obtained previously, read the corresponding table of each functional area label and predefined material identification, match the functional area label and material identification in a one-to-one mapping manner, and then obtain parameters such as thermal conductivity, specific heat capacity, density, heat conduction structure number, and surface smoothness level for each functional area. During this process, a material list needs to be prepared in advance and the specific number and thermophysical data source of each material should be listed in it. For example, the intervals of thermal conductivity and specific heat capacity are obtained through on-site measurement or professional literature data collection, and these data are recorded in the material identification table after collection. Subsequently, retrieve the material number corresponding to the matched functional area from the table, extract the corresponding thermal conductivity and specific heat capacity values, and synchronously read in the density value and heat conduction structure number according to the established index association method. The heat conduction structure number needs to be allocated one by one according to different material categories. For example, the numbers from 101 to 200 are set within the range of alloy components, and the numbers from 201 to 300 are set within the range of composite materials, and so on. The surface smoothness level can be represented by an integer interval from 1 to 10 increasing gradually according to the processing accuracy, where 1 represents relatively rough and 10 represents very flat. To clearly define the segmentation, the value of each level of surface smoothness should correspond to the actual measured surface roughness range. For example, the Ra value between 1.6 μm and 3.2 μm is classified as level 4, and the Ra value between 0.8 μm and 1.6 μm is classified as level 5. Thus, when reading the material corresponding to the functional area, it is allocated to the corresponding smoothness level according to the roughness actually measured of the component. Once the above material property elements are extracted, the material data entries of each functional area will be summarized during the processing to ensure that the thermal conductivity, specific heat capacity, density, heat conduction structure number, and surface smoothness level of each area are recorded in a set of regional material property parameters. Finally, summarize the parameters of all functional areas to form the overall set of regional material property parameters.

[0087] Formula: The advantage of the formula is that it combines multiple factors such as material thermal conductivity, surface smoothness level, specific heat capacity, density, heat conduction structure number, and regional area and gradient characteristics, achieving a comprehensive consideration of material properties and geometric elements when calculating the local grid resolution control value, thus more precisely reflecting the changes in the heat conduction process during the grid division stage;

[0088] k t The steps to obtain the parameter are as follows: First, determine the heat conduction performance of the material in the functional area. Measure the heat conduction performance of the material by arranging heat flux monitoring points in the actual use environment. Each monitoring point records the heat value conducted per unit time and calculates the local heat conductance value in combination with the temperature gradient. Then, perform a weighted average of these local heat conductance values to obtain the overall heat conductivity. To make the weighting process quantifiable, a monitoring point area coefficient α can be defined. i, when the coverage area of the monitoring points is large, α i Take a relatively high value, and vice versa take a relatively low value. The thermal conductivity values of each monitoring point are combined through the following formula: where represents the thermal conductivity measured at the i-th monitoring point, m is the number of monitoring points. In actual projects, when measuring areas with different shapes or thicknesses, α i is usually determined by the ratio of the coverage area of the monitoring points to the total area of the region. For example, if 5 monitoring points are arranged in a certain functional area, and the coverage area of each point is 0.09 square meters, 0.08 square meters, 0.06 square meters, 0.1 square meters, and 0.07 square meters respectively, and the total area of the region is 0.4 square meters, then the α i of these 5 monitoring points can be set to 0.225, 0.2, 0.15, 0.25, and 0.175 respectively. Multiply the thermal conductivity value measured at each monitoring point by the corresponding α i and then sum them up, and then divide by the sum of α i to obtain k t , for example, finally obtain k t = 35.2 W / (m·K);

[0089] The steps to obtain the μ parameter are as follows: Measure the surface roughness of the functional area in actuality, and give a quantification level from 1 to 10 according to the surface flatness. The surface roughness interval corresponding to each level is obtained through measurement by a metal surface profiler. The corresponding relationship can be presented using a mapping table. For example, level 1 corresponds to Ra above 25 μm, and level 10 corresponds to Ra below 0.1 μm. When measuring, select at least four measurement points on the surface of the functional area, detect the roughness values one by one, take the average and compare with this mapping table to determine the integer value of μ. For example, if the measured average value of the surface roughness is 0.8 μm, corresponding to level 5, then μ = 5;

[0090] c p The steps to obtain the parameter are as follows: When measuring the specific heat capacity of the material, it is necessary to use a calorimeter for detection. First, heat a sample with a specified mass (such as 0.5 kg) under a controlled environment, and monitor the relationship between the heat input and the temperature rise in real time. Use to calculate the specific heat capacity, where Q is the heat input, m is the sample mass, and ΔT is the temperature rise value. In order to obtain more stable and accurate data, five heating experiments can be carried out within the temperature range of 30 °C to 80 °C, and the arithmetic mean of the five results is taken and recorded as the final value of c p , for example, the detection results of a certain region are summarized to obtain c p ≈ 900

[0091] J / (kg·K);

[0092] The steps to obtain the ρ parameter are as follows: Calculate the density by weighing and volume measurement on-site. When the regional materials can be directly sampled, put the materials with a known volume on an electronic scale and read its mass value. According to Obtain the density by taking the ratio of the mass m to the volume V, and then complete the calculation in combination with the mass obtained from the electronic scale. For example, a material sample with a volume of about 2.5 liters is weighed to obtain a mass of 2.2 kilograms, then ρ = 2.2 / 0.0025 = 880 kg / m 3 ;

[0093] The steps to obtain the γ parameter are as follows: When determining the heat conduction structure number, it is necessary to first analyze the heat conduction stratification of the regional materials. For the regions with multi-layer composite structures, it is necessary to calculate their overall heat conduction resistance and assign a number corresponding to the resistance range. The number range can be defined as between 100 and 999, and the specific value is determined by the heat conduction resistance interval after comprehensive detection. For example, for single-layer metal components, the numbers are taken from 100 to 199, for multi-layer composite metal materials, the numbers are taken from 200 to 299, for ceramic or non-metal materials, the numbers are taken from 300 to 399, etc. Compare the detection results to determine the corresponding number and record it in the system. For example, determine γ = 215;

[0094] σ s The steps to obtain the parameter are as follows: After scanning the geometric outer surface of the same functional area, count the largest continuous plane or curved surface area among them, and use this area as σ s to represent the largest surface area of the region. The detection of this area can be completed by three-dimensional contour measurement and splitting its projection into multiple triangular units for cumulative summation. If the scanning result shows that the largest plane in this region is about 1.2 square meters, record σ s = 1.2;

[0095] λ s The steps to obtain the parameter are as follows: Measure the temperature or heat flux density gradient at the boundary of the functional area, and take the ratio of the difference between the maximum and minimum gradients in this area to the average gradient. Define this ratio as the boundary gradient change ratio λ s , for example, the measured maximum temperature gradient in the region is 25 K / m, the minimum temperature gradient is 15 K / m, and the average gradient is about 20 K / m, then

[0096] Calculation process:

[0097] After all the parameters are obtained and substituted into the formula, for example, take k t = 35.2, μ = 5, c p = 900, ρ = 880, γ = 215, σ s = 1.2, λ s = 0.5, then first calculate log 10 (k t·μ) part:

[0098] k t ·μ = 35.2×5 = 176;

[0099] log 10 (176) ≈ 2.2455;

[0100] Then calculate part:

[0101] c p ·ρ = 900×880 = 792000;

[0102]

[0103] Add the two together:

[0104] 2.2455 + 890 ≈ 892.2455;

[0105] And take the square:

[0106] (892.2455) 2 ≈ 796106.67;

[0107] Then calculate the denominator

[0108] γ·σ s = 215×1.2 = 258;

[0109]

[0110] ln(1.25) ≈ 0.2231;

[0111] So the denominator part is approximately 258 + 0.2231 = 258.2231, and finally find Then take the square root of it to get:

[0112]

[0113] This result shows that R m The higher the value of, the greater the resolution requirement for the local grid, and a finer grid division can be adopted for the corresponding area. If R m is lower, a relatively coarser grid can be used.

[0114] According to the local grid resolution control value R obtained previously m , establish an index according to the material property parameter set of each functional area, and use R mMatch and assign each item according to the functional area label to the corresponding areas. During the process, it is necessary to record the identifier of each functional area in the management table and compare it with the corresponding heat conduction structure number and thermal property data. Subsequently, a relatively higher grid accuracy is set for the areas belonging to the fluid medium transmission part. Referring to an accuracy level range of 5 to 10, corresponding to the detection data at this stage, if the R of a certain functional area m exceeds 50, it is classified as a high-precision area and its accuracy level is raised to 8 or 9. If it is lower than 30, it is classified as an ordinary precision area and its accuracy level is maintained between 5 and 6. These threshold ranges are obtained based on the statistical records of on-site heat transfer analysis. The grid sensitivity values measured during the operation of multiple batches of the same type of equipment are sorted and the upper and lower percentiles are extracted from them. Then, when the functional area of the phase change part is detected, it is necessary to further retrieve the phase change material category. If the R m exceeds 40, polyhedral meshes are used and the number of mesh elements is appropriately increased to label a larger mesh shape factor. If the R m is lower than 40, relatively simplified tetrahedral meshes are used. Finally, after all areas are labeled, the final material and mesh configuration parameters are obtained by summarization.

[0115] The steps to obtain the FEA preprocessing configuration file are as follows:

[0116] According to the material and mesh configuration parameters and the channel characteristics identified from the TES geometric feature set, locate the geometric boundaries corresponding to the heat flow inlet and heat flow outlet in each functional area. Based on the boundary node normal vector direction, node distribution density, and geometric area type, identify and map the labels of the inlet and outlet units to generate the fluid channel boundary node attribute information;

[0117] According to the fluid channel boundary node attribute information, extract the surface area, normal vector distribution, and mesh element type corresponding to each inlet and outlet unit. Combining the material thermal properties and mesh element distribution patterns in the material and mesh configuration parameters, set the boundary condition types for the inlets and outlets in each functional area, including the dominant mode and action dimension of the heat flux boundary or velocity boundary, to generate the boundary condition setting information;

[0118] According to the boundary condition setting information, call the mesh division details and material thermodynamics properties of each unit area in the material and mesh configuration parameters to construct the node, element, boundary item structure, and constraint type syntax content of the finite element preprocessing input fields, and generate the FEA preprocessing configuration file.

[0119] Specifically, according to the channel features identified by the material, mesh configuration parameters, and TES geometric feature set obtained previously, first find the geometric boundary information corresponding to the heat flow inlet and heat flow outlet in each functional area, and read the relevant boundary node coordinates and functional labels. Check the normal vector orientation of these nodes one by one and the distribution density of the nodes in space. If the normal vectors of the boundary nodes are concentrated within a certain angular range and the distance between the nodes is less than the 8-mm interval threshold obtained through actual equipment measurement and comprehensive analysis, they are grouped into candidate groups for the same inlet or outlet. This 8-mm interval threshold is a value formed by taking the average of the boundary node distributions measured in multiple experiments and adding a 1-mm safety margin. When the geometric region type of the candidate group is identified as a fluid channel, it is necessary to further check whether the normal vector of this group is consistent with the previously recorded channel direction data. For example, when the channel direction data indicates an angle in the range of zero to thirty degrees and the average normal vector angle of the nodes in this group is twenty-five degrees, it can be determined that the matching degree is relatively high. Subsequently, record the specific labels of each inlet group and outlet group during the summary process, and check whether there are multiple candidate groups overlapping. If it is detected that the sum of the number of nodes in multiple groups differs from the actual system calibration number by more than 5 nodes, it is necessary to conduct a review based on the boundary location and shape survey data to confirm whether there is a phenomenon of missed selection or duplicate selection. After all inlet units and outlet units are identified without repetition, perform label mapping on the node indices corresponding to each unit and write them into a structured record, and finally generate the fluid channel boundary node attribute information.

[0120] According to the fluid channel boundary node attribute information obtained previously, read the surface area and normal vector distribution corresponding to its inlet unit and outlet unit one by one. Based on the thermal property parameters such as thermal conductivity, specific heat capacity, and density retrieved from the material and mesh configuration parameters before, compare the possible heat transfer situations at the inlet and outlet with the distribution pattern of the mesh elements. During this process, compare the surface area with the previously collected measured values. If the difference between the two exceeds 3% (this proportional threshold determined based on the measurement data of multiple test devices), a node repeatability check will be performed to confirm whether there are some unmarked nodes in the spatial topology information. Then, specifically set the boundary condition types at the inlets and outlets in the functional area. By comparing the previously recorded heat flux and fluid velocity data, identify the dominant boundary attribute of a region as a heat flux boundary or a velocity boundary, and mark it separately in the three-dimensional coordinate directions according to the difference in the acting dimension. If the fluid velocity value in this region is detected to be within the range of 1 meter per second to 3 meters per second, it is inclined to consider the velocity boundary as the dominant one. If the heat flux is in the range of 500 watts to 2000 watts, it is inclined to consider the heat flux boundary as the dominant one. All the threshold ranges are quantitative standards formed by combining the operating indicators provided by the equipment manufacturer and the statistics of several groups of operating data. Finally, after completing the description of the dominant mode and acting dimension for all inlets and outlets, summarize and generate the boundary condition setting information.

[0121] According to the boundary condition setting information obtained previously, read the mesh division details and material thermodynamic properties of each functional area in the material and mesh configuration parameters. Cross-match the unit division scheme and the regional thermal property values item by item according to the index order of the node numbers. Then, add the corresponding control syntax content for each boundary node that needs to establish a constraint relationship. For example, write the parameter values related to the heat flux or the coordinate directions and calculation dimensions related to the velocity boundary under the corresponding index. Among them, if the thermal conduction structure number of some functional areas exceeds 200 (the starting value of the multi-layer or composite structure divided in the previous detection), more boundary instructions need to be added to cover the heat transfer switching process between the inner layer and the outer layer, and configure more triangular or tetrahedral units according to the corresponding mesh shape factor. Next, arrange the relevant information in a unified format according to the arrangement order of the nodes and units. If it is detected that there are still individual boundaries that cannot find the matching thermodynamic properties after the allocation is completed, the previously measured data records need to be retrieved and compared again to confirm whether some thermal conduction numbers or mesh categories are missed. Once the settings of each unit area are all completed, integrate the node, unit, boundary item structure, and constraint type syntax into the same processing set, and finally generate the FEA preprocessing configuration file.

[0122] The steps to obtain the original FEA result data are as follows:

[0123] Read the finite element analysis result file of the FEA preprocessing configuration file target, retrieve the result index tag corresponding to the analysis task, locate the temperature distribution result path and heat flux density vector result path in the output file through the mapping field, and generate a set of finite element result target paths;

[0124] According to the set of finite element result target paths, sequentially load the temperature field and heat flux density field information within each time step or iteration step, extract all temperature values and heat flux vector values according to the node number, and uniformly organize the extraction results into a structured matrix format to generate a temperature distribution matrix and a heat flux density vector matrix;

[0125] According to the temperature distribution matrix and the heat flux density vector matrix, compare the mesh structure and node numbering rules set in the FEA preprocessing configuration file, perform matrix node order verification and boundary matching detection, process invalid numerical units and perform node structure binding to generate the original FEA result data.

[0126] Specifically, when reading the finite element analysis result file of the FEA preprocessing configuration file target, first retrieve the index tag of the analysis task and match it with the corresponding task identifier stored in the system. After successful matching, reference the index tag to the internal record, and then search for the field path associated with the index tag in the pre-compiled field comparison table and distinguish the temperature distribution data and heat flux density vector data according to the identifier category. If the length of the searched field path is insufficient or the path format does not conform to the retrieval rule stratified by " / ", mark the path as an abnormal state and record it. When manual review is required, compare each record one by one to confirm whether there are spelling mistakes in the path name or incomplete previous index pointers. If the retrieval result is normal, perform a structural analysis on the path content. When analyzing, first verify whether each path contains the corresponding file paragraph description and compare whether the number of sub-fields is within the empirical range from 1 to 10. If the number of sub-fields exceeds this range, check whether there are blank fields or duplicate fields. Once it is confirmed that both the path format and the number of sub-fields are correct, include it in the target path set. Subsequently, a brief file type identification is also required for each target path to ensure that the storage formats of the temperature distribution and heat flux vector data match. If the file type is inconsistent with the expected format, separate the path and uniformly mark it in the abnormal record, waiting for the inspection personnel to judge whether the path points to the wrong or the source data is incomplete by referring to the previously established file list. When all files pass the verification and the number of the target path set is the same as the expected number in the index description, summarize these paths and compile them into a set of finite element result target paths.

[0127] When loading the temperature field and heat flux density field information according to the order of time steps or iteration steps based on the set of finite element result target paths obtained previously, it is necessary to first read the file structure for each time step and compare it with the existing time step configuration before reading. If it is found that the time step label is greater than the upper limit threshold of 50 defined previously, it is prompted that there may be an inconsistency in the task configuration and the specific difference information is recorded. After passing the check, the node numbers and the corresponding temperature values are extracted one by one. Each node number is combined with the temperature field at the current time step, and then the heat flux density vector value is extracted and the vector components are split into three components and respectively matched with the node numbers. This can make all the temperature values and heat flux vector components have a unified index structure. Next, all the nodes together with their temperature values and vector components are put into a matrix structure for unified management. The temperature distribution is defined in the form of a two-dimensional matrix, where each row represents a node and each column represents the data of a time step or iteration step. The heat flux density vector components are also integrated into another set of matrices to form a vector matrix set corresponding to the temperature distribution. During the process, if it is detected that some node numbers are not within the set range of 1 to 20000, these nodes are marked as out of bounds and the corresponding rows are removed. If the temperature value is lower than -50°C or higher than 2000°C, which are two thresholds determined based on the device tolerance limit respectively, it is considered that the temperature fails and it is necessary to supplement and check whether there is a measurement error for this node in the source record. After all the temperature values and vector components are loaded, the complete temperature distribution matrix and heat flux density vector matrix are obtained.

[0128] According to the temperature distribution matrix and heat flux density vector matrix obtained previously, when comparing with the mesh structure and node numbering rules set in the FEA preprocessing configuration file, it is first necessary to query the element topology relationship corresponding to each node in the mesh structure information table, and check whether the arrangement of node numbers corresponds one-to-one with the internal row numbers of the matrix. If the numbering deviation exceeds 3, which is a standard value obtained based on empirical data collection, then it is necessary to count whether the total number of these numbered deviation nodes exceeds 2% of all nodes. When it exceeds 2%, these batches of nodes are uniformly marked as out-of-order and their initial allocation positions are traced. Then, scan for possible null or missing values in the matrix, and perform replacement or deletion processing on these invalid value elements. Replacement can be done by comparing the average temperature of neighboring nodes for patching, or directly removing the element from the data structure if it exceeds the equipment applicable range. Subsequently, perform boundary matching detection to observe whether the entries of all nodes marked as boundary nodes are consistent with the configuration information in terms of temperature and heat flux distribution. If it is found that a boundary node marked as a fluid inlet region does not detect a positive flow component in the heat flux vector information, then perform a re-comparison in combination with the functional region label to check whether there was a regional classification confusion in the early stage. Finally, after all nodes are confirmed to match the preprocessing configuration, perform the final binding of node indices and matrix rows and columns, and incorporate the revised data into the structured result set to generate the original FEA result data.

[0129] The steps for obtaining the surface heat flux density gradient field are as follows:

[0130] According to the original FEA result data, extract the surface mesh node indices and geometric coordinate information of the original CAD model, and make the node coordinates correspond one-to-one with the temperature distribution data and heat flux density vector data through the indices to generate a combined attribute set of the thermophysical and geometric structures of the mesh nodes;

[0131] According to the combined attribute set of the thermophysical and geometric structures of the mesh nodes, calculate the heat flux density gradient intensity value of each node. The calculation formula is:

[0132]

[0133] Where, G x is the heat flux density gradient intensity value, is the heat flux density vector of the i-th node, is the temperature gradient vector of the i-th node, is the node position vector, is the node normal vector, AG i is the mesh area where the node is located, φ i is the angle between the node heat flux vector and the normal vector;

[0134] According to the intensity value of the heat flux density gradient, value interpolation and reconstruction are carried out for each grid node in the order of the number, and the contour distribution map of the node heat flux gradient is constructed by combining the topological connection structure of the geometric patches. The field data is filled in the space position mapped according to the grid node coordinates on the surface of the original CAD model to generate the surface heat flux density gradient field.

[0135] Specifically, based on the original FEA result data, first read the grid node index and geometric coordinate information on the surface of the original CAD model, and proofread the coordinate values of each node from multiple angles to confirm that it is consistent with the previously generated grid structure in the three-dimensional space. If it is found that the difference between the node coordinates and the reference node index is greater than the 0.2 mm threshold summarized from previous tests, the node will be put into the recheck queue and its coordinate difference will be recorded. Then continue to retrieve the corresponding temperature distribution data and heat flux density vector data in the node index list, and compare the temperature value and the three components of the heat flux density of each node one by one to see if the temperature value is between 0 °C and 90 °C. This interval is obtained according to the conventional operating conditions during the operation of the equipment. If the temperature value is lower than 0 °C or exceeds 90 °C, it will be marked as needing to be re-verified. The interval detection is also carried out for the heat flux density components. For example, the heat flux density is set to a reasonable interval of 0 W / m2 to 3000 W / m2 with reference to the on-site monitoring conditions. For nodes exceeding this interval, the values of their adjacent nodes will be compared again to check if there are measurement or recording errors. Next, bind the temperature and heat flux density information to the corresponding node index to ensure that each node has both coordinate data and thermophysical data. Then generate a combined attribute entry of thermophysics and geometric structure through the node index, and mark both the temperature and heat flux components in the same structure. This can directly call the corresponding geometric coordinates and thermophysical properties according to the same node number in the subsequent processing. During this period, if it is found that the heat flux component of an individual node is missing, it is necessary to go back to the extraction process in the previous step to check if there is any missing data, and also check if there is a situation where the node number is written wrong. When all nodes have completed data merging and there are no more out-of-bounds or missing records, the generated node index, coordinates, and the three components of temperature and heat flux density will be stored uniformly to obtain the combined attribute set of grid node thermophysics and geometric structure.

[0136] Formula: The benefit of the formula is that by introducing multi-dimensional elements such as heat flux density vector, temperature gradient vector, node position vector, node normal vector, grid area where the node is located, and included angle, etc., it realizes the coupled consideration of local heat flux distribution and geometric shape, and can quantitatively reflect the interaction between the heat flux density and temperature gradient at the node in the spatial direction;

[0137] The steps for obtaining the parameters are as follows: In the original FEA result data extracted previously, the vector value of the heat flux density, including three components of x, y, and z, is recorded for each node. When obtaining, it is necessary to first perform unit conversion on the measured heat flux data. If the original data is measured in W / m 2 measurement, maintain this unit. If there are different dimensions such as W / cm 2 or kW / m 2 etc., uniformly convert them to W / m 2 and then store. For the sake of accuracy, multiple nodes can be detected through multi-point measurement. Integrate the heat flux measurement values at different positions and take the weighted average. The weight coefficient is determined by the proportion of the actual covered surface area of the node. The heat flux vector of the i-th node can be summarized using the following formula: where is the measured value of the heat flux vector of node i in the j-th sampling, and β j is the node coverage ratio coefficient corresponding to this measurement, and m is the number of samplings. For example, if the heat flux components measured at node i in four samplings are (200, 150, 100), (205, 152, 98), (198, 148, 102),

[0138] (202, 151, 99) W / m 2 , and the corresponding coverage ratio coefficients are 0.25, 0.3, 0.2, and 0.25 respectively, then the sum of β j is 1. Finally, multiply the four groups of data by their respective β j and add them together to obtain a single vector;

[0139] The steps for obtaining the parameters are as follows: Perform a differential solution on the temperature gradient around the same node position. First, collect the temperature values at adjacent grid nodes, and then use the spatial difference formula to calculate the temperature change rate. Specifically, finite differences can be performed in the x, y, and z directions respectively to form During the calculation, divide the average temperature difference between adjacent nodes by the coordinate difference to obtain the gradient component. To exclude abnormal values, a radius d can be set around this node. If the component differs from the adjacent area by more than the set threshold of 10 K / m, then review the temperature data of this node again. Finally, retain the qualified gradient components. For example, at a certain node i, the temperatures of several surrounding nodes are 45 °C, 46 °C, and 47 °C respectively, and the coordinate difference from node i is about 0.5 meters. Through differentiation,

[0140] The steps to obtain the parameters are as follows: Pack the x, y, and z coordinates of node i in the CAD model coordinate system into a vector. If the node coordinates are (1.0, 1.5, 2.0) m, then r i = (1.0, 1.5, 2.0). When obtaining in the scene, it can be directly read from the grid node and coordinate comparison table generated previously, and then verified. If there is a deviation between the coordinate value and the previous calibration, re-compare the node index and the grid structure until they are confirmed to be consistent;

[0141] The steps to obtain the parameters are as follows: In the previous stage, normal vector information has been stored for each node. The normal vector is usually obtained by calculating the average normal of the grid patch where the node is located. For each triangular or tetrahedral element, there will be a normal vector. After collecting and accumulating them and normalizing, the average normal vector of the node can be obtained. For the surface of metal or homogeneous materials, the average of all element normals can be calculated at the geometric center. If the normal vector of the grid surface where the node is located is (0, 1, 0) or (0.1, 0.95, 0.05), etc., it is weighted according to the area ratio and then normalized to obtain the final value. If it is found that the length of the normal vector after weighted normalization differs from 1 by more than 0.001, re-traverse the element normal vectors and correct them, and finally output in the standardized form of (0, 1, 0);

[0142] AG i The steps to obtain the parameters are as follows: In the grid patch where the corresponding node is located, count the area value of the grid patch and allocate it according to the proportion occupied by the node. If it is a triangular meshed patch, the area can be calculated using Heron's formula, and then the proportion is divided according to the angle of the node. If it is a tetrahedral element, it can be divided into several small triangles through the known coordinates and then summarized. After determining the effective coverage area of the node, let it be the value of AG i If the patch area is found to be 0.4 m 2 from the recorded grid structure, and the node occupies approximately 25% of the angle division weight in this patch, the value of AG i = 0.4 × 0.25 = 0.1 m 2 can be obtained;

[0143] φ i The steps to obtain the parameters are as follows: The cosine value of the angle between the heat flux vector of this node and the normal vector needs to be obtained first through the dot product of vectors and then converted into an angle representation. For example, if where should be 1. If is about 300 W / m 2 and its dot product with is 260, then φ i≈30°;

[0144] Calculation process:

[0145] Taking a certain node i as an example, take φ i = 30° and perform cross operation and dot product operation:

[0146] First It can be calculated:

[0147]

[0148] Then It can be calculated as:

[0149]

[0150] Then the dot product ((300, -600, 150)·(-2, 0, 1)) = (300 × -2) + (-600 × 0) + (150 × 1)

[0151] = -600 + 0 + 150 = -450, take the absolute value to get 450, further |…| 2 = 450 2 = 202500;

[0152] Inside Its square is approximately 0.001539;

[0153] cos 2 (φ i ) + 1 where φ i = 30°, After squaring, 0.866 2 ≈ 0.75, add 1 to get 1.75;

[0154] Substitute these values into the numerator of the formula:

[0155] 202500 × 0.001539 ≈ 311.8375;

[0156] The denominator is 1.75, so:

[0157]

[0158] Then take the cube root of it:

[0159]

[0160] This result shows that when G x = 5.63, the coupled distribution of heat flux density and temperature gradient at the node is at this order of magnitude. If G xA larger value indicates a more significant local clustering phenomenon between heat flux and temperature gradient, while a smaller value indicates a weaker coupling;

[0161] When interpolating and reconstructing each grid node in the order of its number according to the heat flux density gradient intensity values obtained previously, it is necessary to first map these gradient values one by one with the node indices and mark their coordinate positions in three-dimensional space in the same list. Subsequently, read the face topological connection information of the grid, and based on the node relationships on the faces, gradually smooth the gradient values of each node. During this process, if it is found that the difference in gradient values between adjacent nodes exceeds 10, which is an empirical threshold obtained from a large number of previous model statistics, calculate their average value and compare it with the previous node's gradient distribution record to confirm whether there is an abnormal jump. If the proportion of the number of nodes in the jump region to the total number of nodes in the current face exceeds 20%, it is necessary to further check whether this face is marked as a large curvature or high heat flux difference region in the original FEA result data. If it is in a high heat flux difference region, retain this part of the larger gradient difference; otherwise, correct it to the intermediate level of the gradient values of the surrounding nodes during the interpolation process. Then, perform linear or quadratic interpolation transitions on all interpolated node gradient values along the topological connections between the faces, so that the gradient distribution can form a continuous change trend on the grid surface. Finally, transfer the gradient values of each node to the visualization mapping program in the order of their numbers and complete the filling in the three-dimensional coordinates, overlay the gradient distribution contour lines on the surface grid of the original CAD model, and display different gradient intensity ranges by differentiating colors, vividly presenting the gradient level corresponding to each node on the model surface to generate the surface heat flux density gradient field.

[0162] The steps to obtain the evaluation results of the heat recovery rate of thermal energy storage are as follows:

[0163] According to the TES geometric feature set, identify the energy input surface, energy output surface, and the interface of the energy storage material domain, mark the boundary edges through node indices and the grid face topological structure and bind the facing directions. At the same time, call the heat flux density vector and node temperature field data in the original FEA result data to form the target surface thermophysical mapping structure;

[0164] According to the target surface thermophysical mapping structure, calculate the evaluation value of the heat recovery rate of thermal energy storage. The calculation formula is:

[0165]

[0166] Among them, Ψ is the evaluation value of the heat recovery rate of thermal energy storage, is the surface heat flux density vector, is the unit vector of the face normal direction, τ s is the node temperature change rate, ξ s is the average temperature gradient of the energy storage material domain interface, S inis the energy input surface, S out is the energy output surface, and t0 and t1 are the start and end boundaries of the preset time window;

[0167] According to the evaluation value of the heat recovery rate of thermal energy storage, combined with the distribution map of the node heat flux directions in the surface heat flux density gradient field, the gradient direction and the heat recovery rate correlation layer are superimposed on the grid node coordinate system of the 3D CAD model through color mapping to generate the evaluation result of the heat recovery rate of thermal energy storage.

[0168] Specifically, after identifying the energy input surface, the energy output surface, and the interface of the energy storage material domain according to the previously obtained TES geometric feature set, it is necessary to locate the coordinates of these surfaces in the grid structure of the 3D model, and sequentially read the numbers of the nodes on each surface and the associated boundary edge data. During the process, the topological relationship between the node numbers and the grid patches will be checked. For example, for each patch, it is confirmed one by one whether the node coordinates coincide with the specified surface edge. The coincident points found will be marked as the input surface, the output surface, or the interface of the energy storage material domain in the number list, and a direction identifier will be assigned to these nodes to record the orientation information of the patch where the node is located. If it is detected that the orientation information is different from the pre-recorded functional attribute, then the difference between the normal vector of the node at this place and the standard surface direction is compared. If the difference in the angle exceeds 15 degrees, which is the empirical threshold obtained during the installation and commissioning stage, the node will be marked as a doubtful record and the items that need to be rechecked will be listed in the record. If the difference in the angle is within 15 degrees, it is considered that it can correspond to the functional attribute. After completing the binding of the nodes and the boundaries, the heat flux density vector and the temperature field data corresponding to the node numbers in the original FEA result data are retrieved, and these heat flux density components and temperature values are mapped to the corresponding node coordinates. If it is detected that some nodes are not within the range of this heat flux density data, the previously generated original FEA result data will be returned for comparison again to check whether there is an error in the number or position input. After all checks are correct, the energy input surface, the energy output surface, and the interface of the energy storage material domain are arranged in the record in sequence, and the node heat flux density and the temperature field values are written into this record together to form the target surface thermophysical mapping structure in tabular form.

[0169] Formula: The benefit of the formula is to compare the heat flux contributions of the energy input end face and the output end face within the same time window, by introducing to consider the influence of the temperature change rate, and further incorporate the average temperature gradient ξ s of the interface of the energy storage material domain into the denominator term, making the evaluation process more comprehensive when measuring the balance relationship between energy storage and heat recovery;

[0170] The steps for obtaining the parameters are as follows: represents the surface heat flux density vector, usually in W / m2 It is expressed in units to represent its magnitude and direction. In order to obtain readings need to be taken at multiple monitoring points distributed on the energy input surface and the output surface, record the heat flux components in the x, y, and z directions, and integrate these monitoring values into discrete surface element nodes through area weighting or density weighting. When the monitoring density is high, sensors can be arranged at the center of each or adjacent grid cells to obtain the corresponding vectors, and then use the following weighted summary formula: where α j is the surface coverage ratio coefficient, obtained through the grid area or the sensor coverage range, m is the number of acquisitions, is the heat flux density vector of the j-th measurement point. Finally, the summary result is bound to each surface element or surface node. For example, when five monitoring points are arranged on the output surface and the average heat flux values measured at each point are in the range of 100 to 200 W / m 2 After weighting and integrating to the cell center, we get

[0171] The steps to obtain the parameters are as follows: is the unit vector in the surface normal direction, which needs to be obtained by normalizing the normal vector of each patch when constructing the grid. If the surface normal vector was originally calculated by cross-multiplying the coordinates of surrounding nodes, then its magnitude needs to be calculated and then divided by the magnitude to obtain the unit vector. When multiple adjacent grid surfaces share a certain boundary, to avoid direction conflicts, the normal direction can be uniformly smoothed during the geometric inspection stage to keep the direction consistent, and then recorded in the surface element information table. For example, if the original value of the normal vector of a grid surface is (0, 3, 4) and the magnitude is about 5, after normalization, it becomes When further combining with surface functions, it will be judged according to the orientation of the input surface or the output surface whether it coincides with the main energy flow direction. If not, it will be corrected to

[0172] τ s The steps to obtain the parameter are: τ s represents the node temperature change rate. When obtaining it, the node temperatures in the energy input and output regions can be collected at different time points, and the difference form is used to divide the temperature change by the time interval to obtain the change rate of each node. To prevent short-term perturbations in the temperature values of individual nodes, the results of multiple samplings need to be smoothed and averaged. For example, a certain node is temperature-read ten times within a period of time, and the sampling interval is 1 second each time. If the temperature rises from 45 °C to 55 °C, then

[0173] ξ s The steps to obtain the parameter are: ξs is the average temperature gradient at the interface in the energy storage material domain. To obtain it, temperature measurements need to be carried out on the boundary cross-section of the energy storage material domain. For example, a measurement point is set at a certain distance (such as 5 mm) in the material thickness direction, and the temperature value is recorded. Then, the overall gradient is obtained through piecewise difference or linear regression methods. For example, for a certain composite energy storage material, the temperature of 6 layers of points is measured from the surface to the inside. The outermost layer is 35 °C, the innermost layer is 25 °C, and the thickness is 0.05 m. Then the average temperature gradient

[0174]

[0175] S in The steps to obtain the parameter are: S in is the energy input surface. During the early identification, all patches that may input heat flux need to be screened, and its range is determined during the engineering design stage. For example, several adjacent patches where the high-temperature fluid inlet is located are collectively called the input surface. The coordinate distribution or topological relationship of these patches can be retrieved from the TES geometric feature set, and then compared with the actual measured data of the sensor to confirm that there is a significant heat inflow. In this way, S in is delimited. If it is detected that the heat flux value of some patches is too low, lower than 10 W / m 2 this lower limit threshold obtained from on-site measurement, it is considered that they do not belong to the typical input surface, and they can be removed to retain the main input surface area. Finally, the confirmed surface set is defined as S in ;

[0176] S out The steps to obtain the parameter are: S out is the energy output surface, opposite to S in Generally, it points to the low-temperature area or the heat dissipation part. When obtaining it, it is also necessary to retrieve those patches in the grid that show heat discharge or a temperature gradient from high to low during monitoring. If the temperature and heat flux monitoring show that the average heat flux direction in this area is towards the external environment, and the value is greater than 10 W / m 2 , then it is considered that this patch belongs to the output part. After comparing the connection topological relationship with the actual sensor distribution range, the boundary node set of S out is finally determined. If the heat flux direction of some patches does not match the overall heat dissipation direction, they are excluded from the output surface;

[0177] The steps to obtain the t0 and t1 parameters are: These two parameters indicate the start and end boundaries of the preset time window. For integration, it is necessary to define a collection period in the engineering plan and record its start and end times in the system. When the system intends to observe the energy storage status within ten minutes, then t0 = 0 s and t1 = 600 s;

[0178] Calculation process:

[0179] In a certain measurement, [t0, t1] = [0, 600] s, S in The area is approximately 1 m 2 , S out The area is approximately 1.2 m 2 , and a number of nodes are arranged on each of these two surfaces to detect the heat flow direction and the temperature change rate. It is assumed that the average projected value obtained by summarizing and weighting is . Then the integral term on the output surface can be approximately expressed as:

[0180]

[0181] If S out ≈1.2 m 2 , this term is approximately:

[0182] 200×0.693×1.2×600≈200×0.693×720=200×499.0=99800;

[0183] The integral term on the input surface is:

[0184]

[0185] If the average temperature gradient ξ s of the energy storage material domain interface is 200 K / m, Substituting these two terms into the numerator and denominator gives:

[0186]

[0187] 0.01229 1 / 3 ≈0.228;

[0188] This result shows that in this time period, after comparing the energy output surface with the input surface, Ψ≈0.228. When the Ψ value is relatively close to 1, it means that the recovery rate is relatively balanced with the input energy. If the Ψ value is much less than 1, it indicates that the recovery is relatively low, and the temperature gradient and transmission efficiency of each part of the system need to be checked subsequently.

[0189] According to the evaluation value of the heat recovery rate of thermal energy storage, first obtain the previously known node numbers and the heat flow direction distribution map of each node in the surface heat flux density gradient field, and correlate the two. When operating, it is necessary to first integrate a list of heat recovery rate values according to the node index. Each record corresponds to the coordinate position of a node and the evaluation value calculated in the early stage. Subsequently, read the established color mapping rules. This rule is usually set according to empirical intervals. For example, the range of heat recovery rate from 0 to 0.3 is designated as the cyan-green color system, 0.3 to 0.6 as the yellow color system, and 0.6 to 1.0 as the orange to red color system. RGB values or HSV values must be specified in the mapping table for each color interval to ensure no repetition. And when the value exceeds 1.0, a purple or magenta mark is additionally set at the end of the color table. The division of these intervals is set after percentile statistics by combining the measurement results of multiple energy storage devices. When counting, the heat recovery rate range of the same type of device is divided into 10 to 20 levels, and the more obvious separation points are selected as the interval boundaries, and the corresponding colors are mapped to these boundaries to form an intuitive difference in the final visualization result. Then, retrieve the flow vector information displayed by each node in the heat flux density gradient field, link its angle with the heat recovery rate value, and hook the color with the vector direction according to the predefined color band order. If the angle between the vector direction and the normal vector of the node falls between 0 and 45 degrees, a gradual change in color depth is presented. If it is greater than 45 degrees and does not exceed 90 degrees, another way of saturation transition is adopted and contour lines are marked around the node. Then, write the processed relationship between color and vector distribution into the visualization database of the corresponding grid node, and wait for overlay according to the grid node coordinate system of the 3D CAD model during rendering. When the colors and vectors of all nodes are rendered, the differences in heat recovery rate and the changes in heat flux density gradient in the energy input and output areas will be directly presented on the surface of the CAD model in different colors and lines. Finally, overlay the gradient direction and heat recovery rate correlation layer onto the grid node coordinate system of the 3D CAD model to generate the evaluation result of the heat recovery rate of thermal energy storage.

[0190] The above is only a preferred embodiment of the present invention, and it is not intended to limit the present invention in other forms. Any person skilled in the art may use the technical content disclosed above to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, as long as it does not depart from the technical solution content of the present invention, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis, characterized in that Including the following steps: Import the CAD geometric model of the thermal energy storage system, analyze the topological relationship, shape attributes, and component naming information of the geometric model, and obtain the original geometric feature data; based on the original geometric feature data, classify and identify the fluid domain, solid domain, energy storage material domain, heating channels, cooling channels, and insulation layer functional areas, and establish a TES geometric feature set; Based on the TES geometric feature set, retrieve material properties from a predefined library and assign them to the corresponding functional areas, set the mesh division parameters for the fluid boundary layer and the phase change region, obtain the material and mesh configuration parameters, and based on the channel features identified by the material and mesh configuration parameters and the TES geometric feature set, set the inlet and outlet boundary condition types, and generate an FEA preprocessing configuration file; Read the finite element analysis result file targeted by the FEA preprocessing configuration file, extract the temperature distribution data and the heat flux density vector data, and obtain the original FEA result data; based on the original FEA result data, map the temperature distribution data and the heat flux density vector data to the surface mesh nodes of the original CAD model to obtain the surface heat flux density gradient field; Based on the TES geometric feature set, determine the energy input surface, energy output surface, and the interface of the energy storage material domain, combine the mapped finite element analysis results, calculate the evaluation value of the thermal recovery rate of the thermal energy storage, and superimpose and display the surface heat flux density gradient field on the 3D CAD model through color mapping to obtain the evaluation result of the thermal recovery rate of the thermal energy storage.

2. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, wherein The steps for obtaining the original geometric feature data are as follows: Import the CAD geometric model of the thermal energy storage system, call the topological relationship information, component shape parameter information, and component naming structure information, generate the geometric boundary layer structure based on the mapping of the geometric body connection node structure, and divide the functional component combinations according to the naming structure to obtain the topological structure and the boundary association information of the component naming; According to the topological structure and the boundary association information of the component naming, identify the structural connections of adjacent geometric bodies on each boundary layer, perform geometric profile fitting of the node and edge combinations, and perform triangulation calculation on the connection areas in the order of functional combinations to generate the boundary feature lattice and the set of triangular elements of the functional areas; According to the boundary feature lattice and the set of triangular elements of the functional areas, sequentially extract the geometric dimension parameters, centroid coordinate positions, side length distributions, and curvature change information of each triangular element to generate the original geometric feature data.

3. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, wherein The steps for obtaining the TES geometric feature set are as follows: According to the original geometric feature data, extract the node coordinate sequences, boundary profile curvature values, normal direction distributions, boundary node densities, and regional surface areas of all closed boundary regions, group the boundary regions into an independent candidate region set based on the structural connection relationship, and obtain the preliminary regional structure division data; According to the preliminary regional structure division data, calculate the morphological distribution dispersion index of each candidate region, and the calculation formula is: Among them, D z is the morphological distribution dispersion index, A i is the surface area of the i-th candidate region, L i is the boundary node density of the i-th candidate region, is the average normal vector of the i-th candidate region, is the boundary center vector of the i-th candidate region, θ i is the maximum boundary normal angle of the i-th candidate region, and k is the number of candidate regions; According to the morphological distribution dispersion index, screen the candidate region set with different structural attributes, and combine the normal vector direction distribution and the node density change trend of each candidate region to perform functional area classification and labeling to generate the TES geometric feature set.

4. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, characterized in that The steps for obtaining the material and mesh configuration parameters are as follows: According to the TES geometric feature set, match the predefined material identifiers corresponding to each functional area label, extract the thermal conductivity, specific heat capacity, density, thermal conduction structure number, and surface smoothness level, and generate a regional material property parameter set; According to the regional material property parameter set, calculate the local mesh resolution control value, and the calculation formula is: Among them, R m is the local grid resolution control value, k t is the regional thermal conductivity, μ is the surface smoothness level, c p is the specific heat capacity, ρ is the density, γ is the heat conduction structure number, σ s is the maximum surface area of the region, λ s is the boundary gradient change ratio; According to the local mesh resolution control value, assign the local mesh resolution control value to the corresponding functional area, and respectively label the mesh accuracy of the fluid boundary layer and the mesh shape factor of the phase change area to form the material and mesh configuration parameters.

5. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, characterized in that The steps for obtaining the FEA preprocessing configuration file are as follows: According to the channel features identified by the material and mesh configuration parameters and the TES geometric feature set, locate the geometric boundaries corresponding to the heat flow inlet and heat flow outlet in each functional area, and based on the boundary node normal vector direction, node distribution density, and geometric area type, perform the identification and label mapping of the inlet unit and the outlet unit to generate the fluid channel boundary node attribute information; According to the fluid channel boundary node attribute information, extract the surface area, normal vector distribution, and mesh cell type corresponding to each inlet unit and outlet unit, and combine the material thermophysical properties and mesh cell distribution form in the material and mesh configuration parameters to set the boundary condition types at the inlets and outlets in each functional area, including the dominant mode and action dimension of the heat flux boundary or velocity boundary, to generate the boundary condition setting information; According to the boundary condition setting information, call the mesh division details and material thermodynamics properties of each unit area in the material and mesh configuration parameters to construct the node, element, and boundary term structures and constraint type syntax content of the finite element preprocessing input fields, and generate the FEA preprocessing configuration file.

6. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, wherein The steps for obtaining the original FEA result data are as follows: Read the finite element analysis result file targeted by the FEA preprocessing configuration file, retrieve the result index label corresponding to the analysis task, and locate the temperature distribution result path and heat flux density vector result path in the output file through the mapping field to generate a set of finite element result target paths; According to the set of finite element result target paths, sequentially load the temperature field and heat flux density field information within each time step or iteration step, and extract all temperature values and heat flux vector values according to the node numbers, and organize the extraction results into a structured matrix format to generate a temperature distribution matrix and a heat flux density vector matrix; According to the temperature distribution matrix and the heat flux density vector matrix, compare the mesh structure and node numbering rules set in the FEA preprocessing configuration file, perform matrix node order verification and boundary matching detection, process invalid numerical cells, and perform node structure binding to generate the original FEA result data.

7. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, characterized in that, The steps for obtaining the surface heat flux density gradient field are as follows: According to the original FEA result data, extract the surface mesh node indexes and geometric coordinate information of the original CAD model, and put the node coordinates in one-to-one correspondence with the temperature distribution data and the heat flux density vector data through the indexes to generate a combined attribute set of the mesh node thermophysics and geometric structure; According to the combined property set of the thermophysical and geometric structures of the grid nodes, calculate the heat flux density gradient intensity value of each node, and the calculation formula is as follows: Among them, G x is the intensity value of the heat flux density gradient, is the heat flux density vector of the i-th node, is the temperature gradient vector of the i-th node, is the node position vector, is the node normal vector, AG i is the area of the grid where the node is located, φ i is the angle between the node heat flux vector and the normal vector; According to the heat flux density gradient intensity value, perform value interpolation reconstruction on each grid node in the order of the numbers, and combine the geometric patch topological connection structure to construct a node heat flux gradient contour distribution map. Fill the field data at the spatial positions mapped according to the grid node coordinates on the surface of the original CAD model to generate a surface heat flux density gradient field.

8. The method for evaluating the heat recovery rate of thermal energy storage based on finite element analysis according to claim 1, wherein The steps for obtaining the evaluation result of the heat recovery rate of the thermal energy storage are as follows: According to the TES geometric feature set, identify the energy input surface, the energy output surface and the interface of the energy storage material domain, label the boundary edges through node indexing and the grid patch topological structure and bind the facing directions. At the same time, call the heat flux density vector and the node temperature field data in the original FEA result data to form a target surface thermophysical mapping structure; According to the target surface thermophysical mapping structure, calculate the evaluation value of the heat recovery rate of the thermal energy storage, and the calculation formula is as follows: where Ψ is the evaluation value of the heat recovery rate of thermal energy storage, is the surface heat flux density vector, is the unit vector in the surface normal direction, τ s is the nodal temperature change rate, ξ s is the average temperature gradient at the interface of the energy storage material domain, S in is the surface of energy input, S out is the surface of energy output, and t0 and t1 are the start and end boundaries of the preset time window; According to the evaluation value of the heat recovery rate of the thermal energy storage, combine the node heat flux direction distribution map in the surface heat flux density gradient field, and overlay the gradient direction and the heat recovery rate correlation layer to the grid node coordinate system of the 3D CAD model through color mapping to generate the evaluation result of the heat recovery rate of the thermal energy storage.

Images