Method for predicting heat treatment residual stress based on microstructure evolution

Abstract

The application belongs to the technical field of numerical simulation of material processing technology, and specifically discloses a heat treatment residual stress prediction method based on microstructure evolution, which comprises the following steps: obtaining a finite element grid and an initial physical field of a metal component with a blind hole and a flange; dividing a surface area based on a node characteristic thickness and a depth-diameter ratio, and respectively applying a differential convection heat transfer coefficient to form a heat transfer boundary; in a time step increment cycle, solving a temperature distribution by using the heat transfer boundary, calling a grain growth kinetics equation to update a real-time average grain diameter distribution; substituting the distribution into a grain size yield strength relationship to calculate a local real-time yield strength, combining a thermal strain to perform stress balance iteration to obtain a stress tensor; and repeating the above steps until a heat treatment is completed, and synchronously outputting a three-dimensional residual stress and a grain size distribution. The application realizes real-time embedded coupling of temperature, organization and stress, and improves the prediction accuracy.

Description

Technical Field

[0001] This invention belongs to the technical field of numerical simulation of material processing technology, and relates to a method for predicting residual stress in heat treatment based on microstructure evolution. Background Technology

[0002] The quenching and cooling process of metal components involves complex interactions between temperature fields, microstructure transformations, and stress-strain fields. Accurately predicting the residual stress distribution after heat treatment is crucial for controlling component deformation, preventing cracking, and ensuring service performance. Computer numerical simulation technology, by discretizing the continuous physical process into a finite number of computational nodes and recursively solving the problem along the time axis using heat conduction equations, grain growth kinetic models, and elastoplastic constitutive relations, has become an important tool for evaluating and optimizing heat treatment processes. Metal components with blind holes and flange structures exhibit abrupt spatial changes in their geometry. The differences in wall thickness in different regions and the semi-enclosed internal cavity structure lead to strong non-uniformity in the flow of the quenching medium and the heat transfer intensity, making the prediction of residual stress more technically challenging than for components with simpler shapes.

[0003] In existing technologies, the prediction of residual stress after heat treatment for components with complex geometries typically employs the finite element method. This method assigns a uniform heat transfer coefficient function to the entire outer surface of the component, which depends only on the surface temperature or the characteristic parameters of the quenching medium. The contribution of grain size to yield strength is often set once in the initial stage of stress calculation using empirical formulas, or a temperature-dependent yield strength degradation curve is used throughout the time history without considering the dynamic evolution of grain size. The heat transfer boundary conditions for semi-enclosed regions such as blind holes and cavities are usually simplified to the same convection conditions as those at openings, or handled by manually setting a constant reduction factor.

[0004] The aforementioned technical solutions suffer from several drawbacks. First, the uniform heat transfer coefficient function fails to reflect the fundamental difference in quenching cooling rates between thin-walled and thick-walled flange regions, leading to discrepancies between calculated and measured temperature histories of different parts of the component. Second, ignoring grain growth during the high-temperature stage of heat treatment results in the local yield strength of the material remaining constant or changing only with temperature throughout the cooling process. This contradicts the physical fact that the microstructure of the material continuously evolves during quenching, leading to distortion in the predicted plastic strain increment in stress tensor calculations. Third, using a constant reduction or equivalent heat transfer setting for the blind hole cavity fails to characterize the physical law of medium flow attenuation along the depth direction of the blind hole, causing the heat transfer coefficients at the bottom and wall of the blind hole to deviate from actual operating conditions, affecting the accuracy of residual stress prediction in this region. Summary of the Invention

[0005] In order to overcome the above-mentioned defects of the prior art and to achieve the above objectives, the present invention proposes the following technical solution: a heat treatment residual stress prediction method based on microstructure evolution, comprising: S1, obtaining a three-dimensional geometric model of a metal component with blind holes and flanges and dividing it into a finite element mesh, and mapping the obtained initial temperature field, initial grain diameter distribution and initial stress-strain field of the component to the finite element mesh nodes as the initial physical field.

[0006] S2. Based on the finite element mesh nodes, calculate the local characteristic thickness of each node on the outer surface and the local depth-to-diameter ratio of each node on the inner surface of the blind hole. Divide the outer surface nodes into thin-walled nodes and thick-walled flange nodes according to the local characteristic thickness, and mark the nodes on the inner surface of the blind hole.

[0007] S3. Based on the node division results, apply the thin-walled and thick-walled convective heat transfer coefficients determined by the corresponding heat transfer characteristic functions to the thin-walled and thick-walled flange nodes respectively. For the blind hole inner cavity surface nodes, the reference heat transfer coefficient is corrected by the attenuation function of the local depth-to-diameter ratio to form the heat transfer boundary.

[0008] S4. In the time step increment loop, the heat conduction equation is solved using the heat transfer boundary. Then, based on the current temperature and the grain diameter of the previous time step, the grain growth kinetic equation is called to calculate the grain growth rate and perform time integration to generate the real-time average grain diameter distribution of the current time step.

[0009] S5. Substitute the current grain diameter distribution into the grain size-yield strength relationship to calculate the local real-time yield strength of each node. Then, combine thermal strain and elastoplastic constitutive model to perform stress balance iteration to obtain the stress tensor of the current time step.

[0010] S6. Determine whether the cumulative solution time has reached the total heat treatment time. If not, use the current temperature distribution, real-time average grain diameter distribution, and stress tensor of the current time step as the initial conditions for the next time step, and return to the time step increment loop. If it has reached the time step, output the three-dimensional residual stress distribution and final grain size distribution of the metal component.

[0011] Compared with the prior art, the beneficial effects of the present invention are as follows: (1) The present invention automatically clusters the outer surface of the component into thin-walled areas and thick-walled flange areas according to the local characteristic thickness and applies the heat transfer characteristic function determined by inversion differently. At the same time, the convective heat transfer coefficient of the blind hole cavity is corrected by the exponential decay function of the local depth-to-diameter ratio, so that the spatial distribution of heat transfer boundary conditions can truly reflect the local influence of complex geometric features on the quenching cooling intensity.

[0012] (2) In each time step, the present invention performs temperature field solution, grain growth dynamics calculation, grain size-yield strength update and thermo-elastic-plastic stress balance iteration in the same time sequence, forming a real-time embedded iterative coupling between temperature history-driven micro-grain evolution and macro-mechanical response, thereby reducing the state lag error introduced by time mismatch in the traditional step-by-step decoupling method.

[0013] (3) The present invention directly transforms the real-time average grain diameter distribution into the local real-time yield strength that varies with time and space through the grain size-yield strength relationship, and embeds the yield strength into the yield criterion of the elastoplastic constitutive model, so that the local strength difference caused by the grain size gradient during the quenching process can immediately affect the evolution of plastic strain increment and stress tensor.

[0014] (4) The present invention simultaneously advances the solution of three types of physical fields, namely heat conduction, grain growth and stress balance, under a unified time step incremental cyclic framework, so that the final grain size distribution and three-dimensional residual stress distribution are output at the same time, ensuring the consistency of the microstructure evolution history and the residual stress formation history on the time axis. Attached Figure Description

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

[0016] Figure 1 This is a schematic diagram of the implementation steps of the method of the present invention.

[0017] Figure 2 This is a curve showing the attenuation of the convective heat transfer coefficient in the blind hole cavity of the present invention.

[0018] Figure 3 This is a graph showing the coupling relationship between grain size and local real-time yield strength in this invention.

[0019] Figure 4 This is a histogram of the residual stress tensor components at the end of the heat treatment process of this invention. Detailed Implementation

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

[0021] Please see Figure 1 As shown, the heat treatment residual stress prediction method based on microstructure evolution proposed in this invention includes: S1, obtaining a three-dimensional geometric model of a metal component with blind holes and flanges and dividing it into a finite element mesh, and mapping the obtained initial temperature field, initial grain diameter distribution and initial stress-strain field of the component to the finite element mesh nodes as the initial physical field.

[0022] Specifically, a three-dimensional geometric model of a metal component with blind holes and a flange structure is obtained. This metal component is a rotating part whose heat treatment deformation is to be predicted. The blind hole refers to a cylindrical internal cavity that is closed at one end and connected to the outside at the other end. The flange refers to a flange structure that extends circumferentially around the metal component. Together, they form a region of abrupt change in geometric features. The three-dimensional geometric model can be obtained by directly constructing it using computer-aided design software, or by reverse reconstructing it after acquiring the point cloud of the actual object's shape using a 3D scanning device. The model data includes complete topological information such as the inner surface of the blind hole, the end face of the flange, and the outer wall of the thin wall. After obtaining the three-dimensional geometric model, it is subjected to finite element mesh generation to generate finite element mesh nodes. Finite element mesh generation refers to discretizing a continuous geometry into a finite number of simple-shaped elements. The elements are interconnected through finite element mesh nodes, and each node carries a unique number and three-dimensional spatial coordinates. During mesh generation, for locations with drastic curvature changes around the blind hole and at the root of the flange, the mesh size is controlled so that the edge length of the elements in these regions does not exceed one-tenth of the local feature size, in order to accurately describe the geometric interface. The mesh type can be an eight-node hexahedral element or a four-node tetrahedral element. The final output is a node coordinate file and an element connection relationship file. It is ensured that the global spatial coordinate system of the current three-dimensional geometric model is consistent with the three-dimensional spatial coordinate system of the simulation model of the previous molding process, so as to eliminate rigid body displacement and rotation errors in the subsequent physical field mapping process.

[0023] After obtaining the finite element mesh nodes, it is necessary to acquire the initial temperature field, initial grain diameter distribution, and initial stress-strain field of the metal component at the end of the forming process before heat treatment. The forming process before heat treatment is typically a hot forming operation such as forging, rolling, or extrusion. At the end of this process, the internal temperature, microstructure, and mechanical state of the metal component constitute the starting conditions for subsequent heat treatment simulation. The initial temperature field refers to the temperature distribution at various points within the metal component, expressed in Kelvin or degrees Celsius. The initial grain diameter distribution refers to the average grain diameter at various points in the component, expressed in micrometers or meters. The initial stress-strain field refers to the stress tensor and equivalent plastic strain at various points in the component. The stress tensor contains six independent components, expressed in Pascals, and the equivalent plastic strain is dimensionless. These initial physical information can come from the finite element simulation results of the previous forming process, or they can be obtained by destructive temperature measurement, metallographic analysis and residual stress measurement of sampled parts manufactured under the same process conditions. Here, based on the measured temperature and grain size data of 200 sets of formed parts of this material, the interface heat transfer coefficient and material constitutive parameters in the forming simulation model are adjusted by the inverse optimization algorithm to complete the model calibration. Then, the calculation results of the calibrated model with a continuous spatial physical field distribution are exported as the initial field.

[0024] The initial temperature field, initial grain diameter distribution, and initial stress-strain field are mapped to the nodes of the finite element mesh to obtain the initial physical field conditions. The mapping process involves transferring these physical quantities, defined on the simulation mesh or discrete point cloud of the previous process, to the nodes of the current heat treatment simulation mesh. Here, a shape function interpolation method based on isoparametric inverse transformation is used to achieve the mapping. For each node of the current mesh, in the spatial domain formed by the simulation mesh of the previous process, a conventional spatial bounding box indexing algorithm is used to quickly locate the candidate source element containing the point. If the source element is an eight-node hexahedral element, a geometric mapping nonlinear equation system based on isoparametric transformation is constructed, that is, the spatial coordinates of the current node are equal to the sum of the products of the spatial coordinates of each node of the source element and their corresponding standard shape functions. The local coordinates of the point in the source element are obtained by solving the nonlinear equation system using a conventional numerical iteration method, with each component ranging from negative one to positive one. The node values ​​after the scalar field (e.g., temperature and grain diameter) mapping are calculated by weighted summation of the physical quantities stored in each node of the source element and the standard shape function corresponding to the element type. The stress tensor mapping is performed by directly applying the weighted summation interpolation scheme based on standard shape functions to each independent component, assuming the source mesh and the current mesh are in the same global coordinate system, to obtain the six-dimensional stress components of the current node. The equivalent plastic strain is also processed using scalar interpolation. After the above node-by-node interpolation operation, each finite element mesh node of the current heat treatment simulation mesh is assigned temperature, grain diameter, stress components, and plastic strain values. The set of these node values ​​constitutes the initial physical field conditions required for subsequent calculations. Shape function interpolation itself has internal element continuity and node value consistency; its dimensions are determined by the transferred physical quantities. The interpolation operation itself does not introduce additional dimensional changes, thus ensuring the dimensional consistency of physical quantities before and after mapping.

[0025] For example, a metal component made of alloy structural steel (e.g., grade 40CrNiMoA) is used as the specific implementation object. This component has a blind hole and flange structure. The blind hole diameter is 20 mm, the blind hole depth is 50 mm, the flange outer diameter is 100 mm, the flange thickness is 25 mm, and the nominal thickness of the thin-walled region is 8 mm. A 3D geometric model of the component is obtained using 3D modeling software and output in a standard geometry exchange format. This geometric model is then read into finite element preprocessing software, and finite element meshing is performed using second-order tetrahedral elements, generating a total of 156,243 finite element mesh nodes and 98,520 elements. By checking the mesh quality indicators, the side length ratio of all elements is less than 5, and the minimum Jacobian determinant is greater than 0.7, meeting the calculation accuracy requirements.

[0026] The previous forming process was closed-die forging. After the forging was completed, the metal component was ejected from the die and air-cooled to room temperature. This moment at the end of air cooling was taken as the initial state. Using a thermo-mechanical coupled finite element model calibrated with grain size measurement data under 200 sets of heating and deformation conditions of this material, the initial temperature field, initial grain diameter distribution, and initial stress-strain field at this moment were obtained and stored as a source mesh. The source mesh contains a total of 132,187 nodes and 81,654 elements. Now, the finite element mesh node numbered 35421 in the target mesh is extracted, and its spatial coordinates are... All units are millimeters. Spatial location was performed using a spatial bounding box indexing algorithm within the source mesh's element set. The point was detected as being contained within an eight-node hexahedral element in the source mesh. The eight nodes of this source element are numbered as follows: to The local coordinate inverse equation of the geometric mapping is constructed. Initial guesses are set to zero, and the local coordinates of the point within the source element are obtained through numerical iteration. , , According to the shape function calculation rules, the shape function values ​​corresponding to each node are as follows: , , , , , , , The sum is 1.000, which satisfies the normalization property of shape functions.

[0027] The physical quantities stored at each node of the source element are extracted: the temperatures of the eight nodes are 781, 784, 779, 776, 765, 768, 761, and 758 degrees Celsius; the grain diameters are 24.8, 25.1, 24.5, 24.2, 23.6, 23.9, 23.2, and 23.0 micrometers; and the normal stress components in the x-direction of the stress tensor are 15.2, 14.8, 13.1, 16.5, 10.4, 9.8, 11.3, and 12.7 MPa. The mapped node physical quantities are calculated by weighted summation of shape functions: temperature is 768.9 degrees Celsius; grain diameter is 23.85 micrometers; and the normal stress component in the x-direction is 12.61 MPa. The remaining stress components and equivalent plastic strain are obtained by interpolation in the same manner. Thus, node 35421 has obtained complete mapped values. The above positioning and interpolation process was repeated for all 156,243 nodes to generate the initial temperature field distribution, initial grain diameter distribution, and initial stress-strain field of the entire metal component. The physical quantities on these finite element mesh nodes constitute the initial physical field conditions, which are stored in a solver-readable binary file for subsequent time step increment cyclic calls. Through this closed-loop verification mapping operation process, the reproducibility of the complete technical process of three-dimensional geometric model, finite element mesh generation, initial field acquisition, and the formation of initial physical field conditions based on shape function mapping was verified.

[0028] S2. Based on the finite element mesh nodes, calculate the local characteristic thickness of each node on the outer surface and the local depth-to-diameter ratio of each node on the inner surface of the blind hole. Divide the outer surface nodes into thin-walled nodes and thick-walled flange nodes according to the local characteristic thickness, and mark the nodes on the inner surface of the blind hole.

[0029] In a preferred embodiment, based on finite element mesh nodes, the local characteristic thickness of each node on the outer surface and the local depth-to-diameter ratio of each node on the blind hole cavity surface are calculated. The outer surface nodes are divided into thin-walled region nodes and thick-walled flange region nodes according to their local characteristic thickness, and the blind hole cavity surface nodes are marked. This includes: extracting the finite element mesh nodes of the outer surface of the three-dimensional geometric model and obtaining the surface unit normal vector at each outer surface node; taking each outer surface node as the starting point, drawing a ray along its surface normal vector direction into the component, finding the intersection point of the ray with the opposite outer surface of the component, calculating the spatial distance between the starting point and the intersection point, and using the spatial distance as the local characteristic thickness of that outer surface node; and further... The finite element mesh nodes of the blind hole cavity are selected, and the projection distance from each node on the surface of the blind hole cavity along the blind hole axis to the blind hole opening plane is calculated as the node depth. The blind hole diameter is obtained, and the ratio of the node depth to the blind hole diameter is used as the local depth-to-diameter ratio of the node. The local feature thickness of all the calculated outer surface nodes is statistically clustered to identify two thickness cluster centers representing thin-walled structures and thick-walled flange structures, respectively. The classification thresholds for thin-walled and thick-walled flange regions are determined based on the thickness cluster centers. Nodes with local feature thickness less than the classification threshold are classified as thin-walled region nodes, and nodes with local feature thickness not less than the classification threshold are classified as thick-walled flange region nodes.

[0030] Specifically, following the finite element mesh nodes and pre-defined 3D geometric model obtained in step S1, the core task of step S2 is to identify and classify the geometric features of the heat transfer surface of the metal component, distinguishing areas with significantly different physical cooling behaviors. The geometric features to be processed include local feature thickness and local depth-to-diameter ratio. Local feature thickness refers to the shortest spatial distance along the wall thickness direction from any point on the outer surface of the component to the other outer surface, reflecting the thickness of the solid material involved in heat conduction at that point. Local depth-to-diameter ratio refers to the ratio of the depth of any point on the inner surface of the blind hole from the opening of the blind hole along the axial direction to the cylindrical inner diameter of the blind hole. This dimensionless quantity characterizes the degree of attenuation of the medium flow and heat transfer capacity at the depth of that point within the blind hole.

[0031] During implementation, the outer surface nodes of the 3D geometric model are first extracted from the finite element mesh nodes. The outer surface refers to all boundaries of the metal component that are in direct contact with the quenching cooling medium, including the visible cylindrical surfaces, flange end faces, and planar rings near the blind hole openings, but excluding the hole walls and bottoms within the blind hole cavity. The outer surface nodes are extracted by determining the facet affiliation of the node to the element. Specifically, all faces located on the surface of the solid element are traversed, and the vertex numbers on each facet are collected and deduplicated to form the outer surface node set. Let the total number of nodes in this outer surface node set be denoted as . .

[0032] For each node in the set of outer surface nodes, obtain the surface unit normal vector at that node. The surface unit normal vector is a unit-length vector perpendicular to the outer surface and pointing outwards from the component at the node location. The normal vector at each node is obtained by area-weighted averaging of the normal vectors of all outer surface patches sharing that node. Specifically, the unit normal vectors of each patch are weighted by their area, and the resulting vector is normalized to obtain the surface unit normal vector at that node. The advantage of area-weighted averaging is that it reduces the interference of smaller or distorted patches on the normal direction, making the starting direction of ray tracing more stable.

[0033] Next, starting from each outer surface node, a ray is drawn into the component along its surface normal vector direction. The intersection point of the ray with the opposite outer surface of the component is then determined. Here, a ray refers to an infinite straight line extending from the starting point along the negative normal vector direction, corresponding to the direction pointing towards the inner cavity of the component or the opposite outer wall. The opposite surface refers to the first outer surface encountered after penetrating the entire solid from the node along the ray direction. The intersection operation is performed using a conventional ray-to-grid intersection algorithm. All outer surface patches are traversed, excluding the patch containing the starting point and common neighbor patches with too small an angle. For each candidate patch, the intersection point of the ray with the plane containing the triangle is calculated, and then it is determined whether the intersection point is located inside the triangle. If a valid intersection point exists, its spatial coordinates are recorded.

[0034] Calculate the spatial distance between the starting point and the intersection point. This spatial distance is defined as the local feature thickness of the outer surface node. Let the local feature thickness of node i be denoted as . Its calculation formula is Euclidean distance: ; In the formula, Let be the local feature thickness of node i, in meters; Let i be the spatial coordinates of the outer surface node i; Let be the coordinates of the intersection point of the ray and the opposite outer surface.

[0035] Parallel to the processing of these surface nodes is the calculation of the geometric features of the surface nodes of the blind hole cavity. The surface nodes of the blind hole cavity refer to the finite element mesh nodes located on the cylindrical wall and bottom surface of the blind hole. The extraction method is based on geometric location; nodes whose cylindrical radius is close to half the diameter of the blind hole and whose axial position is between the hole opening plane and the hole bottom surface are all included in this set. Here, the blind hole opening plane refers to the plane containing the annular end face of the blind hole connecting to the external environment, and the axial direction of the blind hole is the centerline direction of the cylindrical cavity. These geometric references have already been constructed in the three-dimensional geometric model of step S1.

[0036] For each node in the set of nodes on the inner surface of the blind hole, calculate the projected distance of that node along the blind hole axis to the blind hole opening plane as the node depth. Let the unit vector along the blind hole axis be... The coordinates of any reference point on the plane of the blind hole opening are: The coordinates of the position of node k on the surface of the inner cavity are: Then the node depth for: ; In the formula, The depth of node k is in meters; and All are three-dimensional position vectors, with units of meters; It is a unit vector, dimensionless; taking the absolute value ensures that the depth is non-negative, because the bottom node of the blind hole is located in the positive projection direction of the opening plane.

[0037] The blind hole diameter is obtained, referring to the diameter of the cylindrical cross-section of the blind hole's inner cavity. This diameter is derived from the 3D geometric model design parameters in step S1 and is determined during the manufacturing and modeling stages. The ratio of the node depth to the blind hole diameter is used as the local depth-to-diameter ratio of that node. Since the unit of the node depth and the hole diameter are the same, their ratio is a dimensionless pure number. When the node is located at the root of the blind hole opening, the depth is close to zero, and the local depth-to-diameter ratio is close to zero; when the node is located at the bottom of the blind hole, the depth is approximately equal to the hole depth, and the local depth-to-diameter ratio reaches the maximum value of the depth-to-diameter ratio of this blind hole.

[0038] After calculating the local feature thickness of all outer surface nodes, these thickness values ​​need to be statistically clustered to automatically identify which outer surface regions belong to thin-walled structures and which belong to thick-walled flange structures. Statistical clustering refers to an unsupervised learning method that divides a dataset into two subclasses based on the density of the data. Here, a conventional statistical clustering algorithm, such as K-means clustering, is used, with the number of clusters set to 2. Initially, one cluster center is set to the calculated minimum feature thickness value, and the other to the maximum feature thickness value. Then, the following two steps are iteratively executed until the changes in the positions of the two cluster centers are less than... Meters: First, assign each outer surface node to the category of the cluster center that is closest to its local feature thickness value; second, within each category, take the arithmetic mean of the local feature thickness values ​​of all nodes belonging to that category, and update the cluster center of that category.

[0039] After clustering convergence, two thickness cluster centers are obtained, denoted as . and ,in The natural result of clustering is that smaller numerical values ​​form cluster centers. This represents the central tendency of wall thickness categories in thin-walled structural regions, with larger cluster centers. This represents a concentration trend in thick-walled structures such as flanges. The classification threshold between thin-walled structures and thick-walled flange structures is denoted as... The thickness is determined by averaging two centers of concentration. ; In the formula, The classification threshold is in meters. The thickness concentration center corresponding to the thin-walled structure is in meters; This represents the thickness clustering center for the thick-walled flange structure, measured in meters. This threshold is located in the physical space at the midpoint between two characteristic thicknesses, ensuring that the two types of nodes are linearly separated to the greatest extent possible.

[0040] Next, node partitioning is performed: external nodes with local feature thickness less than the classification threshold are classified as thin-walled zone nodes, and nodes with local feature thickness not less than the classification threshold are classified as thick-walled flange zone nodes. Separation markers are written to the integer attribute array corresponding to each node for later use in step S3. Simultaneously, the blind hole inner cavity surface nodes previously filtered by geometric location are independently recorded and marked, forming three sets of boundary nodes: thin-walled zone node set, thick-walled flange zone node set, and blind hole inner cavity surface node set. Thus, step S2 completes the geometric feature quantization and automatic classification and marking of all heat exchange surface nodes.

[0041] For example, using the alloy structural steel metal component given in step S1, the blind hole diameter is 20 mm, the blind hole depth is 50 mm, the flange outer diameter is 100 mm, the flange thickness is 25 mm, and the nominal thickness of the thin-walled region is 8 mm. The finite element mesh has 156,243 nodes, of which 42,518 are on the outer surface and 3,841 are on the inner surface of the blind hole cavity.

[0042] Node number 12087 from the outer surface node set is retrieved for tracing. This node is located on the outer side of the thin-walled cylinder, with spatial coordinates x=18.215, y=7.983, z=42.630 (all in millimeters). This node is shared by six triangular outer surface patches with areas of 1.203, 1.187, 1.224, 1.192, 1.209, and 1.199 square millimeters, respectively. The unit normal vector of each patch has been calculated and normalized from the patch edge vectors. After area-weighted summation and normalization, the surface unit normal vector of this node is calculated, with a magnitude of 1.0. A ray is drawn from this node along the negative normal vector direction, i.e., the vector (-0.412, -0.781, -0.470), into the component.

[0043] The ray intersects with all triangular faces on the outer surface of the component. After excluding the face containing the starting point and its adjacent faces sharing edges with the starting point, the ray-triangle intersection algorithm is performed on each of the 13427 candidate faces. An intersection point is detected between the ray and face number 35891. The coordinates of the intersection point are... , , The intersection point is located on the surface of the inner hole of the thin-walled cylinder. By calculating the spatial Euclidean distance between the starting point and the intersection point, the local feature thickness of this outer surface node is found to be 7.602 mm. The deviation from the nominal thickness of 8 mm comes from the fillet transition and mesh discretization error at this location, which is within a reasonable range.

[0044] After performing ray intersection and distance calculations on all 42,518 outer surface nodes, a local feature thickness dataset was obtained. K-means clustering was used, with initial cluster centers set to the minimum value of 4.532 mm and the maximum value of 49.871 mm in the dataset. After 12 iterations, convergence was achieved, yielding a thickness cluster center of 7.843 mm for thin-walled structures and 24.672 mm for thick-walled flange structures. The classification threshold was 16.258 mm. Traversing all outer surface nodes, 31,206 nodes with a local feature thickness less than 16.258 mm were classified as thin-walled region nodes, and 11,312 nodes with a local feature thickness of 16.258 mm or more were classified as thick-walled flange region nodes. The classification results showed that thin-walled region nodes were concentrated on the thinning areas of the outer wall of the cylinder and the inner edge of the flange, while thick-walled flange region nodes were concentrated on the thicker parts of the flange end face and sides. The partition boundaries closely matched the actual geometric transition zones of the components.

[0045] The process then shifts to processing the surface nodes of the blind hole cavity. Node number 71592 is selected, located on the inner wall of the cylinder within the blind hole. Its spatial coordinates are x=9.882, y=4.126, z=13.574 mm. The unit vector along the blind hole axis is obtained through cylinder fitting, meaning it is parallel to the z-axis of the global coordinate system. The opening plane of the blind hole is selected as the average plane passing through all nodes on the outer end face of the opening, and a reference point is chosen on this plane. The projected distance from this node to the opening plane is calculated to be 13.574 mm. The blind hole diameter is 20.0 mm, and the local depth-to-diameter ratio is 0.679, indicating that the node is located at approximately two-thirds of the blind hole depth. The local depth-to-diameter ratio is calculated for all 3841 surface nodes of the blind hole cavity, with values ​​ranging from 0 to 2.5. The local depth-to-diameter ratio of the bottom nodes, which are close to 50 mm deep, reaches approximately 2.5. The integer marker arrays of the nodes on the inner surface of the blind hole are set independently, and finally, the markers of all three types of nodes are written into the binary interface file. The above data and calculation process fully verify that the technical path of obtaining local feature thickness through the ray distance method, obtaining local depth-to-diameter ratio through the depth-to-diameter ratio, and determining the classification threshold by using the thickness aggregation center to achieve automatic node division can be reproduced.

[0046] S3. Based on the node division results, apply the thin-walled and thick-walled convective heat transfer coefficients determined by the corresponding heat transfer characteristic functions to the thin-walled and thick-walled flange nodes respectively. For the blind hole inner cavity surface nodes, the reference heat transfer coefficient is corrected by the attenuation function of the local depth-to-diameter ratio to form the heat transfer boundary.

[0047] In a preferred embodiment, thin-walled and thick-walled convective heat transfer coefficients determined by corresponding heat transfer characteristic functions are applied to thin-walled nodes and thick-walled flange nodes, respectively. This includes: acquiring measured temperature change data of thin-walled geometric feature samples during quenching and cooling; obtaining a thin-walled heat transfer characteristic function that varies with surface temperature by solving an inverse heat conduction problem; acquiring measured temperature change data of thick-walled geometric feature samples during quenching and cooling; obtaining a thick-walled heat transfer characteristic function that varies with surface temperature by solving an inverse heat conduction problem; for each thin-walled node, acquiring the surface temperature of that node at the current time step; and calculating the thin-walled convective heat transfer coefficient using the thin-walled heat transfer characteristic function; for each thick-walled flange node, acquiring the surface temperature of that node at the current time step; and calculating the thick-walled convective heat transfer coefficient using the thick-walled heat transfer characteristic function.

[0048] In a further preferred embodiment, the reference heat transfer coefficient of the blind hole inner cavity surface node is corrected using a decay function of the local depth-to-diameter ratio, including: obtaining the convective heat transfer coefficient of the node at the blind hole opening and determining it as the reference heat transfer coefficient; constructing a monotonically decreasing exponential decay factor function with the local depth-to-diameter ratio as the independent variable, such that the decay factor function takes a value of one when the local depth-to-diameter ratio is zero, and gradually decays as the local depth-to-diameter ratio increases; for each blind hole inner cavity surface node, inputting its corresponding local depth-to-diameter ratio into the decay factor function to obtain the decay factor, and calculating the product of the reference heat transfer coefficient and the decay factor as the blind hole inner cavity convective heat transfer coefficient of that node.

[0049] Specifically, after completing the division and marking of thin-walled region nodes, thick-walled flange region nodes, and blind hole inner cavity surface nodes in step S2, the objective of step S3 is to assign convective heat transfer boundary conditions reflecting their local geometric characteristics and differences in medium flow to these surface nodes of different categories. The convective heat transfer boundary condition refers to the quantitative relationship between the heat flux density transferred from the surface node to the cooling medium and the surface temperature and medium temperature of that node in the heat conduction equation, determined by the convective heat transfer coefficient. The convective heat transfer coefficient is a dimensional physical quantity, measured in watts per square meter per Kelvin, characterizing the rate of heat exchange per unit area between the solid surface and the contacting quenching medium under a unit temperature difference.

[0050] The convective heat transfer coefficients of the thin-walled and thick-walled flange nodes depend on pre-calibrated thin-walled and thick-walled heat transfer characteristic functions, respectively. The heat transfer characteristic function defines the functional mapping relationship between the surface heat transfer coefficient and the real-time surface temperature at that point. Its calibration process is based on quenching experiments using samples made of the same material as the metal component but with clearly different geometric wall thicknesses. The thin-walled geometric characteristic sample refers to a simple plate or cylindrical specimen with a single wall thickness corresponding to the thickness concentration center of the thin-walled structure obtained in step S2; here, it is taken as 8 mm. The thick-walled geometric characteristic sample refers to a block or thick-walled cylindrical specimen with a wall thickness corresponding to the thickness concentration center of the thick-walled flange structure; here, it is taken as 25 mm. Both types of samples are made of alloy structural steel (e.g., grade 40CrNiMoA), consistent with the analyzed metal component.

[0051] To obtain measured temperature change data of thin-walled geometric feature samples during quenching and cooling, a 1.5 mm diameter armored K-type thermocouple was embedded along the thickness direction at the geometric center of the thin-walled sample, with the measuring point located at half the depth of the sample wall thickness. The sample was heated to 850 degrees Celsius in a resistance furnace and held at that temperature for 30 minutes, then rapidly transferred to 25 degrees Celsius rapid quenching oil for cooling. The electromotive force signal of the thermocouple was continuously acquired using a paperless recorder at a sampling frequency of 10 Hz and converted into temperature values, recorded for 600 seconds, thus obtaining a sequence of temperature changes at the core of the thin-walled sample over time. Similarly, for thick-walled geometric feature samples, under the same heat treatment and quenching conditions, a thermocouple was placed at half the wall thickness, and measured temperature change data for the same duration were obtained.

[0052] After obtaining the measured temperature change data, the thin-walled heat transfer characteristic function, which describes the change of the heat transfer coefficient with surface temperature, is obtained by solving the inverse problem of heat conduction. The inverse problem of heat conduction refers to a mathematical method that inversely calculates the surface heat transfer boundary conditions based on the temperature response at a temperature measurement point inside the solid. Its core is to establish a forward problem model describing the one-dimensional unsteady heat conduction inside the sample, using the surface heat transfer coefficient as the unknown to be optimized, and iterating to make the simulated temperature approximate the measured temperature. The governing equation for the forward problem of the thin-walled sample adopts a conventional one-dimensional unsteady partial differential equation of heat conduction, which describes the energy conservation law of heat transfer within the solid via conduction.

[0053] When solving the inverse problem, the surface heat transfer coefficient is used. Consider as surface temperature The function is parameterized as a piecewise linear function controlled by the heat transfer coefficients at several characteristic temperature nodes. The parameter vector to be optimized is denoted as... , corresponding to Q equally spaced characteristic surface temperature points. The forward problem is discretized on the time and space axes using the finite difference method, given a set of... Based on the current estimated value, the corresponding heart temperature curve can be calculated. The goal of the inversion is to adjust... To minimize the deviation between the simulated and measured temperatures, the deviation is measured by the root mean square error function: ; In the formula, The objective function value is expressed in Kelvin. The number of data points included in the measured temperature time series is dimensionless. Let K be the core temperature calculated by the forward problem model at the p-th time point, in Kelvin. Let be the measured temperature at time point p, in Kelvin. The optimization algorithm uses the Levenberg-Marquardt method, updating parameters until the objective function value falls below 0.5 Kelvin or the relative parameter change is less than one percent. After inversion convergence, the heat transfer coefficient values ​​at the discrete characteristic temperature nodes are obtained. Connecting them using linear interpolation yields the thin-walled heat transfer characteristic function. The function uses surface temperature As input, output the thin-walled convective heat transfer coefficient value at the corresponding temperature.

[0054] The heat transfer characteristic function of thick-walled geometric feature samples follows the same technical process. Simply change the geometric thickness in the forward problem model to 25 mm and use the measured temperature change data of the thick-walled sample as the inversion target to finally obtain the thick-walled heat transfer characteristic function. It is also stored in an interpolated table.

[0055] Within a specific time-step increment loop, when it is necessary to apply the convective heat transfer coefficient to the thin-walled nodes, for each node in the set of thin-walled nodes, the surface temperature of that node at the current time step is first obtained. This temperature is updated from the solution of the heat conduction equation in step S4. Then, using... Using the independent variable, query the thin-walled heat transfer characteristic function. Table, if If the value falls between two table nodes, then linear interpolation is used to derive the thin-walled convective heat transfer coefficient. Similarly, for nodes in the thick-walled flange area, their surface temperature is taken. The thick-walled convective heat transfer coefficient can be obtained by interpolation from the table. At this point, the heat transfer coefficients of all nodes on the outer surface of the component have been assigned values ​​corresponding to the local wall thickness classification.

[0056] For the surface nodes within the blind orifice cavity, their heat transfer environment differs significantly from that of the outer surface exposed to external forced convection due to their semi-enclosed geometry. The internal medium velocity decreases rapidly with increasing depth, thus requiring a correction method to describe this decrease. First, a baseline heat transfer coefficient is determined. The nodes at the blind orifice opening are several outer surface nodes located at the edge of the orifice opening. These nodes belong to the outer surface set and are adjacent to the blind orifice inlet; their heat transfer conditions represent the baseline heat transfer capacity provided by the external flow at that moment. The baseline heat transfer coefficient is obtained by arithmetically averaging the convective heat transfer coefficients assigned to these nodes at the current time step. ; In the formula, The reference heat transfer coefficient is expressed in watts per square meter per Kelvin. This represents the total number of nodes at the blind hole opening; is the convective heat transfer coefficient calculated for the j-th opening node according to its assigned zone, in watts per square meter per Kelvin.

[0057] Subsequently, a monotonically decreasing exponential decay factor function with the local depth-to-diameter ratio as the independent variable was constructed. This decay function needs to reflect the physical law that the heat transfer intensity decreases continuously along the depth direction after the medium enters the blind hole. Its mathematical form must satisfy the following condition: the decay is zero and takes the value of one at the local depth-to-diameter ratio (i.e., at the blind hole opening), and it decreases smoothly as the local depth-to-diameter ratio increases. The constructed decay factor function is as follows: ; In the formula, Let be the attenuation factor at node k, which is dimensionless; Let be the local depth-to-diameter ratio of node k, which is dimensionless and is calculated in step S2. It is an exponential function with the natural constant e as its base; The attenuation coefficient is a positive dimensionless constant used to control the degree of attenuation. The setting is based on the nonlinear least squares fitting of experimental data of the inner wall Nusselt number of 30 blind hole samples with different depth-to-diameter ratios under similar quenching conditions. When the value is 1.2, the sum of squares of the residuals of the ratio curve of the attenuation factor function and the measured dimensionless heat transfer capacity is minimized, while satisfying the physical constraint that the heat transfer coefficient at the bottom of the hole drops to less than 10% of that at the opening.

[0058] For each node in the set of nodes on the surface of the blind hole cavity, the local depth-to-diameter ratio calculated in step S2 is used. Substituting into the attenuation factor function, the attenuation factor is obtained. This node is calculated by multiplying the reference heat transfer coefficient by the attenuation factor. Since the reference heat transfer coefficient... The dimensions of the coefficient are consistent with those of the convective heat transfer coefficient. After multiplying by the dimensionless attenuation factor, the dimensions of the result remain correct. All blind hole internal surface nodes are assigned their respective convective heat transfer coefficients through this step. Together with the convective heat transfer coefficients on the thin-walled and thick-walled flange nodes, these form the complete nodal heat transfer boundary conditions for the entire metal component during the quenching and cooling process. These boundary conditions will be used as the third type of boundary conditions when solving the heat conduction equation, and will be input into the temperature field calculation at each time step of step S4.

[0059] For example, based on the alloy structural steel metal component and its node division results given in steps S1 and S2, the specific implementation and numerical calculation process of step S3 are demonstrated. The thin-walled geometric feature sample uses a flat plate with dimensions of 200 mm long, 200 mm wide, and 8 mm thick, with a material density of 7850 kg / m³. Specific heat capacity and thermal conductivity are retrieved from the material database as functions of temperature. The quenching oil temperature is 25°C, and the sample is cooled starting at 850°C. The thermocouple-recorded core temperature change data contains 6000 time points with a time interval of 0.1 seconds. The forward problem model uses a finite difference scheme with a spatial step of 0.2 mm and a time step of 0.01 seconds. The inverse problem feature temperature nodes are selected from 750°C to 50°C, with 15 nodes spaced 50°C apart. After 71 iterations, the objective function decreases from the initial 127.3 Kelvin to 0.37 Kelvin, satisfying the convergence condition. In the obtained thin-walled heat transfer characteristic function table, the heat transfer coefficient corresponding to a surface temperature of 700 degrees Celsius is 1852.7 watts per square meter per Kelvin, and the heat transfer coefficient corresponding to a surface temperature of 300 degrees Celsius is 1438.3 watts per square meter per Kelvin. The geometric dimensions of the thick-walled geometric feature sample are 200 mm long, 200 mm wide, and 25 mm thick. After inversion, the thick-walled heat transfer characteristic function is obtained, where the heat transfer coefficient is 1721.5 watts per square meter per Kelvin at a surface temperature of 700 degrees Celsius, and 1325.8 watts per square meter per Kelvin at a surface temperature of 300 degrees Celsius.

[0060] At the current time step, with a step size of 0.5 seconds, the surface temperature of node 12087 in the thin-walled region is 634.7 degrees Celsius. In the thin-walled heat transfer characteristic function table, this temperature value lies between the characteristic temperature nodes 650 degrees Celsius and 600 degrees Celsius. The heat transfer coefficient corresponding to 650 degrees Celsius is 1812.4 watts per square meter per Kelvin, and the heat transfer coefficient corresponding to 600 degrees Celsius is 1765.9 watts per square meter per Kelvin. Through linear interpolation, the thin-walled convective heat transfer coefficient is calculated to be 1798.17 watts per square meter per Kelvin. Node 12087 is therefore assigned a thin-walled convective heat transfer coefficient of 1798.17 watts per square meter per Kelvin. The current surface temperature of node number 22341 in the thick-walled flange area is 621.3 degrees Celsius. Similarly, by referring to the thick-walled heat transfer characteristic function table and interpolating, the thick-walled convective heat transfer coefficient is obtained as 1673.5 watts per square meter per Kelvin.

[0061] There are 76 nodes at the blind orifice opening. The convective heat transfer coefficients assigned to these nodes at the current time step are 1795.22, 1801.65, 1796.84, ... , and their arithmetic mean yields a baseline heat transfer coefficient of 1798.34 watts per square meter per Kelvin. Now, node 71592 on the surface of the blind orifice cavity is processed. Its local depth-to-diameter ratio of 0.679 has already been calculated in step S2. The attenuation coefficient is 1.2, and substituting it into the attenuation factor function yields an attenuation factor of 0.4428. The convective heat transfer coefficient of this node within the blind orifice cavity is the product of 1798.34 and 0.4428, which is 796.33 watts per square meter per Kelvin. Next, consider a node located at the bottom of the blind hole, with a local depth-to-diameter ratio of 2.498. The calculated attenuation factor is 0.0499, and the calculated convective heat transfer coefficient is 89.74 watts per square meter per Kelvin. This value is only about five percent of the baseline heat transfer coefficient at the opening, demonstrating the weak heat transfer characteristics at the bottom of the blind hole. Please refer to [link / reference]. Figure 2 As shown, Figure 2 This is a graph showing the decay curve of the convective heat transfer coefficient inside the blind hole cavity of this invention. As shown in the figure, the correction of the convective heat transfer boundary inside the blind hole cavity is achieved by relying on an exponential decay function. For the semi-closed blind hole structure, the reference heat transfer coefficient at the opening is extracted, and the local depth-to-diameter ratio of the node is used as the input decay factor function as the independent variable. Subsequently, through multiplication, the convective heat transfer coefficient deep inside the blind hole decreases exponentially, realistically depicting the physical phenomenon that the heat transfer capacity of the quenching medium weakens sharply with increasing depth inside the blind hole. After all the convective heat transfer coefficients of the three types of nodes are generated, they are written into the solver in the form of node attribute files as the heat transfer boundary input for the heat conduction solution in step S4, thus completing all the operations in step S3.

[0062] S4. In the time step increment loop, the heat conduction equation is solved using the heat transfer boundary. Then, based on the current temperature and the grain diameter of the previous time step, the grain growth kinetic equation is called to calculate the grain growth rate and perform time integration to generate the real-time average grain diameter distribution of the current time step.

[0063] In a preferred embodiment, the grain growth kinetics equation is invoked based on the current temperature and the grain diameter of the previous time step to calculate the grain growth rate and perform time integration, generating the real-time average grain diameter distribution of the current time step. This includes: obtaining the grain growth kinetic parameters of the metal component material, including the grain boundary migration activation energy and the rate constant; for each finite element mesh node, based on the current temperature and grain growth kinetic parameters of the node, calculating the grain boundary migration rate coefficient of the node at the current temperature; substituting the grain boundary migration rate coefficient and the grain diameter of the node in the previous time step into the grain growth kinetics model to obtain the instantaneous grain growth rate of the node; numerically integrating the instantaneous grain growth rate within the step size of the current time step to obtain the grain diameter increment; adding the grain diameter increment to the grain diameter of the previous time step to obtain the real-time average grain diameter of the node in the current time step.

[0064] Specifically, after constructing the node heat transfer boundary conditions for all surface nodes in step S3, step S4 enters the actual execution phase of the time step increment loop. The time step increment loop refers to discretizing the continuous heat treatment cooling process on the time axis into a recursive calculation structure between a finite number of cells, with the length of each cell being the time step. The physical state of the current time step is obtained by adding the physical changes within that step to the state of the previous time step, and this process is repeated segment by segment until the entire heat treatment process is covered. Step S4 requires two core calculations within a time step: first, solving the heat conduction equation using the node heat transfer boundary conditions to obtain the current temperature distribution of each node in the current time step; second, using the current temperature distribution and the grain diameter distribution of the previous time step as input, calling the grain growth kinetics equation to calculate the grain growth rate and performing time integration to generate the real-time average grain diameter distribution of the current time step.

[0065] The heat conduction equation describes the conservation relationship of temperature within a solid as a function of space and time. For the three-dimensional unsteady heat conduction problem, the continuous temperature field is represented as an interpolated set of temperatures at the nodal points of the finite element mesh using the conventional finite element spatial discretization method. Combined with the nodal heat transfer boundary conditions generated in step S3, a semi-discrete set of ordinary differential equations containing the heat capacity matrix and the thermal stiffness matrix is ​​constructed.

[0066] Subsequently, a conventional implicit time integration algorithm (such as the two-time-layer implicit trapezoidal rule) is used at time steps. The above equations are discretized. This implicit scheme is unconditionally stable. For the initial stage of quenching and cooling when the temperature gradient is sharp, a time step of 0.1 seconds is used to ensure convergence of the nonlinear solution. When the surface temperature of the component drops below 50 degrees Celsius of the medium temperature, the time step can be gradually increased to 1.0 second to improve computational efficiency. Since the thermal properties of the material (such as thermal conductivity and specific heat capacity) and the convective heat transfer coefficient change with temperature, the coefficient matrix of the equations has nonlinear characteristics. Therefore, conventional nonlinear iterative updates are required within each time step until the maximum value of the node temperature difference between two adjacent iterations satisfies the set convergence criterion, for example, less than... Kelvin. Once the solution converges, the current temperature distribution at each node in the current time step is obtained.

[0067] Once the current temperature distribution is determined, grain evolution calculations within the same time step are immediately initiated. The grain growth kinetics equation describes the migration of grain interfaces and the increase in grain size in metallic materials at high temperatures. The grain interface refers to the boundary between two adjacent micro-regions with different crystal orientations. Under the influence of driving forces, the interface moves, causing some grains to engulf adjacent grains and grow. Grain growth kinetic parameters include the grain boundary migration activation energy and the rate constant, which are derived from high-temperature metallographic experiments conducted on the material of this metallic component. The specific method involves heating several cylindrical alloy structural steel samples with a size of 15 mm to 900, 950, 1000, 1050, and 1100 degrees Celsius, respectively, and holding them at that temperature for 0.5, 1, 2, 4, and 8 hours, respectively, followed by rapid quenching to fix the microstructure. Subsequently, the samples are ground and polished along the cross-section and then etched with a saturated picric acid solution to reveal the grain boundaries. The intercepts of at least 300 grains in each sample are measured using the intercept method and multiplied by a statistical coefficient of 1.5 to convert them into the average grain diameter, thereby establishing a database of the relationship between temperature, time, and grain size. The activation energy and rate constant for grain boundary migration are obtained by inverse fitting of this database.

[0068] The activation energy for grain boundary migration is The unit is kilojoules per mole; the rate constant is The unit is micrometers to the power of n per second; the grain growth index is These parameters are dimensionless. Together, they define the grain growth kinetics model, whose differential form expresses the relationship between the instantaneous grain growth rate and the current grain diameter and temperature. For each finite element mesh node, based on the current temperature of that node... Using grain growth kinetic parameters, calculate the grain boundary migration rate coefficient of this node at the current temperature. : ; In the formula, The grain boundary migration rate coefficient, whose dimensions depend on The dimensions of the quantity, for When the value is 2, the dimension is micrometers squared per second; The rate constant is obtained through nonlinear regression based on the above metallographic experimental data. The setting is based on the premise that the correlation coefficient between the model-predicted grain diameter after heat preservation and the measured value is greater than 0.95. The activation energy for grain boundary migration is expressed in kilojoules per mole, and is set to 385 kilojoules per mole. The universal gas constant has the following value: kilojoules per mole per Kelvin, this value multiplied by temperature We then obtained the energy term in kilojoules per mole; The current temperature of the node is expressed in Kelvin.

[0069] Grain boundary migration rate coefficient and the grain diameter of that node at the previous time step Substituting into the grain growth kinetics model, the instantaneous grain growth rate at this node is obtained. This represents the time step counter in the current finite element simulation of heat treatment. The grain growth kinetics model uses an exponential growth equation, and its instantaneous rate form is: ; In the formula, For the instantaneous grain growth rate, when and When the unit is micrometers squared per second, The unit is micrometer, and the unit of speed is micrometer per second; The grain diameter; The grain growth index, based on regression analysis of metallographic experimental data, is set at 2.0, which is consistent with the theoretical index for normal grain growth in ideal single-phase alloys. The physical meaning of this equation is that the driving force for grain growth is directly proportional to the interfacial energy and inversely proportional to the current grain size; the larger the grain, the slower the growth.

[0070] The grain diameter increment is obtained by numerically integrating the instantaneous grain growth rate within the current time step. A conventional explicit numerical integration scheme, such as the Euler forward integration scheme, is used here. This assumes the instantaneous rate remains constant within the step, and the grain diameter increment is calculated by multiplying the instantaneous grain growth rate by the time step. Subsequently, this grain diameter increment is added to the grain diameter from the previous time step to obtain the real-time average grain diameter of that node in the current time step. The same calculation process is performed on all finite element mesh nodes to obtain the real-time average grain diameter distribution of the entire metal component in the current time step. Thus, step S4 completes the update of the temperature field and grain size field within one time step, providing the necessary input for the stress analysis in step S5.

[0071] For example, continuing the process of alloy structural steel metal components processed in steps S1 to S3, the current time step is the 320th time step, the cumulative solution time is t = 160.0 seconds, the time step size is 0.5 seconds, and the quenching oil medium temperature is 298.15 Kelvin. Taking the thin-walled region node on the outer surface with node number 12087, the temperature of this node at the end of the previous time step, i.e., the 319th time step, was 638.2 Kelvin, and the grain diameter was 23.85 micrometers.

[0072] To determine the current temperature distribution, based on the function interpolation in step S3, the convective heat transfer coefficient of node 12087 at the current temperature is estimated to be 1803.5 watts per square meter per Kelvin, the material thermal conductivity is 31.2 watts per meter per Kelvin, and the specific heat capacity is 642 joules per kilogram per Kelvin. Substituting these data into the integrated heat capacity matrix and thermal stiffness matrix at the element level, the lumped heat capacity of the element containing this node is calculated as follows: The temperature is measured in Joules per Kelvin, which includes the equivalent volume represented by the node. Solved using an implicit scheme, the temperature of this node at the current time step is 631.5 Kelvin, or 358.35 degrees Celsius. The temperature solution for the entire model converges after three nonlinear iterations.

[0073] Entering the grain evolution calculation stage. The grain growth kinetic parameters for the metallic component material are: grain boundary migration activation energy of 385 kJ / mol, and rate constant of... The growth rate is micrometers squared per second, and the grain growth index is 2.0. The current temperature of node 12087 is 631.5 Kelvin. Substituting these parameters into the grain boundary migration rate coefficient formula, the calculated grain boundary migration rate coefficient is: Micrometers squared per second. This extremely low grain boundary migration rate coefficient indicates that grain growth has almost stopped at a low temperature of 358 degrees Celsius.

[0074] Substituting the grain boundary migration rate coefficient and the grain diameter from the previous time step into the instantaneous grain growth rate equation, the instantaneous rate is calculated as follows: Micrometers per second—this order of magnitude indicates that grain growth has physically terminated at this temperature. The grain diameter increment is... The value is in micrometers, which can be considered zero in double-precision floating-point representation. Therefore, the real-time average grain diameter at the current time step is updated to 23.85 micrometers, while the grain size remains at the level of the previous time step.

[0075] Compared to the 10th time step when the node temperature was still 1050 Kelvin, the calculated grain boundary migration rate coefficient was... The instantaneous grain growth rate is in micrometers squared per second. Micrometers per second, grain diameter increment is In this high-temperature range, within one time step, the grain diameter grows from 23.88 micrometers to 23.880000154 micrometers. Although each increment is small, the cumulative effect over hundreds of high-temperature time steps results in a considerable grain size change. This temperature-driven grain growth process reflects the physical nature of how, as the component cooling temperature decreases, the grain boundary migration rate coefficient decays exponentially, the grain growth rate drops sharply, and eventually freezes. The current temperature distribution and real-time average grain diameter distribution of all 156,243 nodes are calculated and stored according to the above steps, serving as the complete physical state passed to step S5 for the current time step.

[0076] S5. Substitute the current grain diameter distribution into the grain size-yield strength relationship to calculate the local real-time yield strength of each node. Then, combine thermal strain and elastoplastic constitutive model to perform stress balance iteration to obtain the stress tensor of the current time step.

[0077] In a preferred embodiment, the current grain diameter distribution is substituted into the grain size-yield strength relationship to calculate the local real-time yield strength of each node, including: obtaining the grain refinement coefficient and single-crystal friction stress of the metal component material; for each node, the real-time average grain diameter, the grain refinement coefficient, and the single-crystal friction stress are input into the relationship that the yield strength is equal to the single-crystal friction stress plus the ratio of the grain refinement coefficient to the square root of the grain diameter to calculate the local real-time yield strength; the calculated local real-time yield strength is used as the yield criterion parameter of the node at the current time step to update the yield criterion in its elastoplastic constitutive model.

[0078] In a further preferred embodiment, stress equilibrium iteration is performed by combining thermal strain and elastoplastic constitutive model to obtain the stress tensor of the current time step, including: calculating the thermal strain increment of each node using the temperature field difference between the current time step and the previous time step; embedding the local real-time yield strength into the yield function of the elastoplastic constitutive relationship, and forming the element stiffness matrix of each element by combining the thermal strain increment and the thermoelastic parameters of the material; assembling the global stiffness matrix and applying displacement boundary constraints for nonlinear iterative solution until the nodal force residuals meet the convergence condition to obtain the stress tensor of each node in the current time step.

[0079] Specifically, step S5 builds upon the calculation results of the current temperature distribution and real-time average grain diameter distribution completed in step S4 within the same time step, transforming the microscopic grain size information into macroscopic mechanical strength parameters, and on this basis, solving for the stress tensor generated inside the component due to uneven cooling in the current time step. This step achieves cross-physics coupling from thermal and microstructure analysis to mechanical analysis, and is further divided into two sequentially executed sub-stages: first, the local real-time yield strength of each node is calculated using the grain size-yield strength relationship; then, based on this yield strength, an elastoplastic finite element stress balance iterative solution considering thermal strain is performed.

[0080] First, we will explain the process of converting the real-time average grain diameter distribution into local real-time yield strength. Local real-time yield strength refers to the critical stress value at each finite element mesh node, under the current time step's microstructure, at which the material begins to undergo plastic deformation, measured in Pascals. This value is determined by the relationship between grain size and yield strength; the relationship used here is the classic grain size-yield strength equation. The grain size-yield strength equation is an empirical equation describing the linear positive correlation between the reciprocal of the square root of the average grain diameter and the yield strength in polycrystalline metallic materials. Its physical essence lies in the fact that grain interfaces act as barriers to dislocation slip; the finer the grains, the higher the interface density, and the stronger the resistance to dislocation movement, macroscopically manifesting as a higher yield strength.

[0081] Before performing this calculation, it is necessary to obtain the grain refinement coefficient and single-crystal friction stress of the metallic component material. The grain refinement coefficient is denoted as... is a material constant that measures the efficiency of grain interface strengthening, and its unit is the square root of megapascals (MPa) multiplied by micrometers (µm). Single-crystal frictional stress is denoted as , representing the inherent resistance of the lattice within the grain to dislocation motion, can be understood as the yield strength of a single crystal without grain boundary strengthening, measured in megapascals (MPa). These two parameters were calculated by back-calculating room-temperature tensile tests on alloy structural steel samples with different grain sizes obtained under various heat treatment regimes. Specifically, six groups of samples were prepared, and different microstructures with average grain diameters ranging from 8 to 80 micrometers were obtained by controlling the quenching temperature and holding time. Three standard tensile samples from each group were used to determine the yield strength. The yield strength was plotted as the ordinate and the reciprocal of the square root of the grain diameter. A linear regression was performed using the x-axis as the abscissa, and the slope of the regression line was the grain refinement strengthening coefficient. The longitudinal intercept is the single-crystal frictional stress. Based on the regression of this experiment, Values , The value is 245 MPa.

[0082] For each node where the real-time average grain diameter distribution has been generated in step S4, extract the real-time average grain diameter at the current time step. This value is in micrometers. It includes the obtained fine-grain strengthening coefficient. and single crystal friction stress Substituting the values ​​into the grain size-yield strength formula, the local real-time yield strength is calculated. : ; In the formula, This represents the local real-time yield strength, expressed in megapascals (MPa). The frictional stress of a single crystal is expressed in megapascals (MPa). This is the grain refinement strengthening factor, in units of... ; The real-time average grain diameter is expressed in micrometers. The unit is Therefore, the ratio The unit is megapascal (MPa), and both sides of the equation have the same dimensions. This formula clearly shows that the smaller the grain diameter and the larger the reciprocal of the square root, the higher the yield strength; as the grains gradually grow during heat treatment... As the yield strength increases, the yield strength decreases accordingly.

[0083] Please see Figure 3 As shown, Figure 3 This figure shows the coupling relationship between grain size and local real-time yield strength in this invention. As shown, the transformation from microstructure to macroscopic mechanical strength is achieved through the grain size-yield strength relationship. For the grain size distribution that evolves over time, the real-time average grain diameter is substituted into the strengthening model, ensuring that regions in the fine-grained state maintain high yield strength. As the temperature process progresses and grains grow, the local real-time yield strength decreases nonlinearly, ultimately achieving dynamic embedded correction of the macroscopic elastoplastic constitutive yield criterion by microstructure evolution.

[0084] Due to temperature variations, the yield strength of materials decreases further at high temperatures. This temperature effect is reflected by introducing a temperature correction factor. The correction factor is denoted as... Its form is: ; In the formula, This is a temperature correction factor, dimensionless, with a value ranging from zero to one. As the temperature increases, this factor decreases, thus reducing the effective yield strength. The current temperature of the node is in Kelvin and is taken from the current temperature distribution output in step S4. The reference temperature is set at 298.15 Kelvin, or 25 degrees Celsius, which is the baseline condition for the grain size-yield strength coefficient regression at this temperature. The melting point of this material is approximately 1773 Kelvin, or 1500 degrees Celsius. The temperature index is derived by fitting high-temperature tensile yield strength data of the material at 423 Kelvin, 573 Kelvin, and 723 Kelvin, and is set to 1.5. The effective local real-time yield strength after temperature correction is: ; The calculated effective local real-time yield strength The yield criterion in the elastoplastic constitutive model is updated using the yield criterion parameter of this node at the current time step. The elastoplastic constitutive model is a mathematical model describing the stress-strain relationship of a material under stress, where it first undergoes elastic deformation and then plastic yielding. The yield criterion specifies the stress state under which a material transitions from the elastic stage to the plastic stage. Here, the conventional Mises yield criterion is adopted, with the effective local real-time yield strength as the critical equivalent stress limit for plastic yielding. This criterion is parameterized by the updated effective local real-time yield strength of each finite element mesh node, determining the elastoplastic behavior of that node in each iteration step during subsequent stress solutions.

[0085] After updating the yield strength, the process proceeds to the stress equilibrium solution stage for the current time step. The stress tensor is a second-order symmetric tensor containing six independent components: normal stress and shear stress, describing the three-dimensional stress state of a material element. The goal of the stress equilibrium solution is to find a displacement field that satisfies force equilibrium and displacement continuity throughout the entire component, given the thermal strain caused by the temperature distribution and the aforementioned elastoplastic constitutive relations. The stress tensor is then derived from the displacement field using geometric and constitutive equations.

[0086] First, by utilizing the temperature field difference between the current time step and the previous time step, and combining this with the linear expansion coefficient of the material, the thermal strain increment at each node is calculated. Thermal strain is the volume expansion or contraction of a material due to temperature changes. Under unconstrained conditions, it does not generate stress. However, when the temperature distribution within a component is uneven, leading to inconsistent expansion and contraction, and when it is subject to geometric constraints and external boundary conditions, thermal strain becomes the direct source of thermal stress. Thermal expansion is considered isotropic, occurring only in the normal strain direction, and there is no thermal shear strain.

[0087] The corrected effective local real-time yield strength will then be... The element stiffness matrix of each element is formed by embedding the yield function into the elastoplastic constitutive relation, combining the thermal strain increment and the thermoelastic parameters of the material. The thermoelastic parameters of the material include the elastic modulus E and Poisson's ratio. The data on the elastic modulus versus temperature were obtained from high-temperature ultrasonic resonance testing of the material, with a Poisson's ratio of 0.3. The plastic strain increment was determined by the plastic flow law and constrained by the plastic consistency condition. A conventional elastoplastic integral algorithm (e.g., radial return algorithm) was used to perform nonlinear solutions at the integration point level, updating the stress and outputting a consistent tangent stiffness matrix. The element stiffness matrix was then obtained by assembling the elements using the standard formula for isoparametric elements.

[0088] Using a conventional finite element assembly method, the element stiffness matrices of all elements are accumulated according to the nodal degrees of freedom to assemble the global stiffness matrix, and the external force load vector is assembled simultaneously. When applying displacement boundary constraints, simply supported constraints are applied at three points on the central axis of symmetry of the bottom surface of the metal component, that is, restricting its vertical displacement along the axial direction and two tangential displacements, so as to avoid rigid body displacement while constraining the thermal deformation degrees of freedom of the component.

[0089] Nonlinear iterative solutions involve repeatedly performing stress updates and stiffness reassembly iterations at each time step due to the path dependence and nonlinearity of the constitutive relations in elastoplastic constitutive models, until the nodal force residuals satisfy the convergence condition. Conventional iterative algorithms for solving nonlinear equations, such as the Newton-Raphson iterative method, have a convergence condition where the ratio of the norm of the nodal force residuals to the norm of the external loads is less than... Furthermore, the ratio of the norm of the displacement increment to the norm of the current displacement is less than [the value of the norm]. Typically, convergence requires 3 to 6 iterations within a single cooling time step. After convergence, the complete stress tensor for each node at the current time step is obtained, including... The stress tensor consists of six components, measured in megapascals (MPa). This stress tensor reflects the transient thermal stress state at the current moment, caused by temperature inhomogeneity and grain size gradient. As the temperature cools to room temperature, this transient stress will eventually evolve into residual stress permanently retained within the component.

[0090] For example, continuing the calculation for time step 320 in step S4, the outer surface node of the thin-walled region with node number 12087 is retrieved. Its current temperature is 631.5 Kelvin, the temperature in the previous time step was 638.2 Kelvin, and the real-time average grain diameter is 23.85 micrometers. First, the yield strength is calculated. The fine-grain strengthening coefficient... Single crystal friction stress Megapascals (MPa). Substituting the above parameters into the grain size-yield strength relationship, the local real-time yield strength is calculated. The value is 248.79 MPa. Next, the temperature correction factor is calculated, taking the current temperature and the reference temperature as factors. Kelvin, melting point Substituting Kelvin and the temperature index of 1.5 into the formula, the temperature correction factor is calculated. The value is 0.8926, and then the effective local real-time yield strength is obtained. The value is 222.07 MPa. This value is updated to the yield criterion corresponding to node 12087, serving as the plastic initiation condition for that node at the current time step.

[0091] Then, the stress equilibrium solution is obtained. The temperature change at node 12087 is considered. Kelvin, a negative value indicates cooling, resulting in contraction and thermal strain. Coefficient of linear expansion. For each Kelvin, combined with the coefficient of linear expansion, the positive components of the thermal strain increment are calculated as follows: The tangential thermal strain component is zero. This cooling contraction strain, under the condition that the cooling rates of the thin-walled region and the thick-walled flange region are different, further develops the differential stress between the tensile stress in the thin-walled region and the compressive stress in the flange region.

[0092] At this time step, the elastic modulus of the element containing node 12087 near 400 degrees Celsius is E = 187 gigapascals. Megapascals, Poisson's ratio Combining The megapascals and thermal strain increments are plastically corrected at the integration point level using a conventional elastoplastic integration algorithm. After initializing the elastic test stress bias, the Mises yield criterion is checked. When the equivalent elastic test stress exceeds 222.07 MPa, plastic correction is performed to pull the stress state back to the yield surface. After constitutive integration, the stress tensor components and uniform tangent modulus at the integration point are output. After integration at all 8 Gaussian points within the element, a 24×24 element stiffness matrix is ​​formed, with the element value in the first row and first column being... The stiffness matrix of all 98,520 elements was assembled into a global stiffness matrix containing approximately 460,000 degrees of freedom. After applying displacement boundary constraints, a conventional large sparse matrix solver (such as the preprocessed conjugate gradient method) was used to solve the matrix. The solution converged after four Newton iterations. The ratio of the nodal force residual norm to the load norm was... After convergence, the stress tensor at node 12087 is... Megapascal, Megapascal, Megapascal, Megapascal, Megapascal, The stress is measured in megapascals (MPa). The normal stress component in the x-direction is close to but does not exceed the yield limit of 222.07 MPa, indicating that the node is in an elastic unloading state at this time step. The stress caused by cooling contraction is insufficient to trigger further plastic flow. This stress will gradually accumulate as a residual value, eventually becoming part of the residual stress remaining at room temperature. After the stress tensors of all nodes are calculated, the mechanical state update for this time step is complete, and the process awaits step S6 to determine whether to continue the cycle.

[0093] S6. Determine whether the cumulative solution time has reached the total heat treatment time. If not, use the current temperature distribution, real-time average grain diameter distribution, and stress tensor of the current time step as the initial conditions for the next time step, and return to the time step increment loop. If it has reached the time step, output the three-dimensional residual stress distribution and final grain size distribution of the metal component.

[0094] In a preferred embodiment, it is determined whether the cumulative solution time has reached the total heat treatment time. If not, the current temperature distribution, real-time average grain diameter distribution, and stress tensor of the current time step are used as the initial conditions for the next time step, and the process returns to the time step increment loop. If the time has reached the total heat treatment time, the three-dimensional residual stress distribution and final grain size distribution of the metal component are output, including: obtaining the total heat treatment time specified by the heat treatment process of the metal component as the total heat treatment time; after the calculation of each time step is completed, the current cumulative solution time is compared with the total heat treatment time; in response to the cumulative solution time being less than the total heat treatment time, the current temperature distribution is assigned to the temperature distribution of the previous time step, the real-time average grain diameter distribution is assigned to the grain diameter distribution of the previous time step, the stress tensor is assigned to the stress field of the previous time step, the time step is updated, and the calculation of the next time step is returned; in response to the cumulative solution time being equal to or greater than the total heat treatment time, the loop is terminated, the final obtained stress tensor is used as the residual stress, the final obtained real-time average grain diameter distribution is used as the final grain size distribution, and the corresponding three-dimensional visualization cloud map is generated and output.

[0095] In a further preferred embodiment, the time-step increment loop includes the following coupled calculation process: within each time step, the currently solved temperature distribution is input into the grain growth kinetic equation to calculate and update the real-time average grain diameter distribution; the updated real-time average grain diameter distribution is input into the grain size-yield strength relationship to calculate and update the local real-time yield strength; the updated local real-time yield strength is used as the yield criterion of the elastoplastic constitutive model to calculate the plastic strain increment and stress tensor of the current time step; when the cumulative solution time of the time-step increment loop reaches the total heat treatment time, the final grain size distribution and three-dimensional residual stress distribution are output synchronously.

[0096] Specifically, step S6 is the control center of the time-step increment loop. It is responsible for determining time progression and passing state variables generation by generation after each time step calculation, and outputting the final three-dimensional residual stress distribution and final grain size distribution at the end of the entire heat treatment simulation. The time-step increment loop is a progressive solution framework that runs through steps S3 to S5, and step S6 provides the logical gate for continuing or terminating the loop. The total heat treatment duration, cumulative solution time, and the variable substitution relationships of each physical field between time steps constitute the core operational elements of this step.

[0097] The total heat treatment time refers to the duration specified in the heat treatment process document for the metal component, from the start of quenching to the complete cooling of the component and its removal from the bath after reaching equilibrium with the medium temperature. For the alloy structural steel metal component with blind holes and flange structures in this embodiment, the process specification is as follows: austenitization at 850 degrees Celsius, followed by immersion quenching in rapid quenching oil at 25 degrees Celsius, and holding until the component surface temperature drops below the medium temperature. The process specifies a total duration of 3600 seconds from the moment of medium immersion. This value is determined by the process engineer based on the maximum wall thickness of the component, the hardenability of the material, and the characteristics of the cooling medium to ensure complete transformation of the core structure and the attainment of thermal equilibrium. It is obtained by parsing the process specification parameters and writing them into the solver input file.

[0098] After the calculation of each time step is completed, i.e., after the stress tensor of the current time step is obtained in step S5, step S6 compares the current cumulative solution time with the total heat treatment time. The current cumulative solution time is the total time value accumulated from the start of the heat treatment to the end of the current time step, in seconds. It is updated by adding the step size of the current time step to the cumulative solution time at the end of the previous time step at the end of each time step.

[0099] The comparison operation is implemented through a conditional judgment logic. When the cumulative solution time is less than the total heat treatment duration, it is determined that the heat treatment process has not yet ended, and all the physical fields just calculated in the current time step need to be used as the starting state for the solution of the next time step. This variable substitution operation involves the complete transfer of three physical fields. The first is to assign the current temperature distribution to the temperature distribution of the previous time step, that is, to completely overwrite the array storing the temperature field of the previous time step with the temperature values ​​of each node obtained by solving the heat conduction equation in step S4 for the current time step.

[0100] The second step is to assign the real-time average grain diameter distribution to the grain diameter distribution of the previous time step. This operation takes the real-time average grain diameter value of the current time step obtained by numerical integration of the grain growth kinetic equation in step S4 and passes it to the array storing the grain diameter of the previous time step.

[0101] The third step is to assign the stress tensor to the stress field of the previous time step. The stress tensor contains six independent components. This operation will gradually cover the complete symmetric stress tensor of the corresponding node in the previous time step with the stress components of each node in the current time step output by the finite element stress equilibrium iterative solution in step S5.

[0102] After completing the node-by-node assignment of the three physical fields mentioned above, the time step is updated, that is, the step size of the next time step is determined according to the current temperature state, and the time step count is incremented by one. Then, the process returns to step S3 to re-establish the node heat transfer boundary conditions corresponding to the new moment and executes all calculations in steps S4 and S5. This behavior constitutes a single complete iteration of the time step increment loop.

[0103] Conversely, when the cumulative solution time is equal to or greater than the total heat treatment time, the heat treatment process is considered to have been fully executed, and the loop terminates. The stress tensor at this time step is considered the final stress state remaining inside the component under room temperature conditions, i.e., the three-dimensional residual stress distribution. The three-dimensional residual stress distribution refers to the stress field at each point inside the metal component after all heat treatment processes are completed and it is not subjected to any external mechanical loads. Its unit is megapascals (MPa), and it consists of six complete stress components. It can be scalarized and characterized using the conventional Mises equivalent stress. The equivalent stress calculation formula is: The unit is megapascal (MPa). Simultaneously, the real-time average grain diameter distribution at this time step is considered the final grain size distribution that no longer changes after the component cools to room temperature, and the unit is micrometer (µm).

[0104] Generating and outputting the corresponding 3D visualization cloud map is a graphical rendering process that transforms discrete numerical data distributed across all finite element mesh nodes into an intuitive continuous color chart. The specific process involves reading the nodal stress tensor file and grain diameter file output from the final time step, extracting the equivalent stress value and grain diameter value of each node as a scalar field, and using conventional finite element post-processing rendering techniques to map the values ​​within the element to points within the element. In the user-defined section view, the numerical ranges are mapped to a color space according to a color scale mapping table. The equivalent stress gradually transitions from blue (representing zero stress) to red (representing the highest stress region), and the grain diameter transitions from light color (representing fine-grained regions) to dark color (representing coarse-grained regions). Finally, the result is rendered as a 2D cross-sectional cloud map or a 3D semi-transparent volume rendering cloud map, and output as image files and data reports.

[0105] In the complete time-step incremental loop, each round of advancement from S3 to S6 follows a tightly coupled real-time embedded iterative mechanism. Whenever step S4 solves the heat conduction equation to obtain the current temperature distribution, the change in temperature distribution immediately updates the real-time average grain diameter distribution directly through the grain growth kinetics equation. This update does not wait for the entire time step to complete but is performed continuously within the same step. Furthermore, the change in the real-time average grain diameter distribution directly updates the local real-time yield strength through the grain size-yield strength relationship. The change in the local real-time yield strength directly affects the plastic strain increment and stress tensor of the current time step in step S5. This forms a real-time embedded iterative coupling between the temperature history-driven microscopic grain evolution and the macroscopic mechanical response. The characteristic of this coupling is that the final grain size distribution and the three-dimensional residual stress distribution are not obtained independently step-by-step, but are updated sequentially within each time step through the same time-step incremental loop and ultimately reach the heat treatment endpoint synchronously. This ensures that the microstructure evolution history and the macroscopic stress evolution history are consistent on the time axis, eliminating the errors introduced by time point mismatch in traditional step-by-step decoupled calculations.

[0106] For example, following the complete processing of the alloy structural steel metal component at time step 320 in steps S1 to S5, the execution logic of step S6 is substituted. The total duration of the heat treatment process for this metal component is 3600 seconds. After time step 320, the current cumulative solution time is 160.0 seconds, which is obtained by adding the cumulative time of 159.5 seconds at the end of time step 319 to the current step size of 0.5 seconds. Comparing this cumulative solution time with the total heat treatment duration of 3600 seconds, it is clear that 160.0 < 3600, the condition is met, and it is determined that the heat treatment process has not yet ended.

[0107] Therefore, variable substitution is performed: the current temperature distribution is assigned to the temperature distribution of the previous time step, the real-time average grain diameter distribution of the current time step is assigned to the grain diameter distribution of the previous time step, and each component of the stress tensor of the current time step is assigned to the stress field of the previous time step, completing the assignment operation for all nodes one by one. Then, the time step size is updated. The highest temperature of the current node is approximately 634 Kelvin, still in the initial rapid cooling phase. The step size of the next time step is determined to be 0.5 seconds, the time step count m is updated to 321, and the process returns to step S3 to re-establish the node heat transfer boundary conditions and continue the simulation of the cooling process.

[0108] The loop continues in this pattern. When the cumulative solution time reaches 3600 seconds, the time step count is approximately 9400 steps, and the temperature of node 12087 has dropped to 299.3 Kelvin (26.15 degrees Celsius), approaching the quenching oil medium temperature of 298.15 Kelvin (25 degrees Celsius). At this point, the real-time average grain diameter of node 12087 is 23.85 micrometers, consistent with the size after freezing at the end of the high-temperature stage. Its stress tensor has six components as follows: Megapascal, Megapascal, Megapascal, Megapascal, Megapascal, The equivalent stress was calculated to be 193.6 MPa. This residual tensile stress value is consistent with the expected distribution of the external surface tensile residual stress caused by the thin-walled region's initial cooling and contraction, and its constraint by the thick-walled flange region.

[0109] Please see Figure 4 As shown, Figure 4 This is a histogram of the residual stress tensor components at the end of the heat treatment process according to the present invention. As shown in the figure, the quantized output of the three-dimensional residual stress consists of six independent tensor components obtained from the elastoplastic constitutive integral. For the thin-walled region of complex components, the incremental loop is terminated when the accumulated solution time reaches the total heat treatment time. The three normal stress components and three shear stress components of the spatial node are directly extracted and comprehensively converted into equivalent stress, intuitively presenting the final three-dimensional stress state of the component under the combined action of uneven cooling and geometric constraints.

[0110] The cumulative solution time equals the total heat treatment time of 3600 seconds, the condition criterion is met, and the loop terminates. The total stress tensor at the final moment constitutes the three-dimensional residual stress distribution, and the real-time average grain diameter distribution at the final moment constitutes the final grain size distribution; both are simultaneously written into the results file. In the post-processing stage, this results file is read, the equivalent stress field is extracted, and a color contour map is drawn in the longitudinal cross-sectional plane passing through the blind hole axis. Blue corresponds to low-stress areas with equivalent stress below 50 MPa, and red corresponds to high-stress areas with equivalent stress above 300 MPa. The grain size field is extracted, and the grain size variation from 20 micrometers to 45 micrometers is represented by a grayscale gradient within the same cross-sectional view. The contour map visually shows stress concentration in the transition fillet area at the flange root, a lower residual stress level in the bottom area of ​​the blind hole, and a circumferential tensile stress band at the edge of the hole. The grain size exhibits a gradient distribution between the fine-grained region on the outer surface and the relatively coarsened region in the core due to slow cooling. Finally, the grain size distribution and the three-dimensional residual stress distribution are synchronously output at the end of 3600 seconds through the same time step increment cycle, achieving the coupling of the microstructure evolution history and the macro-stress development history.

[0111] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A method for predicting residual stress after heat treatment based on microstructure evolution, characterized in that, include: S1. Obtain the three-dimensional geometric model of the metal component with blind holes and flanges and divide it into finite element meshes. Map the obtained initial temperature field, initial grain diameter distribution and initial stress and strain field of the component to the finite element mesh nodes as the initial physical field. S2. Based on the finite element mesh nodes, calculate the local characteristic thickness of each node on the outer surface and the local depth-to-diameter ratio of each node on the inner surface of the blind hole. Divide the outer surface nodes into thin-walled nodes and thick-walled flange nodes according to the local characteristic thickness, and mark the nodes on the inner surface of the blind hole. S3. Based on the node division results, apply the thin-walled and thick-walled convective heat transfer coefficients determined by the corresponding heat transfer characteristic functions to the thin-walled and thick-walled flange nodes respectively. For the blind hole inner cavity surface nodes, the reference heat transfer coefficient is corrected by the attenuation function of the local depth-to-diameter ratio to form the heat transfer boundary. S4. In the time step increment loop, the heat conduction equation is solved using the heat transfer boundary. Then, based on the current temperature and the grain diameter of the previous time step, the grain growth kinetic equation is called to calculate the grain growth rate and perform time integration to generate the real-time average grain diameter distribution of the current time step. S5. Substitute the current grain diameter distribution into the grain size-yield strength relationship, calculate the local real-time yield strength of each node, and then combine thermal strain and elastoplastic constitutive stress balance iteration to obtain the stress tensor of the current time step. S6. Determine whether the cumulative solution time has reached the total heat treatment time. If not, use the current temperature distribution, real-time average grain diameter distribution, and stress tensor of the current time step as the initial conditions for the next time step, and return to the time step increment loop. If it has reached the time step, output the three-dimensional residual stress distribution and final grain size distribution of the metal component.

2. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The method involves calculating the local characteristic thickness of each node on the outer surface and the local depth-to-diameter ratio of each node on the inner surface of the blind hole, based on the local characteristic thickness. The outer surface nodes are then divided into thin-walled and thick-walled flange-area nodes, and the nodes on the inner surface of the blind hole are marked. Extract the finite element mesh nodes of the outer surface of the 3D geometric model and obtain the surface unit normal vector at each outer surface node; Starting from each outer surface node, draw a ray into the component along its surface normal direction, find the intersection point of the ray with the opposite outer surface of the component, calculate the spatial distance between the starting point and the intersection point, and use the spatial distance as the local feature thickness of the outer surface node. Extract the finite element mesh nodes of the blind hole cavity, calculate the projection distance of each node on the surface of the blind hole cavity along the blind hole axis to the blind hole opening plane as the node depth, obtain the blind hole diameter size, and take the ratio of the node depth to the blind hole diameter size as the local depth-to-diameter ratio of the node. Statistical clustering is performed on the local feature thickness of all calculated outer surface nodes to identify two thickness clustering centers that represent thin-walled structures and thick-walled flange structures, respectively. The classification thresholds for thin-walled and thick-walled flange regions are determined based on the thickness clustering centers. Nodes with local feature thickness less than the classification threshold are classified as thin-walled nodes, and nodes with local feature thickness not less than the classification threshold are classified as thick-walled flange nodes.

3. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The application of thin-walled and thick-walled convective heat transfer coefficients, determined by the corresponding heat transfer characteristic functions, to the thin-walled nodes and thick-walled flange nodes respectively includes: The measured temperature change data of thin-walled geometric feature samples during the quenching and cooling process are obtained. The thin-walled heat transfer characteristic function of the heat transfer coefficient as a function of surface temperature is obtained by solving the inverse problem of heat conduction. The measured temperature change data of thick-walled geometric feature samples during the quenching and cooling process are obtained. The thick-walled heat transfer characteristic function of the heat transfer coefficient as a function of surface temperature is obtained by solving the inverse problem of heat conduction. For each thin-walled node, the surface temperature of that node at the current time step is obtained, and the thin-walled convective heat transfer coefficient is calculated using the thin-walled heat transfer characteristic function; for each thick-walled flange node, the surface temperature of that node at the current time step is obtained, and the thick-walled convective heat transfer coefficient is calculated using the thick-walled heat transfer characteristic function.

4. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The correction of the reference heat transfer coefficient for the surface nodes of the blind hole cavity using a local depth-to-diameter ratio attenuation function includes: Obtain the convective heat transfer coefficient at the node at the blind hole opening and determine it as the reference heat transfer coefficient; Construct a monotonically decreasing exponential decay factor function with the local depth-to-diameter ratio as the independent variable, such that the decay factor function takes a value of one when the local depth-to-diameter ratio is zero, and gradually decays as the local depth-to-diameter ratio increases. For each node on the surface of the blind hole cavity, the corresponding local depth-to-diameter ratio is input into the attenuation factor function to obtain the attenuation factor. The product of the reference heat transfer coefficient and the attenuation factor is calculated as the convective heat transfer coefficient of the blind hole cavity of that node.

5. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The step of calling the grain growth kinetic equation based on the current temperature and the grain diameter of the previous time step, calculating the grain growth rate, performing time integration, and generating the real-time average grain diameter distribution for the current time step includes: The grain growth kinetic parameters of the metallic component material are obtained, including the grain boundary migration activation energy and the rate constant. For each finite element mesh node, the grain boundary migration rate coefficient at the current temperature is calculated based on the node's current temperature and grain growth kinetic parameters. Substituting the grain boundary migration rate coefficient and the grain diameter of the node in the previous time step into the grain growth kinetics model, the instantaneous grain growth rate of the node is obtained. The instantaneous grain growth rate is numerically integrated within the current time step to obtain the grain diameter increment. This grain diameter increment is then added to the grain diameter of the previous time step to obtain the real-time average grain diameter of the node at the current time step.

6. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The step of substituting the current grain diameter distribution into the grain size-yield strength relationship to calculate the local real-time yield strength of each node includes: To obtain the fine-grained strengthening coefficient and single-crystal frictional stress of metallic component materials; For each node, the real-time average grain diameter, the grain refinement coefficient, and the single-crystal frictional stress are input into the formula that yield strength is equal to the single-crystal frictional stress plus the ratio of the grain refinement coefficient to the square root of the grain diameter, and the local real-time yield strength is calculated. The calculated local real-time yield strength is used as the yield criterion parameter of the node at the current time step, and the yield criterion in its elastoplastic constitutive model is updated.

7. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The stress equilibrium iteration, combining thermal strain and elastoplastic constitutive model, yields the stress tensor for the current time step, including: The thermal strain increment at each node is calculated by using the temperature field difference between the current time step and the previous time step. The local real-time yield strength is embedded into the yield function of the elastoplastic constitutive relation, and combined with the thermal strain increment and the thermoelastic parameters of the material, to form the element stiffness matrix of each element. Assemble the global stiffness matrix and apply displacement boundary constraints for nonlinear iterative solution until the nodal force residuals meet the convergence condition, and obtain the stress tensor of each node at the current time step.

8. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The determination is made as to whether the cumulative solution time has reached the total heat treatment time; if not, the current temperature distribution, real-time average grain diameter distribution and stress tensor of the current time step are used as the initial conditions for the next time step, and the cycle returns to the time step increment. If achieved, the output will show the three-dimensional residual stress distribution and final grain size distribution of the metal component, including: The total heat treatment time specified in the heat treatment process for the metal component is taken as the total heat treatment time. After each time step is completed, the current cumulative solution time is compared with the total heat treatment time. In response to the cumulative solution time being less than the total heat treatment time, the current temperature distribution is assigned to the temperature distribution of the previous time step, the real-time average grain diameter distribution is assigned to the grain diameter distribution of the previous time step, the stress tensor is assigned to the stress field of the previous time step, the time step size is updated, and the calculation of the next time step is returned. In response to the cumulative solution time being equal to or greater than the total heat treatment time, the loop is terminated, the final stress tensor is taken as the residual stress, the final real-time average grain diameter distribution is taken as the final grain size distribution, and the corresponding three-dimensional visualization cloud map is generated and output.

9. The method for predicting residual stress after heat treatment based on microstructure evolution according to claim 1, characterized in that, The time step increment loop includes the following coupled calculation processes: At each time step, the temperature distribution obtained from the current solution is input into the grain growth kinetics equation to calculate and update the real-time average grain diameter distribution; Input the updated real-time average grain diameter distribution into the grain size-yield strength relationship to calculate and update the local real-time yield strength. The updated local real-time yield strength is used as the yield criterion for the elastoplastic constitutive model to calculate the plastic strain increment and stress tensor at the current time step. When the cumulative solution time of the time step increment cycle reaches the total heat treatment time, the final grain size distribution and three-dimensional residual stress distribution are output synchronously.

Images