Method for optimizing comprehensive performance of chisel tooth heat sink and chisel tooth heat sink

Abstract

The application relates to a method for optimizing the comprehensive performance of a notched fin and the notched fin, and the method comprises the following steps: according to a preliminary stress amplification distribution map, a particle swarm optimization algorithm is used to determine a uniform pressure application path, and an optimized pressure distribution scheme is obtained; key node coordinate data is obtained from the optimized pressure distribution scheme, a crack initiation probability is predicted through a Monte Carlo simulation method, and a high-risk area list is generated; for the high-risk area list, a response surface method is used to construct a target response surface, and refined strengthening parameters are obtained; according to the refined strengthening parameters, the durability improvement effect is simulated, and an updated pressure distribution scheme is obtained; a final uniform pressure application path is obtained from the updated pressure distribution scheme, and the overall heat dissipation index optimization path is determined by combining the Monte Carlo simulation result; and through the overall heat dissipation index optimization path, fatigue damage evolution is simulated, and a complete strengthening treatment protocol is output.

Description

Technical Field

[0001] This invention relates to the field of heat sinks, and more particularly to a method for optimizing the overall performance of toothed heat sinks. Background Technology

[0002] In modern industry, heat sinks are core components that ensure stable equipment operation, and their performance directly affects the safety and reliability of critical systems such as electronic equipment, aerospace equipment, and communication base stations. In particular, serrated heat sinks, due to their efficient heat dissipation capabilities and compact structural design, have become the preferred choice for many high-load applications.

[0003] However, despite their wide application, the long-term service life of heat sinks is often limited, especially when facing the challenges of vibration and thermal stress under complex operating conditions. Current solutions for improving heat sink performance mostly focus on material upgrades or overall structural optimization, but these methods often overlook the specific needs of critical local areas. Many traditional strengthening methods are inadequate when dealing with the narrow and vulnerable area at the fin root of the heat sink. They either cannot be applied precisely due to limited operating space, or the thin-walled structure deforms due to process characteristics, potentially introducing new damage risks. This lack of adaptability to local areas makes improving the durability of heat sinks a persistent bottleneck. A deeper challenge lies in the complex stress environment experienced by the fin root. As the transition area connecting the substrate and the fins, the fin root experiences significant tensile force during processing due to tool compression. This tensile force weakens the material strength, making it highly susceptible to cracking under repeated stress or temperature changes. More importantly, the drastic shape changes in this region often amplify local stresses, and existing technologies struggle to effectively apply uniform and controllable pressure within such a confined space to mitigate this issue. This unfavorable stress distribution, coupled with limitations in spatial operation, makes the durability of the fin root a critical technological hurdle that urgently needs to be overcome.

[0004] Therefore, how to implement precise and effective reinforcement treatment on the root of the fins in a narrow space, while improving its stress state to resist fatigue damage, has become a key issue in improving the overall lifespan of the spade-tooth heat sink. Summary of the Invention

[0005] This invention provides a method for optimizing the overall performance of a spade-shaped heatsink, mainly including:

[0006] Geometric model data and historical stress distribution records of the root of the spade-shaped heat sink fins are obtained. Finite element analysis is used to simulate the weakening state of material strength under vibration load and thermal stress field, resulting in a preliminary amplified stress distribution map. Based on this preliminary amplified stress distribution map, a particle swarm optimization algorithm is used to determine the uniform pressure application path, obtaining an optimized pressure distribution scheme. Key node coordinate data are obtained from the optimized pressure distribution scheme, and the crack initiation probability is predicted using Monte Carlo simulation to generate a list of high-risk areas. For the high-risk area list, a target response surface methodology is used to construct a target response surface, obtaining refining and strengthening parameters. Based on these refining and strengthening parameters, the durability improvement effect is simulated, resulting in an updated pressure distribution scheme. The final uniform pressure application path is obtained from the updated pressure distribution scheme, and combined with the Monte Carlo simulation results, an overall heat dissipation index optimization path is determined. Using this overall heat dissipation index optimization path, fatigue damage evolution is simulated, and a complete strengthening treatment protocol is output. Furthermore, the step of acquiring the geometric model data and historical stress distribution records of the root of the spade-shaped heat sink fins, and simulating the material strength weakening state under vibration load and thermal stress field using the finite element analysis method to obtain a preliminary stress magnification distribution map, includes: extracting the root contour from the geometric model data; performing mapping processing on the geometric parameters and material properties of the root contour with the historical stress distribution records to obtain an initial state matrix; constructing a finite element mesh model based on the initial state matrix; obtaining the finite element mesh elements contained in the finite element mesh model; applying vibration load and thermal stress field to the finite element mesh model; and using the finite element method to determine the finite element mesh elements under the combined effects of vibration load and thermal stress field. The nodal displacement during application; the strain gradient of the finite element mesh element is calculated using the nodal displacement, and if the strain gradient is greater than a preset strain gradient threshold, the finite element mesh element is determined to have fatigue damage; for the finite element mesh element that has fatigue damage, a stiffness degradation coefficient is obtained to characterize stiffness degradation, the strength reduction of the finite element mesh element is calculated based on the stiffness degradation coefficient, and the material properties in the initial state matrix are updated based on the strength reduction to obtain a degradation state matrix; the updated stress distribution of the root profile is recalculated using the degradation state matrix, and the ratio of the updated stress distribution to the historical stress distribution record is calculated to obtain a preliminary stress magnification distribution map; further, the step of determining the uniform pressure application path using a particle swarm optimization algorithm based on the preliminary stress magnification distribution map to obtain an optimized pressure distribution scheme includes: obtaining the stress concentration region in the target space from the preliminary stress magnification distribution map, extracting tensile features based on the stress direction and magnitude of the stress concentration region, and calculating the spatial curvature of the target space based on the deformation data corresponding to the tensile features;By matching the spatial curvature with the preset surface curvature of the pressure-applying component, the contact area between the target pressure-applying component and the root of the spade-shaped heat sink fin is determined. An initial pressure angle is set based on the contact area and the force direction of the tensile characteristic. The load gradient is calculated using the initial pressure angle and a preset initial pressure value. At least one discrete pressure node is generated within the contact area based on the load gradient. A pressure path is constructed by connecting the pressure nodes. Using the pressure path as particles and the path coverage rate as the fitness function, a particle swarm optimization algorithm is used to optimize the pressure path, obtaining a path coverage set. The path coverage rate corresponding to each path in the path coverage set is calculated. If the path coverage rate is lower than a preset coverage rate threshold, the initial pressure angle is iteratively adjusted. The load gradient is recalculated based on the iteratively adjusted pressure angle, and particle optimization is performed again to obtain an optimized pressure distribution scheme. Furthermore, the step of obtaining key node coordinate data from the optimized pressure distribution scheme, predicting crack initiation probability using Monte Carlo simulation, and generating a high-risk area list includes: obtaining fatigue damage indicators corresponding to the key node coordinate data; extracting stress alternation amplitude and corresponding load cycle number based on the fatigue damage indicators; constructing a stress probability distribution function using the stress alternation amplitude and the load cycle number; generating random sampling stress samples using the stress probability distribution function; obtaining a preset material fracture toughness; calculating the initial crack length based on the random sampling stress samples and the material fracture toughness; calculating the crack propagation rate based on the initial crack length and the random sampling stress samples; obtaining a preset number of simulation iterations; using Monte Carlo simulation, iteratively accumulating the crack propagation rate based on the number of simulation iterations to obtain the crack initiation probability corresponding to the key node coordinate data; determining if the crack initiation probability is greater than a preset probability threshold, then extracting the corresponding key node coordinate data and generating a high-risk area list based on the key node coordinate data. Furthermore, for the list of high-risk areas, the method of constructing a target response surface using response surface methodology to obtain refined enhancement parameters includes: obtaining the region nodes corresponding to the list of high-risk areas; extracting geometric features and load boundaries based on the region nodes; obtaining preset initial enhancement parameters; constructing an initial response surface using response surface methodology based on the geometric features, the load boundaries, and the initial enhancement parameters; obtaining a preset verification sample set; calculating the predicted response value corresponding to the verification sample set using the initial response surface; calculating the fitting error based on the predicted response value and the actual response value of the verification sample set; determining if the fitting error is greater than a preset error threshold; then obtaining historical enhancement samples; and correcting the initial response surface based on the historical enhancement samples to obtain the target response surface.The parameter sensitivity matrix and response gradient are extracted from the target response surface. The parameter compensation vector is calculated using the parameter sensitivity matrix and the response gradient. The iteration step size is determined based on the parameter compensation vector. The initial enhancement parameters are updated based on the iteration step size and the parameter compensation vector to obtain the refined enhancement parameters. Furthermore, the step of simulating durability improvement effects and obtaining an updated pressure distribution scheme based on the refining and strengthening parameters includes: constructing a fatigue damage assessment model based on the refining and strengthening parameters and a preset load history sequence; calculating the cumulative fatigue damage value through the fatigue damage assessment model; extracting stress concentration areas based on the cumulative fatigue damage value; obtaining preset benchmark test indicators; extracting stress amplification and mitigation indicators for the stress concentration areas; calculating an improvement ratio based on the stress amplification and mitigation indicators and the benchmark test indicators; determining that if the improvement ratio is less than a preset improvement threshold, triggering a population backtracking mechanism; obtaining an initial population corresponding to a preset particle swarm optimization algorithm through the population backtracking mechanism; extracting individual extreme values ​​and global extreme values ​​for the initial population; calculating the particle update rate based on the individual extreme values ​​and the global extreme values; performing position mutation processing on the initial population based on the particle update rate to obtain a mutated particle swarm; extracting structural node displacement vectors through finite element simulation using the mutated particle swarm; reconstructing surface contact forces based on the structural node displacement vectors; and obtaining an updated pressure distribution scheme based on the surface contact forces. Furthermore, the step of obtaining the final uniform pressure application path from the updated pressure distribution scheme and determining the overall heat dissipation index optimization path based on Monte Carlo simulation results includes: extracting surface pressure nodes for the updated pressure distribution scheme; extracting pressure values ​​of adjacent nodes based on the surface pressure nodes; calculating local pressure differences based on the pressure values ​​of adjacent nodes; constructing the final uniform pressure application path based on the local pressure differences; extracting contact thermal resistance parameters for the final uniform pressure application path; calculating heat flux density and temperature gradient based on the contact thermal resistance parameters; constructing a joint probability distribution function based on the heat flux density and temperature gradient; and performing Monte Carlo simulation on the joint probability distribution function. The probability distribution function performs random sampling to obtain the spatial heat distribution simulation results; the heat concentration distribution coordinates are extracted based on the spatial heat distribution simulation results, and the root region of the spade-tooth heat sink fins is located based on the heat concentration distribution coordinates; the local thermal impedance and target heat transfer cross-sectional area are extracted for the root region of the spade-tooth heat sink fins, and the equivalent thermal conductivity is calculated based on the local thermal impedance and the target heat transfer cross-sectional area; if the equivalent thermal conductivity is less than a preset thermal conductivity threshold, a geometric feature vector is extracted for the root region of the spade-tooth heat sink fins, and the geometric feature vector is processed using a pre-established enhancement strategy matching model to output the spade-tooth heat sink fin enhancement processing sequence;The heat dissipation fin reinforcement sequence is processed using a pre-established heat conduction evaluation model to calculate and predict heat conduction efficiency. Based on the predicted heat conduction efficiency, an overall heat dissipation index optimization path is determined. Further, the process of simulating fatigue damage evolution and outputting a complete reinforcement protocol through the overall heat dissipation index optimization path includes: extracting the service load spectrum based on the overall heat dissipation index optimization path; extracting the cyclic stress amplitude based on the service load spectrum; calculating the cumulative plastic strain based on the cyclic stress amplitude and a preset material constitutive relation to obtain the fatigue damage evolution trajectory; extracting the damage accumulation rate based on the fatigue damage evolution trajectory; fitting the life decay curve using the least squares method based on the damage accumulation rate to obtain a durability quantification index; determining if the durability quantification index is greater than a preset damage tolerance threshold; matching the corresponding reinforcement parameter combination in a preset reinforcement database based on the durability quantification index; constructing a protocol mapping matrix based on the reinforcement parameter combination; extracting observation feature vectors from the protocol mapping matrix; decoding the observation feature vectors using a pre-established hidden Markov model to obtain a state transition sequence; and outputting a complete reinforcement protocol based on the state transition sequence.

[0007] The spade-tooth heat sink is manufactured using a method that optimizes the overall performance of spade-tooth heat sinks.

[0008] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects:

[0009] This invention aims to address the challenges of weakened material strength, high crack initiation probability, and insufficient long-term service durability at the fin root under vibration and thermal stress. The invention simulates stress distribution using finite element analysis, determines a uniform pressure application path using particle swarm optimization, and predicts crack probability using Monte Carlo simulation to accurately identify high-risk areas. Subsequently, a reinforcement model is constructed using response surface methodology, parameters are iteratively refined, and the durability improvement effect is simulated. If the effect is insufficient, the optimization algorithm is backtracked and adjusted, ultimately generating a reinforcement processing sequence and an overall optimized heat dissipation performance path. Through multi-algorithm collaboration and iterative optimization, this invention ensures significant improvement in stress relief and fatigue durability at the fin root, achieving stable performance and extended lifespan for heat sinks in long-term service scenarios. Attached Figure Description

[0010] Figure 1 This is a flowchart of the method for optimizing the overall performance of the spade-shaped heat sink of the present invention;

[0011] Figure 2 This is a schematic diagram of the module framework for step S104 in the method for optimizing the overall performance of the spade-shaped heat sink of the present invention.

[0012] Figure 3This is a schematic diagram of the module framework for step S107 in the method for optimizing the overall performance of the spade-shaped heat sink of the present invention. Detailed Implementation

[0013] The technical solutions of the embodiments of the present invention will be clearly and thoroughly described below with reference to the accompanying drawings. The described embodiments are merely some embodiments of the present invention.

[0014] like Figures 1-3 The method for optimizing the overall performance of the shovel-tooth heat sink in this embodiment may specifically include:

[0015] Step S101: By obtaining the geometric model data and historical stress distribution records of the root of the shovel-tooth heat sink fin from the preset database, the finite element analysis method is used to simulate the weakening of material strength under vibration challenges and thermal stress, and a preliminary stress magnification distribution map is obtained.

[0016] A geometric model and historical stress distribution records of the root of the spade-shaped heat sink fins are obtained. The root contour is extracted from the geometric model, and its geometric parameters and material properties are mapped to the historical stress distribution records to obtain an initial state matrix. The mapping process uses nearest-neighbor interpolation to interpolate the stress values ​​at measurement points in the historical stress distribution records to the node positions of the root contour. Each node corresponds to a set of parameters, including node coordinates, local thickness, elastic modulus, Poisson's ratio, and historical peak stress, forming an n-row, 7-column initial state matrix, where n is the total number of discrete nodes in the root contour. A finite element mesh model is constructed based on the initial state matrix. The mesh topology is generated based on the node coordinates in the matrix, and the material properties of each node are assigned to the corresponding mesh elements. The finite element mesh elements contained in the finite element mesh model are obtained. Vibration loads and thermal stress fields are applied to the finite element mesh model, and the nodal displacements of the finite element mesh elements under the combined action of the vibration loads and thermal stress fields are determined using the finite element method. The strain gradient of the finite element mesh element is calculated using the nodal displacements. The finite difference method is employed, calculating the difference in displacement between adjacent nodes divided by the distance between nodes to obtain the strain components εx, εy, and εz. The strain components are then differentiated along the spatial direction to obtain the components of the strain gradient tensor. The Frobenius norm of the tensor is taken as the scalar value of the strain gradient. If the strain gradient is greater than a preset strain gradient threshold of 0.003, fatigue damage is determined to have occurred in the finite element mesh element. For the finite element mesh element that has experienced fatigue damage, the corresponding stiffness degradation coefficient is obtained. This coefficient is defined as D = 1 - e^(-0.5N), where N is the number of accumulated load cycles for the element. The strength reduction ΔE of the finite element mesh element is calculated based on the stiffness degradation coefficient. The strength reduction ΔE is equal to the initial elastic modulus E0 multiplied by the stiffness degradation coefficient D. The material properties in the initial state matrix are updated based on the strength reduction, and the elastic modulus of the corresponding element is updated to E0 minus ΔE to obtain the degradation state matrix. The finite element method is re-solved using the degraded state matrix. The degraded material properties are substituted into the finite element equations, and the updated stress distribution of the root profile is recalculated under the same load boundary conditions. The updated stress distribution is compared with the historical stress distribution record. First, the historical stress distribution record is expanded to the same spatial resolution as the updated stress distribution through cubic spline interpolation. Then, the stress values ​​at the same spatial location are divided point by point to obtain a preliminary stress magnification distribution map. This distribution map is used for subsequent crack propagation path prediction and life assessment.

[0017] For example, this embodiment provides a stress analysis method for the root of the fins of a spade-shaped heat sink. This method is mainly applied to the reliability assessment of high-power server heat dissipation modules. Through multi-physics coupling analysis, it accurately locates the weak points of the heat sink under complex operating conditions.

[0018] First, geometric model data and historical stress distribution records of the fin roots of the spade-shaped heat sink are obtained. The geometric model data includes three-dimensional spatial topological information such as the heat sink substrate thickness, fin spacing, fin height, and fin root radius. This data is exported to STEP or IGES format using computer-aided design software. STL format can be used when rapid modeling is required. STEP format retains complete surface information and is suitable for high-precision analysis, IGES format has better compatibility, and STL format has a smaller file size but lower accuracy. Historical stress distribution records refer to the coordinates of stress concentration areas and corresponding peak stress data accumulated during past actual operation tests of this heat sink model or similar models with fin height differences of less than 10% and substrate thickness differences of less than 15%, obtained through strain gauge testing or infrared thermal imaging. This historical data is used to set initial load boundary conditions when building the finite element model. After the simulation calculation is completed, the calculation results are compared with the historical peak stresses to verify the model's accuracy. Based on the obtained data, a discretized mathematical model for finite element analysis is further constructed.

[0019] Specifically, the finite element method (FEM) is a numerical computation technique that divides a continuous geometric entity into a finite number of interconnected micro-elements and approximates the physical response of the overall structure by solving the mechanical equations of each element. In this embodiment, a static pre-analysis after initial mesh generation is performed to identify the stress gradient in the fin root region. When the stress difference between adjacent elements exceeds 15% of the average stress value, the region is identified as a region with a severe stress gradient. For such regions, local mesh refinement is performed within a range of 2 to 3 millimeters extending outward from the fin root, controlling the element size to 0.3 to 0.5 millimeters, which is 3 to 5 times denser than the 1 to 2 millimeter element size in conventional regions. For regular geometric regions, hexahedral elements are preferred to obtain higher computational accuracy, while tetrahedral elements are used for complex curved surfaces such as root fillets to adapt to the geometry. The material property parameters are assigned in accordance with the principle of consistency with the heat sink substrate material. When the substrate material is pure copper, the elastic modulus E is assigned to be 110 to 130 GPa, the Poisson's ratio μ to be 0.33 to 0.35, and the coefficient of thermal expansion α to be 16 to 18 × 10⁻⁶ per degree Celsius. When the substrate material is aluminum alloy, the elastic modulus E is assigned to be 68 to 72 GPa, the Poisson's ratio μ to be 0.32 to 0.34, and the coefficient of thermal expansion α to be 23 to 25 × 10⁻⁶ per degree Celsius. Yield strength curve data as a function of temperature from 20 degrees Celsius to 150 degrees Celsius are also input to ensure that the discretized model can truly reflect the mechanical characteristics of the physical entity.

[0020] The finite element model constructed above was used to simulate vibration challenges and the effects of thermal stress. Vibration challenges mainly originated from the periodic excitation force generated by the high-speed cooling fan inside the server chassis, as well as random vibrations encountered during equipment transportation. Vibration signals were measured at the heatsink mounting location using a triaxial accelerometer. The frequency domain response was obtained through Fast Fourier Transform, and the acceleration power spectral density curve within the range of 10Hz to 2000Hz was extracted. Periodic excitation force corresponds to discrete peak frequencies, while random vibration corresponds to a continuous spectrum distribution. This curve was used as the excitation load applied to the mounting holes of the heatsink substrate. The dynamic displacement and alternating stress response of the fins at different resonant frequencies were calculated using the modal superposition method. The simulation of thermal stress involved dividing the power consumption of the high-power heating chip by the chip contact area to obtain the heat flux density, which was applied as a heat source boundary condition to the bottom of the heatsink. The convective heat transfer coefficient h was calculated based on the empirical correlation h = C × Re^m × Pr^n, where Re is the Reynolds number, Pr is the Prandtl number, and C, m, and n are coefficients related to the fin geometry. For a typical aluminum finned heat sink, the convective heat transfer coefficient ranges from 50 W / m²·K to 120 W / m²·K when the wind speed is between 2 m / s and 5 m / s. The temperature gradient distribution inside the heat sink is obtained by solving the steady-state heat conduction equation. Since the fin root connects the high-temperature substrate to the relatively cooler distal end of the fin, this temperature gradient causes uneven thermal expansion deformation in different parts of the material, inducing thermal stress at the fin root. When the stress-strain relationship of the material remains linear and the deformation is less than 5% of the fin thickness, a linear superposition method is used to synthesize the vibrational alternating stress and thermal stress. When the temperature gradient exceeds 80℃ / mm or the dynamic displacement amplitude exceeds 5% of the fin thickness, the material enters the plastic deformation stage. A nonlinear superposition method is used, considering the coupling effect of geometric nonlinearity and material nonlinearity, to realistically reproduce the complex stress state of the heat sink in actual service environments.

[0021] Understandably, under the long-term action of the aforementioned combined load, the material will experience a decrease in strength. This decrease in strength refers to the physical process by which a material's ability to resist deformation and fracture gradually declines over time under the combined effects of alternating stress and high temperature.

[0022] Specifically, the simulation process for material strength reduction includes fatigue damage accumulation analysis and high-temperature softening assessment. First, the stress time history data at the fin root obtained from multiphysics coupling calculations is extracted as input. The complex random stress waveform is decomposed into a series of stress cycles using the rainflow counting method, each cycle having a defined stress amplitude and average stress value. Then, combining the material's SN curve (i.e., the correspondence between stress amplitude and cycle number), the damage component Di = ni / Ni caused by each stress cycle is calculated using Miner's linear fatigue damage accumulation rule, where ni is the actual number of cycles at that stress amplitude, and Ni is the fatigue life cycle number of the material at that stress amplitude. The damage components of all stress cycles are summed to obtain the total damage degree D = ΣDi. When the total damage degree reaches 1.0, local microcracks are considered to have initiated in the material, macroscopically manifested as a decrease in load-bearing capacity. On the other hand, the high-temperature softening assessment, based on the local temperature field obtained from thermal analysis, retrieves the yield strength σs(T) at the corresponding temperature from a pre-established material property database. At high temperatures, the lattice slip resistance of metallic materials decreases, and their yield strength decreases significantly. The yield strength σs(T) retrieved in real time is compared with the current composite stress σ to calculate the remaining strength margin M = σs(T) - σ / σs(T). By combining the assessment of fatigue microcrack initiation and high-temperature yield strength reduction, the degree of strength degradation of the fin root material under harsh working conditions is quantified.

[0023] Based on the simulation results of the aforementioned material strength reduction, the system extracts the equivalent stress value σ_eq of all finite element nodes at the fin root, and simultaneously calculates the average stress σ_nom of the heat sink substrate as the nominal stress reference. The stress amplification factor is defined as K_t=σ_eq / σ_nom. For example, if the equivalent stress of a node is 180 MPa and the average stress of the substrate is 60 MPa, then the stress amplification factor for that node is 3.0. After the system completes the calculation by traversing all nodes, it maps the stress amplification factor to a color gradient according to its value. Amplification factors of 1.0 to 2.0 are displayed as blue to green, 2.0 to 3.5 as yellow to orange, and above 3.5 as red. This is then rendered on the 3D geometric model to form a stress amplification distribution map. This distribution map identifies areas of high stress concentration, such as the transition zone at the fin root. The red areas with amplification factors exceeding 3.5 are the key target areas for subsequent structural optimization. By adjusting the radius of the transition fillet at the fin root or optimizing the tooth profile geometry parameters, the stress amplification factor in this area can be reduced to below 2.5, thereby improving the overall reliability of the heat sink.

[0024] Step S102: Based on the preliminary stress magnification distribution diagram, and considering the tensile force formation characteristics in the narrow space at the fin root, a particle swarm optimization algorithm is used to determine the uniform pressure application path. If the path coverage is lower than a preset threshold, the parameters are iteratively adjusted to obtain the optimized pressure distribution scheme.

[0025] A preliminary amplified stress distribution map is obtained and input into the finite element analysis module for subsequent stress concentration identification. Regions with stress values ​​exceeding 0.6 times the material's yield strength are identified from the preliminary amplified stress distribution map as stress concentration areas, located within the narrow space at the fin root. The principal stress directions and values ​​are obtained from the stress tensor decomposition of these stress concentration areas. The tensile stress components and their directions constitute the tensile force characteristics, including the magnitude of the tensile stress, the tensile stress direction angle, and the coordinates of the point of application. The nodal displacement data of the stress concentration area under tensile force is obtained through finite element analysis as deformation data. A surface fitting method is used to perform cubic spline interpolation on the deformed nodal coordinates. The Gaussian curvature K of the fitted surface is calculated to be equal to k1 × k2 as the spatial curvature of the narrow space, where k1 and k2 are the two principal curvatures at that point. The pressure-applying component is predefined as a pressure block with an arc-shaped contact surface, and its surface curvature is measured using a coordinate measuring machine. The spatial curvature K is numerically compared with the surface curvature Kp of the pressure-applying component. When the absolute value of the curvature difference |K - Kp| is less than 0.05 per square millimeter, curvature matching is determined. At this time, the arc surface of the pressure-applying component is in contact with the deformed surface at the root of the fin. The elastic contact area A is calculated using Hertzian contact theory, which is equal to π multiplied by a multiplied by b, where a and b are the major and minor semi-axes of the contact ellipse, respectively. This is obtained through iterative solutions using the material's elastic modulus and contact load. The initial pressure angle θ0 is determined based on the angle between the geometric center of the contact area and the tensile stress direction angle in the tensile characteristic, ensuring that the deviation angle between the pressure direction and the principal tensile stress direction is less than 15 degrees. The load gradient G is calculated using the initial pressure angle θ0 and a preset initial pressure value P0 equal to 50 MPa. The load gradient G is equal to the square root of P0 divided by A, in MPa per millimeter. Based on the load gradient G, multiple discrete pressure nodes are generated within the contact area A in a hexagonal close-packed arrangement. The node spacing d is equal to 2 mm, and the number of nodes N is equal to A divided by the square of d and rounded up. The pressure nodes are connected according to the nearest neighbor principle to construct a pressure path, forming a spiral path progressing from the edge of the contact area towards the center. The node coordinate sequence of the pressure path is used as the particle position vector. The path coverage C is defined as the area covered by the path divided by the total contact area A, and C is used as the fitness function. A particle swarm optimization algorithm is used to optimize the pressure path. Particle velocity updates use the standard particle swarm optimization velocity formula, and position updates are achieved by superimposing the velocity vector on the current node coordinates. The inertia weight is set to 0.7, and the learning factors c1 and c2 are both set to 1.5. After 100 iterations, the path coverage set is obtained.Calculate the path coverage rate C for each path in the path coverage set. If the path coverage rate C is lower than a preset coverage threshold of 0.85, iteratively adjust the initial pressure angle θ0, with each adjustment step being 5 degrees. The adjustment direction is to make the pressure path closer to the tensile stress direction. Based on the iteratively adjusted pressure angle, recalculate the load gradient and regenerate the pressure nodes to perform particle optimization. During the optimization process, simultaneously calculate the pressure distribution uniformity index U, where U equals 1 minus the pressure standard deviation divided by the average pressure value. When U is greater than 0.9 and C is greater than 0.85, the optimized pressure distribution scheme is obtained.

[0026] After obtaining a preliminary stress magnification distribution map, the system identifies stress concentration areas within the narrow space at the root of the fins of the spade-shaped heat sink based on this map. Specifically, areas with a stress magnification factor K greater than or equal to 2.5 are extracted from the distribution map as high-risk stress concentration areas, corresponding to the peak tensile stress locations at the junction of the fin root and the substrate. The system transmits the coordinate information and stress values ​​of these high-risk areas to the subsequent structural optimization module to determine the specific locations where fillet transitions or thickness adjustments are needed, thereby reducing the tensile stress level at the fin root.

[0027] The narrow space at the fin root refers to the tiny gap region close to the substrate between two adjacent spade teeth. The spade teeth are arrayed along the substrate surface, with a spacing of 0.8 to 1.2 mm between adjacent teeth. The substrate is a metal base plate supporting the fins, with a thickness of 2 to 3 mm. Tensile stress formation characteristics within this region are determined by the following indicators: a local tensile stress exceeding 60 MPa and a strain exceeding 0.3% in the normal direction indicate a significant tensile tendency. When the tensile stress is below 40 MPa or the strain is less than 0.15%, no significant tensile stress formation characteristics are considered.

[0028] Specifically, the system traverses high-stress nodes from the preliminary amplified stress distribution map obtained in step S101. This distribution map records the principal tensile stress value, the three-dimensional vector coordinates of the stress direction, and the stress gradient value along the fin height direction for each grid node at the fin root. The system constructs a five-dimensional feature vector for each high-stress node, including the node's spatial coordinates x, y, and z, the principal tensile stress value σ, and the stress gradient change rate dσ / dh. The K-means clustering algorithm is used to group these feature vectors. The number of clusters K is determined by the silhouette coefficient method. The silhouette coefficient S under different K values ​​is calculated, and the K value corresponding to the maximum value of S is selected as the optimal number of clusters. The system marks the centroid position of each group of nodes in the clustering results as the tensile core region, and defines the envelope formed by connecting the boundary nodes of each group as the diffusion boundary, thereby reconstructing the spatial topology of tensile concentration within the narrow space. Based on the extracted tensile formation features, the system further uses a particle swarm optimization algorithm to determine the uniform pressure application path. In real-world business scenarios, external pressure needs to be introduced to perform surface strengthening treatment on the fin roots or to simulate reverse load tests. Planning the movement trajectory for applying pressure is a crucial step.

[0029] The system employs a particle swarm optimization algorithm to plan the pressure application path. The pressure application path is discretized into N consecutive control nodes, with each candidate path represented as a particle. The particle's position vector is X = (x1, y1, z1, x2, y2, z2, ..., x...). n y n , z n This contains the three-dimensional coordinates of N control nodes, where (xᵢ, yᵢ, zᵢ) are the coordinates of the i-th node. x1 ,vᵧ1,vᵨ1,...。 ,v xn ,vᵧ n ,vᵨ n The value represents the iterative adjustment amount of the coordinates of each node. The system initializes 50 particles within the boundary of the narrow space at the root of the fin. Each particle's N control nodes are assigned random initial coordinates and velocity components. The node coordinates are restricted within the boundary of the workspace, and the initial range of the velocity components is set to ±2 mm per iteration.

[0030] To guide the particle swarm in finding the optimal combination of application positions, the system constructs a fitness function to evaluate the uniformity of pressure distribution. During calculation, each particle position is simulated as a point of application of 50 Newtons of pressure. The principal stress σ1 values ​​of the stress field generated at this point are extracted for each grid node within the narrow space. These are algebraically added to the tensile stress values ​​of the corresponding nodes in the previously extracted tensile features to obtain the superimposed stress value σr for each node. The variance Var of σr for all grid nodes is calculated as an indicator of the non-uniformity of the residual stress distribution. A Var less than 2 MPa² is considered to indicate uniform pressure distribution, and the fitness value is defined as f = 1 divided by Var plus 0.01. In each iteration, the system records the global optimal position set of the particle swarm, which contains 8 to 12 discrete application points. After the algorithm converges, these points are connected according to their angular order in the circumferential direction within the narrow space to form a closed or semi-closed pressure application path profile.

[0031] Understandably, obtaining only the initial path is insufficient for direct application in practice; its spatial coverage capability also needs rigorous evaluation. Therefore, the system determines whether the coverage of the initial path is below a preset threshold.

[0032] In one embodiment, the path coverage calculation process is as follows. First, the pressure application radius is determined. This radius is set according to the physical size of the grinding head, typically 1.2 to 1.5 times the grinding head radius; in this embodiment, it is 12 mm. The system extends the initial path in the normal direction along this radius. Specifically, a sampling point is taken every 5 mm along the path, and a sphere is constructed with this point as the center and a radius of 12 mm. The union of all spheres forms a continuous three-dimensional application surface. Then, this surface is orthogonally projected onto the physical surface of the narrow space at the fin root. The number of triangular mesh patches N_c in the projected area is counted, and the coverage C = N_c / N_t is calculated, where N_t is the total number of mesh patches in the narrow space. If C is less than 0.9, it indicates that there are obvious omissions in the currently planned path, and some tension concentration points are not effectively covered. In response to this situation, the system automatically triggers a parameter iteration adjustment mechanism.

[0033] For example, during the iterative parameter adjustment process, the system focuses on modifying the control variables of the particle swarm optimization algorithm. When the improvement in the global optimal fitness value is less than 0.01 after five consecutive iterations, the system determines that the algorithm has fallen into a local optimum. At this point, the system linearly increases the inertia weight parameter w from its current value to 1.2 times its initial value. The calculation formula is w = w_current + 0.04 multiplied by the current iteration number, where w_current is the current inertia weight value. This allows particles to maintain a higher flight speed in subsequent iterations, thereby enhancing the algorithm's exploration ability in the global space and encouraging particles to fly out of the currently dense local area to explore uncovered blank areas. Simultaneously, the system adjusts the self-learning factor c1 from its initial value of 2.0 to 2.5 and the social learning factor c2 from its initial value of 2.0 to 1.5. By increasing the ratio of c1 to c2, the system reduces the tendency of particles to converge towards the swarm's optimal position and encourages particles to search independently. Using these adjusted parameters, the system restarts the particle swarm optimization iteration process, generates new candidate paths, and performs coverage calculation and judgment again.

[0034] The parameter adjustments and path generation employ an iterative optimization approach, calculating the coverage rate of candidate paths after each iteration. Iteration stops when the coverage rate reaches 95%, and the resulting path scheme is output as the final result. If the coverage rate remains below 90% after 50 iterations, the process terminates and is marked as convergence failure, requiring adjustment of the initial parameters and recalculation. The output path scheme includes a three-dimensional spatial coordinate sequence. The velocity control strategy is determined based on the ratio of the distance between adjacent nodes to the set total duration, and the pressure application vector is calculated from the tissue deformation threshold corresponding to the node location.

[0035] Step S103: Extract the coordinates of key nodes from the optimized pressure distribution scheme, obtain the corresponding fatigue damage risk index, predict the crack initiation probability using the Monte Carlo simulation method, and determine the list of high-risk areas.

[0036] Locations with stress concentration factors greater than 2.5 are identified as key nodes from the optimized pressure distribution scheme, and their coordinate data is extracted. Fatigue damage indices corresponding to these key node coordinates are obtained through finite element analysis, including equivalent stress amplitude, average stress, and cumulative damage degree. Rainflow counting is used to statistically analyze these fatigue damage indices, extracting the stress alternation amplitude Δσ and the corresponding load cycle number n. A stress probability distribution function is constructed using the statistical frequency of the stress alternation amplitude and the load cycle number. Randomly sampled stress samples are generated using this stress probability distribution function, and a preset material fracture toughness KIC is obtained. Based on fracture mechanics theory, with the stress intensity factor K equal to the material fracture toughness KIC as a critical condition, the initial crack length a0 is inferred from the randomly sampled stress samples. The crack propagation rate is calculated using the Paris formula da / dN = C × ΔK^m, where da / dN is the crack propagation rate, C and m are material constants, and ΔK is the stress intensity factor amplitude, calculated based on the initial crack length and the randomly sampled stress samples. The preset number of simulation iterations is set to 10,000. Using the Monte Carlo simulation method, in each iteration, the crack length increment Δa is calculated based on the crack propagation rate and accumulated to the current crack length. When the accumulated crack length reaches the critical crack length ac, it is recorded as one failure event. The ratio of the total number of failure events within the preset number of iterations to the number of iterations is calculated to obtain the crack initiation probability corresponding to the key node coordinate data. If the crack initiation probability is greater than a preset probability threshold of 0.05, the corresponding key node coordinate data is extracted, and a high-risk zone list is generated based on the key node coordinate data.

[0037] Based on the optimized pressure distribution scheme obtained from the aforementioned iterative optimization process, the system further conducts an in-depth assessment of the long-term service reliability of the fin root.

[0038] Specifically, the system extracts the coordinates of key nodes from the optimized pressure distribution scheme. These key nodes typically correspond to stagnation points, points of abrupt changes in direction, and extreme points in the initial stress field along the pressure application path. Their geometric positions in three-dimensional space form the basic data framework for subsequent fatigue analysis.

[0039] The system extracts the coordinates of key nodes and obtains the corresponding fatigue damage risk index for heat sinks with fins made of 6061 aluminum alloy or T2 copper. The fatigue damage risk index is quantified through cumulative damage measurement, and the specific calculation process is as follows: The system first extracts the stress amplitude S_a and mean stress S_m at the key node. The stress amplitude is half the difference between the maximum and minimum stresses in the thermo-coupling cycle, and the mean stress is the arithmetic mean of the two. The system calls the SN curve of the corresponding material in the material database, which records the number of failure cycles N_f under different stress amplitudes. For cases where the mean stress is non-zero, the system uses the Goodman correction formula to equate the actual stress state to the equivalent stress amplitude S_eq under zero mean stress. The calculation formula is S_eq equal to S_a divided by 1 minus the ratio of S_m to the tensile strength σ_b. Substituting S_eq into the SN curve yields the number of failure cycles N_f corresponding to that stress level. The system counts the actual number of cycles n experienced by the key node within the expected service life and calculates the cumulative damage degree D according to Miner's linear cumulative damage theory. D equals n divided by N_f. When D is greater than or equal to 1, it indicates that the node has reached the critical failure state. Considering the inhomogeneity of microstructure during manufacturing and random disturbances in the service environment, the system introduces Monte Carlo simulation to predict the crack initiation probability.

[0040] Specifically, in performing Monte Carlo simulations to predict crack initiation probability, the system first constructs a probability distribution model containing multiple random variables. At the microscale of the narrow space at the fin root, the initial fracture toughness of the material follows a Weibull distribution with a shape parameter of 2.5 and a scale parameter of 85 MPa² / m². The equivalent diameter of the micropores follows a log-normal distribution with a mean of 3.2 μm and a standard deviation of 0.8 μm. The peak fluctuation of the external thermal stress cycle follows a normal distribution with a mean of 450 MPa and a standard deviation of 35 MPa. The microscopic fluctuation of the material's yield strength follows a normal distribution with a mean of 520 MPa and a standard deviation of 28 MPa. The system uses a Mason rotation pseudo-random number generation algorithm to perform 10,000 independent samplings within the above probability distribution interval, generating 10,000 sets of virtual samples representing different microstates and load conditions. For each generated virtual sample, the system substitutes it into the Griffith energy balance criterion for verification. The criterion first extracts the fracture toughness K, micropore diameter d, and peak thermal stress σ of the current sample as input parameters. Then, it calculates the strain energy density release rate at the crack tip, G = (σ² × π × d) divided by Young's modulus E, where E is taken as 180 GPa. G is then compared with the critical strain energy density release rate Gc corresponding to the sample. Gc is calculated based on the fracture toughness K of the current sample using the formula Gc = (K²) divided by E. If G is greater than or equal to Gc, the system determines that the sample meets the thermodynamic conditions for crack propagation and marks this sampling as a crack initiation event in the system record. After completing all 10,000 simulation calculations, the system summarizes the total number of crack initiation events.

[0041] Understandably, after completing simulations and verifications on a massive sample size, the system divides the total number of crack initiation events recorded at each key node by the total number of samplings in the Monte Carlo simulation, thereby accurately calculating the crack initiation probability corresponding to each key node. Based on these probability values, the system sorts all key nodes in descending order. The system uses a tiered threshold method to set safety warning standards, specifically: first, based on industry failure criteria, the crack initiation probability is divided into three levels: a probability greater than 0.15 is considered high-risk, between 0.08 and 0.15 is considered medium-risk, and below 0.08 is considered acceptable. Then, combining the fatigue limit data of the fin material and historical failure case statistics, the safety warning threshold is set to 0.15. The system automatically classifies nodes with a crack initiation probability higher than 0.15 and their adjacent small areas as high-risk, ultimately generating a detailed list of high-risk areas. This list clearly indicates the spatial location at the fin root where structural failure is most likely to occur under long-term complex stress, providing precise targeting for subsequent structural reinforcement or process improvement.

[0042] Step S104: For the list of high-risk areas, a precise reinforcement model is constructed using the response surface methodology. If the model fitting error exceeds a preset threshold, sample points are added from historical data to obtain a refined reinforcement parameter set.

[0043] Obtain the region nodes corresponding to the high-risk region list, extract geometric features and load boundaries based on the region nodes, obtain preset initial reinforcement parameters, and construct an initial response surface using the response surface methodology based on the geometric features, load boundaries, and initial reinforcement parameters. Extract 20 independent samples from the design space using the Latin hypercube sampling method as a validation sample set. Calculate the predicted response values ​​corresponding to the validation sample set using the initial response surface. Calculate the fitting error based on the predicted response values ​​and the true response values ​​of the validation sample set. If the fitting error is greater than 0.05, obtain historical reinforcement samples, which are sample points with errors exceeding 0.08 in the previous iteration and their corresponding true response values. Correct the initial response surface by adding the historical reinforcement samples to the training set and refitting the response surface to obtain the target response surface. The partial derivatives of the response function R in the target response surface with respect to each enhancement parameter p are calculated to obtain the parameter sensitivity matrix S = ∂R / ∂p. The response gradient g = ∇R is calculated, and the parameter compensation vector Δp ​​= -α·g is calculated using the gradient descent method, where α is the learning rate. The iteration step size λ = min(0.1 / ||g||, 0.5) is adaptively determined based on the current gradient magnitude. The initial enhancement parameters are updated according to the iteration step size and the parameter compensation vector, with the update formula p_new = p_old + λ·Δp, to obtain refined enhancement parameters. The iteration is terminated when the parameter change is less than 0.01 for three consecutive iterations or when the number of iterations reaches 50. The target enhancement scheme is generated based on the refined enhancement parameters.

[0044] In one possible implementation, the system obtains a list of high-risk areas generated in the aforementioned steps and initiates an optimization process for structural reinforcement parameters for the spatial locations in the list.

[0045] Specifically, the system uses the response surface methodology to construct a precise reinforcement model in order to establish a mathematical mapping relationship between process adjustment parameters and the degree of local stress relief.

[0046] The precision reinforcement model reduces the risk of fatigue damage in high-risk areas by fine-tuning the local geometry or surface treatment of the heat sink. The root transition region at the junction of the heat sink fins and the substrate experiences periodic thermal stress concentration. The local substrate refers to the base plate area directly transmitting heat flow below the fins. The model's input variables include the radius of curvature of the fin root transition fillet, the thickness of the local substrate, and the compressive stress depth of the surface hardening treatment. The system uses these physical parameters as independent variables Ri, where R1 is the radius of curvature of the fillet, ranging from 0.5 to 3.0 mm; R2 is the thickness of the substrate, ranging from 0.2 to 1.5 mm; and R3 is the compressive stress depth, ranging from 0.1 to 0.8 mm. The target response value is the percentage decrease in stress amplitude η in the high-risk area, calculated as η = (σ0 - σ1) / σ0, where σ0 is the stress amplitude before reinforcement and σ1 is the stress amplitude after reinforcement. This reduction ratio directly reflects the degree of stress relief. When the radius of the fin root fillet increases from 0.5 mm to 2.0 mm, the stress amplitude can decrease by 25% to 40%.

[0047] For example, in constructing a precise reinforcement model using the response surface methodology, the system first selects representative initial sample points within the value space of the reinforcement parameters using experimental design theory. The system typically employs a central composite design strategy, extracting several sets of parameter combinations within the upper and lower limits of each reinforcement parameter according to specific geometric distribution rules. The number of sample points is usually 15 to 30 sets. For each extracted parameter combination, the system calls the underlying finite element simulation module for virtual calculation. This module converts the reinforcement parameters into boundary condition inputs, obtains the local stress distribution in the high-risk area after applying the set of reinforcement parameters by solving the stress balance equation, and calculates the corresponding stress relief index, defined as the percentage decrease in stress amplitude after reinforcement relative to stress amplitude before reinforcement. Furthermore, based on the initial sample points and their corresponding simulation results, the system uses a second-order polynomial regression algorithm to fit and generate a response surface function. This function includes linear, quadratic, and cross terms, and is in the form that η equals a0 plus the sum of the products of each strengthening parameter pi and its coefficient ai, plus the sum of the squares of each strengthening parameter square term pi and its coefficient bii, plus the sum of the products of each strengthening parameter cross term pi multiplied by pj and its coefficient bij, where η represents the stress relief index, pi and pj represent different strengthening parameters, and a0, ai, bii, and bij are regression coefficients. This function, in the form of a continuous mathematical surface, approximately expresses the complex nonlinear relationship between the set of strengthening parameters and the stress relief index. By constructing this continuous surface, the system can quickly predict the structural strengthening performance under arbitrary parameter combinations with lower computational cost compared to direct simulation.

[0048] To ensure the reliability of the model's predictions, the system requires rigorous accuracy verification of the fitted response surface function. The system uses leave-one-out cross-validation to calculate the root mean square error (RMSE) as an evaluation metric for model fitting error, where RMSE = √(Σ(yi-ŷi)² / n), yi is the actual value, ŷi is the predicted value, and n is the number of samples. The system compares the calculated RMSE value with a preset threshold. When the RMSE is less than 0.05, the model is considered to have sufficient generalization ability; when the RMSE is greater than or equal to 0.05, additional sample points are needed to reconstruct the response surface model. This threshold is set according to engineering accuracy requirements; for high-precision prediction scenarios, it can be set to 0.02 to 0.03, while for general-precision scenarios, it can be relaxed to 0.08 to 0.1.

[0049] In one embodiment, if the model fitting error exceeds a preset threshold, it indicates that the current number of initial sample points is insufficient or their distribution is unreasonable, making it impossible to accurately capture local nonlinear abrupt changes in the parameter space. At this time, the system triggers a sample supplementation mechanism to supplement sample points from historical data.

[0050] The system accesses the enterprise's internal heatsink process simulation database, which contains historical analysis records. Each record stores the heatsink's geometric parameters, reinforcement parameter combinations, and corresponding temperature field test data. The geometric parameters include four dimensions: heatsink height, thickness, spacing, and substrate size. During retrieval, the system calculates the Euclidean distance between the current high-risk area's geometric parameter vector and the historical geometric parameter vectors in the database. The distance D is calculated as D = (h1-h2)² + (t1-t2)² + (s1-s2)² + (b1-b2)², where h is the height, t is the thickness, s is the spacing, b is the substrate size, subscript 1 indicates the current parameter, and subscript 2 indicates the historical parameter. The system extracts the reinforcement parameter combinations and their corresponding actual temperature field test data from these historical records. The temperature field data is stored in the form of a heatsink surface temperature distribution matrix. The system identifies areas in the current response surface where the relative error between the predicted and simulated temperatures exceeds 8% as high-error areas. It calculates the Mahalanobis distance between the parameter vector corresponding to the center point of this area and the historical parameter vectors, filtering historical parameter combinations with a Mahalanobis distance less than 2.5. The system uses these selected historical parameter combinations as supplementary sample points and injects them directly into the current training set.

[0051] Preferably, after supplementing historical sample points, the system re-executes multivariate polynomial regression fitting using the expanded overall sample set to update the coefficients of the response surface function. The system calculates the root mean square error (RMSE) of the updated model, which is the square root of the average of the squares of the differences between the predicted and measured values, and determines whether this error is lower than a preset threshold of 0.05. If it is not, the iterative process of error judgment and sample supplementation is repeated until the model fitting error is strictly lower than the threshold. This dynamic adaptive modeling mechanism effectively avoids prediction distortion caused by insufficient samples. Based on the above iteratively converged accurate reinforcement model, the system uses an adaptive particle swarm optimization algorithm to perform global optimization in the continuous parameter space. By initializing the particle swarm and setting the position and velocity of each particle, the system iteratively updates the position of each particle to track the individual optimal and global optimal solutions, eventually converging to the parameter coordinates that maximize the stress relief index. The system maps this optimal coordinate to the corresponding physical parameter values, including process parameters such as reinforcement pressure, holding time, mold temperature, and feed rate, thereby obtaining a refined reinforcement parameter set. This parameter set provides accurate data support for subsequent mold trimming and processing program development.

[0052] Step S105: Based on the refined set of enhancement parameters, simulate the durability improvement effect, obtain the stress amplification and mitigation index under virtual testing, and determine if the index improvement rate is insufficient to the preset threshold. Then, backtrack and adjust the initial population of the particle swarm optimization algorithm to obtain an updated pressure distribution scheme.

[0053] After obtaining the refining and strengthening parameters, a load history sequence is extracted from the historical operation database. This sequence contains four types of data: timestamp, load amplitude, load frequency, and stress ratio, with a sampling frequency of 100Hz. The residual stress value σr from the refining and strengthening parameters and the stress amplitude σa from the load history sequence are substituted into the Miner linear cumulative damage criterion to construct a fatigue damage assessment model: D = Σni divided by Ni, where ni is the number of stress cycles at stress level i, and Ni is the fatigue life corresponding to that stress level. The cumulative fatigue damage value D is calculated using this model. Regions where D is greater than 0.8 are identified as stress concentration areas. The stress concentration factor Kt0 of the original structure without refining and strengthening is obtained as a benchmark test index. This index is obtained by extracting the ratio of maximum stress to nominal stress through finite element static analysis. The stress concentration factor Kt1 after refining and strengthening is extracted for stress concentration areas as a stress amplification and mitigation index. The improvement ratio R is calculated as: R = Kt0 minus Kt1 divided by Kt0 multiplied by 100%. If the improvement ratio R is less than a preset improvement threshold of 15%, a population backtracking mechanism is triggered. This mechanism extracts particle positions and fitness values ​​from the first five generations of the population from historical optimization records. The top 30% of particles by fitness are selected as elite individuals and retained. A Gaussian perturbation with a standard deviation of 10% of the current search space is applied to the position coordinates of the remaining 70% of particles, generating a new initial population for use by the particle swarm optimization algorithm. For each particle in this initial population, the fitness function is the reciprocal of the cumulative fatigue damage value D. The fitness corresponding to each particle's historical best position is extracted as the individual extreme value pbest, and the fitness corresponding to the best position among all particles is extracted as the global extreme value gbest. The particle update velocity v is calculated according to the standard particle swarm optimization algorithm's velocity update formula, where the inertia weight ω is 0.7 and the learning factors c1 and c2 are both 1.5. The new position x' = x + v is calculated based on the particle update velocity v and the current position x. A ±5% random perturbation is applied to each dimension of the new position x' to complete the position mutation process, resulting in a mutated particle swarm. The position vector of each particle in the mutated particle swarm contains three optimized parameters: shot peening diameter, shot peening velocity, and coverage. These three parameters are mapped to three input parameters for finite element simulation: material reinforcement layer thickness, surface residual stress value, and reinforcement region range. Finite element models are established using the parameter combinations corresponding to each particle in the mutated particle swarm for dynamic contact simulation. Cyclic loads from the load history sequence are applied, and the displacement vector u of each node on the reinforced region surface is extracted within the load cycle. Based on the structural node displacement vector u, the surface contact force F = K multiplied by u is calculated using the principle of virtual work, where K is the structural stiffness matrix. The surface contact force F is normalized according to the area of ​​action to obtain the surface pressure value p = F divided by A, where A is the area of ​​the node influence domain.The pressure values ​​and spatial coordinates of all nodes are aggregated to form pressure distribution field data. Compared with the initial unoptimized pressure distribution field, this data is an updated pressure distribution scheme. This scheme contains four types of information: node number, three-dimensional coordinates, pressure value, and direction of action, which is used to guide subsequent structural optimization design.

[0054] The system obtains the refined set of enhancement parameters output from the aforementioned steps, and initiates a simulation process for improving the durability of high-risk areas of the heat sink based on this set of parameters.

[0055] Specifically, the simulation of durability improvement aims to evaluate the reliability of strengthening parameters under long-term service conditions. The system imports the refined set of strengthening parameters into the underlying finite element analysis module to reconstruct the three-dimensional geometric model and surface physical properties of the heat sink. Subsequently, the system loads a predefined thermal alternating load spectrum, which simulates the frequent start-stop cycles and extreme temperature fluctuations experienced by the heat sink in actual operation. The system uses a transient thermo-mechanical coupling solver to calculate the stress evolution process within the material in high-risk areas under multiple consecutive thermal cycles, thus completing the construction of the virtual test environment. Furthermore, the system extracts the peak stress data and stress distribution gradient of the high-risk areas in each thermal cycle from the virtual test calculation results.

[0056] The system calculates a stress amplification mitigation index based on extracted stress data. This index is defined as the relative decrease rate between the initial stress amplification factor and the current stress amplification factor, quantifying the ability of strengthening parameters to suppress local stress concentration as a percentage. The specific calculation process is as follows: First, without applying strengthening parameters, the system extracts the local maximum stress σmax0 in the high-risk area and the nominal stress σnom in the uniform stress area far from the stress concentration area, calculating the initial stress amplification factor K0 = σmax0 ÷ σnom. Then, after applying strengthening parameters, the system extracts the current local maximum stress σmax1 at the same location, calculating the current stress amplification factor K1 = σmax1 ÷ σnom. The stress amplification mitigation index R = (K0 - K1) ÷ K0 × 100%. For example, if the initial local maximum stress of a welded joint is 320 MPa, the nominal stress is 160 MPa, and the initial stress amplification factor is 2.0, after applying strengthening parameters, the local maximum stress drops to 240 MPa, and the current stress amplification factor is 1.5. Therefore, the stress amplification mitigation index is 25%, directly reflecting the degree of improvement in the local stress concentration phenomenon.

[0057] The system compares the calculated stress amplification mitigation index with a preset threshold to determine whether the current set of strengthening parameters meets stringent durability requirements. The stress amplification mitigation index characterizes the degree of improvement in stress concentration in the weld area. Specifically, it extracts the peak stress σmax and nominal stress σnom at the weld before and after strengthening using finite element simulation. It calculates the stress concentration coefficients before strengthening (Kt1 = σmax1 / σnom1) and after strengthening (Kt2 = σmax2 / σnom2). The stress amplification mitigation index is defined as the relative reduction rate η = (Kt1 - Kt2) / Kt1 × 100%. When η ≥ 15%, the durability requirements are met; when 10% ≤ η < 15%, the requirements are basically met; and when η < 10%, the requirements are not met and the strengthening parameters need to be re-optimized.

[0058] In one embodiment, the system uses the fatigue life improvement rate η as an evaluation index, which is calculated as η = Noptimized - Noriginal / Noriginal × 100%, where Noriginal is the number of fatigue cycles of the original sample, and Noptimized is the number of fatigue cycles after applying the current set of strengthening parameters. The set of strengthening parameters includes a combination of process parameters such as residual stress distribution depth, surface roughness Ra, hardened layer thickness h, and shot peening intensity A. If the fatigue life improvement rate η is determined to be less than 15%, it indicates that the current set of strengthening parameters has failed to achieve the expected fatigue life improvement effect and has fallen into a local optimum. At this time, the system triggers a backtracking mechanism, which expands the search range by reinitializing the positions and velocities of some particles in the particle swarm and intervening in the particle swarm optimization algorithm used for global extremum search to escape the local optimum region.

[0059] For example, the system backtracks and adjusts the initial population of the particle swarm optimization algorithm. In the conventional optimization process, the initial population is usually uniformly and randomly distributed throughout the parameter space. When backtracking is triggered, the system introduces a Logistic chaotic mapping sequence to regenerate the position coordinates of the initial particles. This mapping is generated by the iterative formula xn+1=μ·xn·(1-xn), where x0 is a random number between 0 and 1, and μ is between 3.8 and 4.0. The generated chaotic sequence values ​​are mapped to the value ranges of each parameter dimension through linear transformation to break the original uniform distribution pattern and enhance the traversal diversity of the population in the parameter space. At the same time, the system expands the initialization range of particle velocity. If the velocity range is -10% to +10% of the value range of each parameter dimension during conventional initialization, it is expanded to -25% to +25% after backtracking, giving the particles a stronger ability to escape local extreme regions. In addition, the system will extract particles that ranked in the top 30% of fitness during the historical optimization process but were not retained due to population updates, and inject them as elite individuals into the new initial population. The number of particles injected is 12% of the size of the new population, thereby guiding the population to gather in potentially better areas while ensuring the breadth of the search.

[0060] Preferably, the system first defines a stress amplification mitigation index, which comprehensively evaluates the effect of the effective depth of the residual compressive stress layer and the peak compressive stress amplitude on the suppression of fatigue crack propagation. The calculation formula is η = d multiplied by the square root of σ divided by 100, where d is the compressive stress layer depth in millimeters, and σ is the absolute value of the maximum residual compressive stress in megapascals. A larger η value indicates a better strengthening effect. The system uses the adjusted initial population to restart the iterative optimization process of the particle swarm optimization algorithm, searching again for parameter coordinates that maximize the stress amplification mitigation index within the continuous parameter space of the strengthening parameter set. This strengthening parameter set includes three types of process parameters: shot peening intensity ranging from 0.2 to 0.6 mm Almen intensity, coverage ranging from 100% to 200%, and blasting angle ranging from 60 degrees to 90 degrees. Finally, the system obtains an updated pressure distribution scheme through iterative convergence. This scheme refers to the residual compressive stress sequence at sampling points every 0.05 mm from the surface to a depth of 0.5 mm in high-risk areas of the heat sink, namely the fin roots and weld junctions where the stress concentration factor is greater than 2.5. Compared to the unoptimized uniform distribution assumption in the initial scheme, the updated scheme exhibits a gradient distribution characteristic, gradually transitioning from -600 MPa at the surface to -200 MPa at a depth of 0.3 mm. The system outputs this updated pressure distribution scheme as core control data to guide the dynamic adjustment of impact energy and path planning in subsequent surface shot peening or laser shock peening processes.

[0061] Step S106: Extract the final uniform pressure application path from the updated pressure distribution scheme, combine it with the Monte Carlo simulation results, generate the fin root strengthening treatment sequence, and determine the overall heat dissipation performance optimization path.

[0062] An updated pressure distribution scheme is obtained, which includes the pressure values ​​of each node, the spatial coordinates of the nodes, and the application sequence. This scheme is obtained through iterative optimization of the initial uniform pressure application scheme. Surface pressure nodes are extracted based on the updated pressure distribution scheme. Nodes with pressure values ​​greater than 1.2 times the average pressure are selected as surface pressure nodes. Based on the spatial coordinates of these surface pressure nodes, neighboring nodes within a distance of less than 3 mm are found in the 3D mesh. The pressure values ​​of these neighboring nodes are extracted, and the absolute value of the pressure difference between each surface pressure node and its neighboring nodes is calculated as the local pressure difference. Nodes are sorted from largest to smallest according to their local pressure differences and connected sequentially to form a pressure application path, resulting in the final uniform pressure application path. Contact thermal resistance parameters are extracted for each node along the path of the final uniform pressure application. These parameters include contact surface roughness, contact pressure, and material hardness. Based on the contact thermal resistance parameters, the heat flux density q is calculated using Fourier's law, which is equal to the negative thermal conductivity k multiplied by the temperature gradient gradT. The temperature gradient is equal to the temperature difference divided by the heat transfer distance. The heat flux density and temperature gradient are used as random variables, and a two-dimensional normal distribution is fitted based on historical measurement data as the joint probability distribution function. Monte Carlo simulation is used to perform 10,000 random samplings on the joint probability distribution function. Convergence is determined when the standard deviation of 500 consecutive sampling results is less than 0.01 of the mean, thus obtaining the simulation results of spatial heat distribution. Based on the spatial heat distribution simulation results, the center coordinates of the regions with heat flux densities greater than 1.5 times the average heat flux density are extracted as the heat concentration distribution coordinates. Based on the heat concentration distribution coordinates, the regions with a vertical distance of less than 5 mm from the surface of the fin substrate are located in the fin geometric model and positioned as the fin root region. For the fin root region, the local thermal impedance is calculated by measuring the temperature distribution and heat flux density, which is equal to the temperature difference divided by the heat flux. The cross-sectional area perpendicular to the heat flux direction is measured as the effective heat transfer cross-sectional area. According to the formula, the equivalent thermal conductivity keq is equal to the heat transfer distance L multiplied by the heat flux q divided by the product of the effective heat transfer cross-sectional area A and the temperature difference deltaT. If the equivalent thermal conductivity is less than a preset thermal conductivity threshold of 80 W / m / Kelvin, a geometric feature vector is extracted for the fin root region. This geometric feature vector includes five feature parameters: fin root thickness, fin spacing, fin height, substrate thickness, and radius of curvature. The geometric feature vector is input into a pre-established reinforcement strategy matching model. This reinforcement strategy matching model adopts a deep Q-network structure. The input layer receives the 5-dimensional geometric feature vector, the hidden layer contains three fully connected layers, and the output layer outputs a priority ranking of five processing actions, including increasing fin root thickness, decreasing fin spacing, adding heat-conducting pillars, changing materials, and optimizing substrate contact surfaces. The reinforcement strategy matching model outputs a fin root reinforcement processing sequence.The fin root strengthening treatment sequence is input into a pre-established heat conduction evaluation model. This heat conduction evaluation model uses the finite element simulation method. After inputting the geometric and material parameters of the strengthening treatment sequence, the model calculates the predicted heat conduction efficiency as the ratio of input heat to output heat through mesh generation, boundary condition setting, and steady-state thermal analysis. When the predicted heat conduction efficiency is greater than 0.85, the combination of the current strengthening treatment sequence and the corresponding fin geometric and material parameters is determined as the overall heat dissipation performance optimization path. This optimization path includes specific fin root thickness adjustment, fin spacing adjustment, and heat conduction pillar arrangement parameters.

[0063] The system generates an updated pressure distribution scheme based on the residual stress 3D distribution data output in step S104. This scheme is a dataset containing residual compressive stress values ​​and their spatial coordinates at different depth levels in various regions of the heat sink, specifically including grid node coordinates, corresponding depth values, and compressive stress levels. The system first identifies high-risk areas, i.e., areas where the residual tensile stress exceeds 0.3 times the material's yield strength or the stress gradient change rate is greater than 50 MPa per millimeter; these areas are prone to fatigue cracking. The system reads the 3D coordinates and compressive stress values ​​of each grid node in the scheme, analyzing them layer by layer from 0 mm at the surface to 5 mm at the depth, extracting the set of nodes with an absolute compressive stress value greater than 200 MPa as the effective processing area. For this area, the system uses a boundary tracing algorithm to identify the contour, connecting points with a spacing of less than 0.5 mm between adjacent nodes to form a closed curve, creating a boundary path for surface hardening treatment. The system arranges this path in descending order of radius of curvature to generate a laser scanning trajectory sequence, with sharp bends with a radius of curvature less than 10 mm designated as key scanning segments. Finally, it outputs a uniform pressure application path file containing start and end coordinates, scanning speed, and energy density parameters.

[0064] The process of extracting the final uniform pressure application path involves complex three-dimensional spatial mapping and trajectory planning algorithms. High-risk areas on the heatsink refer to localized regions identified through finite element analysis where the stress concentration factor is greater than 2.5 and the temperature gradient exceeds 15 degrees Celsius per millimeter. These areas are typically located at the fin root, the connection between the heatsink and the substrate, and at the corners of coolant channels. The system first unfolds the curved geometric features of these high-risk areas into a two-dimensional parametric plane and constructs high-resolution mesh nodes on this plane. The pressure distribution scheme includes three components: the target compressive stress value, the application sequence, and the duration for each mesh node. The initial scheme is generated from a historical process database. The system compares the measured stress field of the current heatsink with the target stress field of the initial scheme, calculates nodes with deviations exceeding 10%, and redistributes their target compressive stress values ​​to obtain an updated pressure distribution scheme. The system calculates the rate of change of stress gradient between adjacent nodes based on the target compressive stress values ​​of each node in the updated scheme. To ensure the uniformity of pressure application and avoid local over-processing or processing dead zones, the system uses an equidistant offset algorithm to generate the initial processing trajectory line, setting the offset distance to 0.8 millimeters, and the number of iterations is adaptively determined based on the area of ​​the region. The system then remaps these two-dimensional trajectory lines back onto the three-dimensional surface model and dynamically adjusts the spacing between the trajectory lines according to the curvature changes. Specifically, when the surface curvature radius R is less than 20 mm, the spacing d between adjacent trajectory lines is calculated using the formula d = 0.5 + 0.02 × R. When the curvature radius R is greater than or equal to 20 mm, the spacing is fixed at 0.9 mm, thus obtaining the final uniform pressure application path. This path clearly defines the movement direction and spatial coordinates of the machining tool on the heat sink surface. The system has a built-in Monte Carlo simulation module. Inputs include the updated pressure distribution scheme, heat sink material parameters, and ambient temperature. The output is the stress distribution statistical results under 1000 random operating conditions. Based on this, the system obtains the Monte Carlo simulation results for the heat sink strengthening process.

[0065] It is understandable that in actual surface shot peening or laser shock peening processes, there are unavoidable random fluctuations in equipment operating conditions, initial material microstructure, and environmental factors. The system introduces Monte Carlo simulation to quantify the potential impact of these uncertainties on the final strengthening quality.

[0066] During the Monte Carlo simulation, the system constructed probability distribution models for three key process parameters. The positioning deviation of the machining tool follows a uniform distribution, ranging from -0.05 mm to +0.05 mm, as this deviation occurs with equal probability in all directions within the equipment's accuracy range. The energy output fluctuation amplitude follows a normal distribution with a mean of 1200 watts and a standard deviation of 48 watts, i.e., a fluctuation rate of 4%, which conforms to the physical characteristics of laser energy output. The material yield strength follows a normal distribution with a mean of 450 MPa and a standard deviation of 22.5 MPa, i.e., a dispersion of 5%. This distribution is obtained by fitting measured data of this batch of materials. The system uses the Latin hypercube sampling method to generate 50,000 samples within the above parameter space, ensuring uniform coverage across all dimensions while maintaining low correlation between samples. Each set of parameters is input into an elastoplastic finite element model based on the von Mises yield criterion. This model receives three process parameters as boundary conditions and material properties, outputting the residual stress tensor components σx, σy, and equivalent stress σeq at the nodes on the heat sink surface. After the calculation, the system statistically processes 50,000 sets of stress results, calculates the arithmetic mean of the equivalent stress as the expected value E(σeq), which is equal to 274 MPa, calculates the standard deviation s(σeq), which is equal to 31 MPa, and calculates the 95% confidence interval based on the normal distribution assumption, from E(σeq) minus 1.96 times s(σeq) to E(σeq) plus 1.96 times s(σeq), which is 213 MPa to 335 MPa. This interval indicates that under the current process conditions, the probability of stress control within this range after strengthening is 95%, thus quantifying the process stability and fault tolerance.

[0067] In one embodiment, the system combines the extracted final uniform pressure application path with the aforementioned Monte Carlo simulation results to generate a fin root reinforcement treatment sequence. The fin root, as the weakest point in the heat sink where thermal stress is most concentrated, requires extremely high precision and customized process arrangements for reinforcement. The system first performs threshold screening on the stress distribution data in the Monte Carlo simulation results, marking areas with stress deviations exceeding 15% of the design value as high-risk deviation areas, and recording their spatial coordinates and deviation amplitude at the fin root. Then, the system compensates for and adjusts the pressure parameters at the corresponding positions in the uniform pressure application path according to the deviation amplitude. The compensation coefficient C is calculated as 0.8 times the ratio of the actual deviation value to the design value, where 0.8 is an empirical safety factor. For areas with deviation amplitudes between 15% and 25%, the pressure increase is controlled within the range of 10% to 18%. For extremely high-risk areas with deviation amplitudes exceeding 25%, the pressure increase is controlled within the range of 18% to 30%, and a secondary reinforcement process is added. The corrected path data was re-verified through stress simulation. Only after the corrected stress deviation was reduced to within 8% could the final fin root strengthening treatment sequence be generated.

[0068] Specifically, the generation of the fin root strengthening treatment sequence aims to clarify the timing arrangement and dynamic parameter adjustment strategy during laser peening. The system employs laser peening to form a uniform residual compressive stress layer on the fin root surface. Based on the process-sensitive areas revealed in Monte Carlo simulations, the final uniform pressure application path is divided into multiple independent processing segments. A specific set of strengthening parameters, including impact angle, single-pulse energy density, and spot overlap rate, is assigned to each segment. These segments are arranged in an ordered sequence, progressing from the fin root center towards the edges. This sequence is based on thermal stress concentration analysis, where the central region experiences 27% to 35% higher cyclic thermal stress amplitude than the edge regions. Prioritizing strengthening the central region maximizes fatigue life improvement. Simultaneously, the system incorporates an energy attenuation compensation mechanism into the sequence. By monitoring the laser output power every five segments processed via an online power meter, when a power drop exceeding 3% is detected, the system increases the pulse energy density of subsequent segments by a corresponding percentage. The compensation amount equals the detected power drop percentage multiplied by a correction factor of 1.08, ensuring a consistent strengthening depth across the entire fin root region. Preferably, the system determines an overall heat dissipation performance optimization path based on the generated fin root strengthening treatment sequence. This optimization path is a complete process plan guiding the heat sink manufacturing process. Unlike the aforementioned final uniform pressure application path, which only describes the pressure distribution of a single processing stage, the optimization path covers the entire process parameter configuration from material selection, fin forming, laser shot peening strengthening to final assembly. The system comprehensively considers the local strengthening process with the global heat conduction requirements of the heat sink. The system analyzes the changes in surface roughness and micromorphology of the fin root surface caused by the strengthening treatment sequence and evaluates the impact of these changes on the fluid boundary layer thickness and convective heat transfer coefficient. Specifically, the three-dimensional morphology data of the strengthened surface are measured using a white light interferometer, and the surface roughness Ra value and micro-pit depth distribution are extracted. The Ra value is substituted into the Prandtl boundary layer thickness correction formula δ=5x / √Rex×1+0.34Ra / x to calculate the change in boundary layer thickness, where x is the distance along the flow direction and Rex is the Reynolds number. Then, based on the change in boundary layer thickness, the rate of change of convective heat transfer coefficient is calculated using the Nusselt number correlation Nu=0.664×Rex^0.5×Pr^0.33×δ0 / δ, where δ0 is the boundary layer thickness of the unstrengthened surface. For example, when determining the overall heat dissipation performance optimization path, the system constructs a multiphysics coupled model incorporating fluid dynamics and solid heat transfer. This model uses the Navier-Stokes equations to describe airflow and the Fourier heat conduction equations to describe heat transfer within the solid. Conjugate heat transfer boundary conditions are used to couple these two processes on the fin surface. The model inputs include wind speed, inlet air temperature, heat source power, and enhanced surface morphology data. The outputs are the overall thermal resistance of the heat sink and a temperature field distribution cloud map. The system uses the surface morphology changes caused by the fin root enhancement treatment sequence as boundary condition inputs to calculate the overall thermal resistance and temperature field distribution of the heat sink under rated operating conditions.The system controls the size and distribution of surface micropits by adjusting the overlap rate of the light spots in the strengthening process sequence. For every 10% increase in the overlap rate, the average diameter of the micropits decreases by 8 to 12 micrometers, and the spacing between micropits decreases by 15 to 20 micrometers. The system sets the overlap rate within the range of 50% to 70% to keep the micropit diameter between 30 and 50 micrometers, which, while improving the fatigue life of the material, can also act as a micro-turbulence generator to disrupt the thermal boundary layer, thereby enhancing local heat transfer efficiency. The system records the combination of process parameters that achieves a balance between structural strength and thermal conductivity. The balance criteria are that the residual compressive stress depth at the fin root reaches more than 0.8 mm and the overall thermal resistance of the heat sink is reduced by no less than 8% compared to the baseline. The system uses a genetic algorithm to search for parameter combinations that simultaneously satisfy these two criteria in the parameter space. After 150 to 200 iterations, the optimal solution set on the Pareto front is selected and used as the overall heat dissipation performance optimization path to guide the final heat sink manufacturing and processing flow.

[0069] Step S107: Through the overall heat dissipation performance optimization path, simulate the fatigue damage evolution under long-term service scenarios, obtain the durability improvement quantitative index, and determine if the quantitative index meets the preset standard, then output the complete reinforcement processing protocol.

[0070] An overall heat dissipation performance optimization path is obtained, which includes temperature field distribution data, heat flux density distribution data, and temperature time history curves of key nodes. Based on the temperature load history in this path, a service load spectrum is extracted. Rainflow counting is used to perform cyclic counting of the load spectrum, and the average stress σm and stress amplitude Δσ of each cycle are statistically analyzed. For the extracted cyclic stress amplitude, the plastic strain increment Δεp of each cycle is calculated using the Ramberg-Osgood material constitutive relation σ=Eε+K(ε)^n, where E is the elastic modulus, K is the strength coefficient, and n is the hardening exponent. The total damage D=Σ(ni / Ni) is calculated using the Miner linear cumulative damage rule, where ni is the actual number of cycles under the i-th stress level, and Ni is the fatigue life under that stress, thus obtaining the fatigue damage evolution trajectory. The damage accumulation rate dD / dt is extracted from the fatigue damage evolution trajectory according to the time step. The damage value D and the number of cycles N are fitted into a power function form D=aN^b using the least squares method, where a and b are fitting parameters. The durability quantification index, i.e., the predicted failure cycle number Nf, is obtained from the fitted curve. If the damage value corresponding to the durability quantification index is greater than the damage tolerance threshold of 0.3, then the corresponding combination of strengthening parameters is matched in the strengthening database. This database contains parameters such as shot peening intensity range of 0.15 to 0.35 mmA, coverage rate of 80% to 98%, and laser power density of 500 to 1500 W / cm², along with their corresponding fatigue life improvement factors. The nearest neighbor combination of strengthening parameters is matched based on the damage value range. The matched shot peening intensity, coverage rate, and laser power are weighted to construct a protocol mapping matrix M, where the matrix element Mij = wi × pi, and wi is the weight coefficient of the i-th parameter, and pi is the normalized parameter value. The protocol mapping matrix is ​​summed column-wise and normalized to extract the observed feature vector O = {o1, o2, o3}, with a vector dimension of 3. The observed feature vectors are input into a pre-established Hidden Markov Model (HMM) for Viterbi decoding. This model defines five hidden states corresponding to five process stages: pretreatment, surface strengthening, heat treatment, coating deposition, and quality inspection. Three observed values ​​correspond to three parameter intensity levels: low, medium, and high. The state transition probability matrix A and the observation probability matrix B are obtained through training with historical strengthening data. The initial state probability π is set to 1 for pretreatment and 0 for the rest. Decoding yields the state transition sequence S = {s1, s2, ..., s5} with the highest probability. A mapping dictionary from state to process is established based on the state transition sequence: state 1 corresponds to sandblasting and rust removal, state 2 to ultrasonic impact strengthening, state 3 to solution aging, state 4 to thermal barrier coating spraying, and state 5 to X-ray residual stress detection. The process parameters corresponding to each state are combined in sequence to output a complete strengthening protocol.

[0071] In one implementation, the system uses the final uniform pressure application path determined in the aforementioned steps as a core component of the overall heat dissipation performance optimization path. This optimization path includes three types of key data: the spatial trajectory coordinate sequence of the ultrasonic impact gun, the impact pressure values ​​corresponding to each trajectory point, and temperature field distribution parameters. The system executes subsequent reliability verification procedures based on this optimization path. Specifically, the system extracts the impact pressure value sequence from the optimization path and inputs it into a fatigue damage evolution model. This model uses the linear cumulative damage rule to calculate the cumulative damage D of the heat sink under cyclic thermal loads, where D equals the sum of each stress cycle ratio ni divided by the corresponding failure cycle number Ni. When D reaches one, fatigue failure is determined. The system sets the service life to 50,000 hours and the thermal cycling frequency to 0.2 Hz, simulating fatigue damage evolution under long-term service scenarios through the optimized path to evaluate the actual contribution of the strengthening process to the heat sink's service life. Simulation results show that the fatigue life of the heat sink treated with the optimized path is improved by more than 30 percentage points compared to the untreated sample.

[0072] Specifically, simulating long-term service scenarios requires the system to construct a dynamic load spectrum that closely matches actual operating conditions. The system collects historical operating data of the heat sink in the target application environment to extract key load features. For alternating thermal stress characteristics, the system uses rainflow counting to identify cycles in the temperature time history, considering cycles with temperature fluctuation amplitudes greater than 5 degrees Celsius as effective thermal load cycles, and statistically obtaining the thermal stress amplitude distribution and cycle number. For mechanical vibration characteristics, the system performs a fast Fourier transform on the vibration acceleration signal, extracting the dominant frequency components with amplitudes greater than 0.5g and frequencies concentrated in the 50-500 Hz range as key vibration loads. The system transforms these discrete load data into a continuous time series function to construct a composite load spectrum including temperature gradient cycles and multi-directional vibration acceleration. During the transformation process, the temperature load uses cubic spline interpolation to generate a continuous temperature change curve T=f(t), with the interpolation node interval set to one-tenth of the actual sampling interval to ensure curve smoothness. The vibration load was reconstructed into a continuous acceleration function A=Σ(An×sin(2πfn×t+φn)) by using the Fourier series superposition method, where An is the amplitude of the nth frequency component, fn is the corresponding frequency, and φn is the phase angle. The first 20 components with a cumulative energy contribution rate of over 95% were retained. This composite load spectrum was used as a boundary condition in the system's mechanical simulation module to reproduce the complex stress state experienced by the heat sink during tens of thousands of hours of service.

[0073] The fatigue damage evolution simulation employs the Smith-Watson-Topper multiaxial fatigue damage model. The model's input is the initial residual compressive stress field formed at the fin root after the aforementioned shot peening treatment. This stress field has a surface compressive stress amplitude of -380 to -420 MPa and an influence depth of 0.15 to 0.25 mm, exhibiting an exponentially decaying distribution. The system performs tensor superposition of this residual compressive stress field with the transient stress field generated by thermal shock and vibration loads during the engine's operating cycle. After superposition, the fin root is in a three-dimensional composite stress state. The system searches for potential physical planes at each grid node at the fin root with a step size of 10 degrees, calculating the maximum shear strain amplitude γmax and normal stress σn on each plane. The Fatemi-Socie criterion is used to determine the critical surface. This criterion is defined as FS equal to γmax multiplied by (1) plus k times σn divided by (σy), where k is a material constant taken as 0.5, and σy is the yield strength. The plane with the largest FS value is the critical surface where fatigue crack initiation is most likely. The system calculates the plastic strain energy density ΔWp for each load cycle on the critical surface using the Ramberg-Osgood cyclic stress-strain relationship. Based on Miner's linear damage accumulation rule, the damage increment ΔD is calculated as 1 divided by Nf, where Nf is the fatigue life cycle number corresponding to that stress level, which can be obtained from the material's ε-N curve. The system accumulates the damage increment for each cycle, and when the accumulated damage D reaches 1, it is considered crack initiation. During the damage accumulation process, the system dynamically updates the local elastic modulus of the material based on the damage variable D, with the update relationship being E = E0 multiplied by (1) minus (D), where E0 is the initial elastic modulus of 210 GPa. Simultaneously, the yield strength σy is updated as σy0 multiplied by (1) minus 0.8 times (D), where σy0 is the initial yield strength of 850 MPa, thus simulating the decay of the material's load-bearing capacity.

[0074] In the damage evolution simulation, the system simultaneously calculates the degradation of the material's microstructure. The system employs the Kocks-Mecking dislocation density evolution model, which describes the change in dislocation density ρ with strain. The evolution equation is dρ divided by dε, which equals k1 multiplied by ρ to the power of 0.5 minus k2 multiplied by ρ. k1 is 0.02, representing the dislocation multiplication coefficient, and k2 is 8 multiplied by 10 to the power of -3, representing the dislocation annihilation coefficient. Under alternating thermodynamic loading, the system calculates the recovery process of a high-density dislocation network. The initial dislocation density is 2.5 x 10^14 per square meter after shot peening, decreasing to 1.8 x 10^14 per square meter after 100,000 thermal cycles. The system employs the Zener-Wert-Avrami equation to calculate residual compressive stress relaxation. The relaxation equation is σr = σr0 multiplied by exp(-t) divided by τ(n), where σr0 is the initial residual compressive stress of -400 MPa, τ is the relaxation time constant (5000 hours at an operating temperature of 850°C), and n is the relaxation exponent (0.6). The system updates the residual compressive stress field every 1000 cycles, feeding the relaxed stress values ​​back to the Smith-Watson-Topper model to correct the stress state in multiaxial fatigue damage calculations, ensuring the accuracy of damage evolution prediction under long-term service conditions. After 50,000 hours of service simulation, the cumulative damage at the critical location at the fin root reached 0.85, with a predicted remaining service life of approximately 9000 hours.

[0075] After completing fatigue damage evolution simulation, the system obtains quantitative indicators for durability improvement. These indicators are a set of specific values ​​used to measure the effectiveness of the strengthening process. First, the system runs a comparative simulation under the same working conditions without residual compressive stress to obtain the baseline fatigue crack initiation life. The system extracts the ratio of the fatigue crack initiation life to the baseline life from the simulation results after strengthening treatment, calculating the life extension rate; a ratio greater than 1 indicates a life extension. Simultaneously, the system analyzes the crack propagation rate in the residual compressive stress layer and calculates the crack propagation hysteresis coefficient by comparing it to the propagation rate of the same crack length in the baseline simulation; a coefficient less than 1 indicates suppressed propagation, with smaller values ​​indicating stronger suppression. The system compares the final fracture cycle count obtained from the simulation with the minimum cycle count required in the heat sink design specifications; the ratio of the difference to the minimum cycle count is defined as the reliability margin index. These indicators together constitute the quantitative indicators for durability improvement, serving as a quantitative basis for evaluating the degree to which a uniform pressure application path improves the durability of the heat sink structure.

[0076] In one embodiment, the system compares the acquired durability improvement quantification indicators with preset standards. The preset standards include the minimum fatigue life threshold of 200,000 cycles required by industry standards and reliability targets customized for specific application scenarios. The system checks each quantification indicator one by one; if all indicators meet or exceed the preset standards, the system determines that the current overall heat dissipation performance optimization path has engineering application value. The overall heat dissipation performance optimization path includes the final uniform pressure application path, the spatial coordinate sequence of the fin root corresponding to the path, and the laser strengthening process parameter set for each coordinate point. The fin root spatial coordinate sequence is automatically extracted and generated by the system based on the uniform pressure application path, and each coordinate point is labeled with three-dimensional Cartesian coordinates Xi, Yi, Zi, where i represents the i-th processing point on the path. The laser strengthening process parameter set includes the impact angle αi, pulse energy density Ei, and spot overlap rate Oi corresponding to the coordinate point. Preferably, if the quantification indicators meet the preset standards, the system outputs a complete strengthening processing protocol. The complete strengthening processing protocol is a standardized document that transforms the theoretical optimization path into actual production instructions. The system uses the ISO 6983 standard to convert spatial coordinate sequences into G-code. Coordinates Xi, Yi, and Zi are directly mapped to X, Y, and Z axis displacement commands in the G-code. The laser power P is calculated based on the pulse energy density Ei and the spot diameter D using the formula P = Ei × π × (D / 2)² × f, where f is the pulse frequency. The system inserts M-commands into the generated G-code to achieve real-time switching of laser parameters. By strictly synchronizing the arrival time ti of each coordinate point with the effective time of the corresponding process parameters, the time deviation is controlled within 10 milliseconds. The system also encapsulates equipment calibration parameters, environmental temperature and humidity control requirements (maintaining 23°C ± 2°C and relative humidity between 45% and 55%), and a sampling frequency of once every 100 processing points for online quality monitoring. This complete enhanced processing protocol not only specifies the kinematic characteristics of the processing tool but also clarifies the dynamic output logic of the laser energy source, providing detailed execution guidelines for the large-scale manufacturing of heat sinks.

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

Claims

1. A method for optimizing the overall performance of a spade-shaped heatsink, characterized in that, include: Geometric model data and historical stress distribution records of the root of the spade-shaped heat sink fins are obtained. Finite element analysis is used to simulate the weakening state of material strength under vibration load and thermal stress field, resulting in a preliminary amplified stress distribution map. Based on this preliminary amplified stress distribution map, a particle swarm optimization algorithm is used to determine the uniform pressure application path, obtaining an optimized pressure distribution scheme. Key node coordinate data are obtained from the optimized pressure distribution scheme, and the crack initiation probability is predicted using Monte Carlo simulation to generate a list of high-risk areas. For the high-risk area list, a target response surface methodology is used to construct a target response surface, obtaining refined strengthening parameters. Based on these refined strengthening parameters, the durability improvement effect is simulated, resulting in an updated pressure distribution scheme. The final uniform pressure application path is obtained from the updated pressure distribution scheme, and combined with the Monte Carlo simulation results, an overall heat dissipation index optimization path is determined. Through this overall heat dissipation index optimization path, fatigue damage evolution is simulated, and a complete strengthening treatment protocol is output.

2. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The process of acquiring geometric model data and historical stress distribution records of the root of the spade-shaped heat sink fins, and simulating the material strength weakening state under vibration load and thermal stress field using finite element analysis to obtain a preliminary stress magnification distribution map includes: extracting the root contour from the geometric model data; performing mapping processing on the geometric parameters and material properties of the root contour with the historical stress distribution records to obtain an initial state matrix; constructing a finite element mesh model based on the initial state matrix; obtaining the finite element mesh elements contained in the finite element mesh model; applying vibration load and thermal stress field to the finite element mesh model; determining the nodal displacement of the finite element mesh elements under the combined action of the vibration load and thermal stress field using the finite element method; calculating the strain gradient of the finite element mesh elements based on the nodal displacements; determining that the finite element mesh elements have undergone fatigue damage if the strain gradient is greater than a preset strain gradient threshold; for the finite element mesh elements that have undergone fatigue damage, obtaining a stiffness degradation coefficient to characterize stiffness degradation; calculating the strength weakening amount of the finite element mesh elements based on the stiffness degradation coefficient; updating the material properties in the initial state matrix based on the strength weakening amount to obtain a degradation state matrix. The updated stress distribution of the root profile is recalculated using the degradation state matrix, and the ratio of the updated stress distribution to the historical stress distribution record is calculated to obtain a preliminary stress magnification distribution map.

3. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The step of determining a uniform pressure application path and obtaining an optimized pressure distribution scheme based on the preliminary amplified stress distribution map using a particle swarm optimization algorithm includes: obtaining stress concentration regions within the target space from the preliminary amplified stress distribution map; extracting tensile features based on the stress direction and magnitude of the stress concentration regions; calculating the spatial curvature of the target space based on the deformation data corresponding to the tensile features; determining the contact area between the target pressure application component and the root of the spade-tooth heat sink fin by matching the spatial curvature with the surface curvature of a preset pressure application component; setting an initial pressure angle based on the contact area and the force direction of the tensile features; calculating the load gradient using the initial pressure angle and a preset initial pressure value; generating at least one discrete pressure application node within the contact area based on the load gradient; constructing a pressure path by connecting the pressure application nodes; using the pressure path as particles and the path coverage rate as a fitness function, performing particle optimization on the pressure path to obtain a path coverage set; Calculate the path coverage rate for each path in the path coverage set. If the path coverage rate is lower than a preset coverage rate threshold, iteratively adjust the initial pressure angle. Recalculate the load gradient based on the iteratively adjusted pressure angle and perform particle optimization to obtain an optimized pressure distribution scheme.

4. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The process of obtaining key node coordinate data from the optimized pressure distribution scheme, predicting crack initiation probability using Monte Carlo simulation, and generating a high-risk area list includes: obtaining fatigue damage indices corresponding to the key node coordinate data; extracting stress alternation amplitude and corresponding load cycle number based on the fatigue damage indices; constructing a stress probability distribution function using the stress alternation amplitude and load cycle number; generating randomly sampled stress samples using the stress probability distribution function; obtaining a preset material fracture toughness; calculating the initial crack length based on the randomly sampled stress samples and the material fracture toughness; calculating the crack propagation rate based on the initial crack length and the randomly sampled stress samples; obtaining a preset number of simulation iterations; using Monte Carlo simulation, iteratively accumulating the crack propagation rate based on the number of simulation iterations to obtain the crack initiation probability corresponding to the key node coordinate data; and determining if the crack initiation probability is greater than a preset probability threshold, then extracting the corresponding key node coordinate data and generating a high-risk area list based on the key node coordinate data.

5. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The method for constructing a target response surface using response surface methodology to obtain refined enhancement parameters for the high-risk area list includes: acquiring region nodes corresponding to the high-risk area list; extracting geometric features and load boundaries based on the region nodes; acquiring preset initial enhancement parameters; and constructing an initial response surface using response surface methodology based on the geometric features, load boundaries, and initial enhancement parameters; acquiring a preset validation sample set; calculating the predicted response value corresponding to the validation sample set using the initial response surface; calculating the fitting error based on the predicted response value and the actual response value of the validation sample set; determining if the fitting error is greater than a preset error threshold; acquiring historical enhancement samples; correcting the initial response surface based on the historical enhancement samples to obtain the target response surface; extracting a parameter sensitivity matrix and response gradient based on the target response surface; calculating a parameter compensation vector using the parameter sensitivity matrix and response gradient; determining the iteration step size based on the parameter compensation vector; and updating the initial enhancement parameters based on the iteration step size and the parameter compensation vector to obtain refined enhancement parameters.

6. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The step of simulating durability improvement effects and obtaining an updated pressure distribution scheme based on the refining and strengthening parameters includes: constructing a fatigue damage assessment model based on the refining and strengthening parameters and a preset load history sequence; calculating the cumulative fatigue damage value through the fatigue damage assessment model; extracting stress concentration areas based on the cumulative fatigue damage value; obtaining preset benchmark test indicators; extracting stress amplification and mitigation indicators for the stress concentration areas; calculating an improvement ratio based on the stress amplification and mitigation indicators and the benchmark test indicators; determining if the improvement ratio is less than a preset improvement threshold; triggering a population backtracking mechanism; obtaining an initial population corresponding to a preset particle swarm optimization algorithm through the population backtracking mechanism; extracting individual extreme values ​​and global extreme values ​​for the initial population; calculating the particle update rate based on the individual extreme values ​​and the global extreme values; performing position mutation processing on the initial population based on the particle update rate to obtain a mutated particle swarm; extracting structural node displacement vectors through finite element simulation using the mutated particle swarm; reconstructing surface contact forces based on the structural node displacement vectors; and obtaining an updated pressure distribution scheme based on the surface contact forces.

7. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The step of obtaining the final uniform pressure application path from the updated pressure distribution scheme and determining the overall heat dissipation index optimization path based on Monte Carlo simulation results includes: extracting surface pressure nodes for the updated pressure distribution scheme, extracting pressure values ​​of adjacent nodes based on the surface pressure nodes, calculating local pressure differences based on the pressure values ​​of adjacent nodes, and constructing the final uniform pressure application path based on the local pressure differences; extracting contact thermal resistance parameters for the final uniform pressure application path, calculating heat flux density and temperature gradient based on the contact thermal resistance parameters, constructing a joint probability distribution function based on the heat flux density and temperature gradient, and performing random sampling on the joint probability distribution function through Monte Carlo simulation to obtain the spatial heat distribution simulation results; Based on the spatial heat distribution simulation results, heat concentration distribution coordinates are extracted. The root region of the spade-shaped heat sink fins is located using these coordinates. Local thermal resistance and target heat transfer cross-sectional area are extracted for this region. The equivalent thermal conductivity is calculated based on these coordinates. If the equivalent thermal conductivity is less than a preset thermal conductivity threshold, a geometric feature vector is extracted for the root region of the spade-shaped heat sink fins. A pre-established enhancement strategy matching model is used to process this geometric feature vector, outputting a spade-shaped heat sink fin enhancement processing sequence. A pre-established heat conduction evaluation model is then used to process this enhancement processing sequence, calculating the predicted heat conduction efficiency. Based on the predicted heat conduction efficiency, an overall heat dissipation index optimization path is determined.

8. The method for optimizing the overall performance of the spade-shaped heat sink as described in claim 1, characterized in that, The process of optimizing the overall heat dissipation index, simulating fatigue damage evolution, and outputting a complete strengthening protocol includes: extracting the service load spectrum based on the overall heat dissipation index optimization path; extracting the cyclic stress amplitude from the service load spectrum; calculating the cumulative plastic strain based on the cyclic stress amplitude and a preset material constitutive relation to obtain the fatigue damage evolution trajectory; extracting the damage accumulation rate from the fatigue damage evolution trajectory; fitting the life decay curve using the least squares method based on the damage accumulation rate to obtain a durability quantification index; determining if the durability quantification index is greater than a preset damage tolerance threshold; matching the corresponding strengthening parameter combination in a preset strengthening database based on the durability quantification index; constructing a protocol mapping matrix based on the strengthening parameter combination; extracting observation feature vectors from the protocol mapping matrix; decoding the observation feature vectors using a pre-established hidden Markov model to obtain a state transition sequence; and outputting a complete strengthening protocol based on the state transition sequence.

9. A spade-shaped heat sink, characterized in that... The spade-shaped heat sink is produced by the method for optimizing the overall performance of the spade-shaped heat sink as described in any one of claims 1-8.

Images