Flange type sliding sleeve material withdrawal structure stress dispersion topology optimization method

Abstract

This invention relates to the fields of computer-aided design and computer-aided engineering manufacturing, specifically to a topology optimization method for stress distribution in a flange-type sliding sleeve ejection structure. Addressing the problems of inaccurate demolding loads, uneven component stress, and the inability to directly machine the model in traditional optimization methods, this invention extracts dynamic frictional resistance and thermal shrinkage stress fields through thermo-mechanical coupling analysis and maps them to time-varying boundary loads. It establishes a variable-density topology optimization model with the goal of minimizing the stress squared difference, setting the lower punch and intermediate ejector as a dual-solid design domain. Manufacturing constraints are incorporated into the iteration process to drive adaptive material redistribution and calculate the optimal load-sharing ratio. Finally, isosurfaces are extracted and fitted as standard geometric primitives, and B-spline curves are used to smooth the stress concentration areas with chamfer transitions, outputting a reconstructed solid. This invention achieves coordinated and balanced stress distribution in the ejection structure, improves mold life, and generates a solid model that can be directly CNC machined.

Description

Technical Field

[0001] This invention relates to the fields of computer-aided design, computer-aided engineering, and computer-aided manufacturing, and specifically to a method for topology optimization of stress distribution in a flange-type sliding sleeve ejection structure. Background Technology

[0002] Flange-type sliding sleeves are key components in fields such as new energy vehicle transmission systems. Their production often employs warm forging to achieve excellent mechanical properties and dimensional accuracy. After warm forging, the high-temperature sleeve needs to be ejected from the mold through a ejector mechanism. Due to the large flange diameter and thin shell of these flange-type sleeves, the design space for the ejector pins in the ejector mechanism is severely limited, making them prone to fatigue fracture under harsh demolding stress. Therefore, optimizing the ejector structure using computer-aided design methods is crucial for ensuring mold life and improving production efficiency.

[0003] In existing technologies, the design of ejection structures largely relies on the experience of design engineers. When using computer-aided design to optimize ejection structures, several significant technical shortcomings exist: First, the input load deviates from the actual working conditions. Traditional stress analysis typically applies a static, uniformly distributed thrust to the ejection structure based on a simplified mechanical model for strength verification or topology optimization of a single part. This estimation is often based on empirical formulas or selects only a maximum instantaneous force during the entire demolding process as a constant load, failing to fully consider the spatial dynamic frictional resistance field and thermal contraction stress field generated by the combined effects of thermal expansion and contraction and material sealing extrusion during warm forging. The inability to accurately reproduce this complex time-varying load field leads to inaccurate load input in the optimization design, affecting the reliability of the final structural design. Second, the isolated design of assembled components prevents optimal stress distribution. For assembled ejection structures requiring coordinated operation, such as the peripheral punch with an opening and the intermediate ejector, although engineering practices employ the method of adding an intermediate ejector with an opening in the punch, the design process is often separated and carried out independently. The load-bearing distribution ratio among components relies heavily on engineers' experience and trial-and-error, making it impossible to accurately calculate the stress relationship between the lower punch and the intermediate ejector. This non-cooperative optimization fails to achieve a reasonable distribution of local stress from the slender ejector rod to the outer mold, easily leading to excessive stress on a certain component, making it a weak point, resulting in material waste or premature mold damage. Third, traditional optimization models cannot be directly implemented in manufacturing. In existing technologies, the 3D models generated by pure topology optimization algorithms often have irregular surfaces, with numerous rough patches, sharp corners, and unmachineable internal holes. These models lack engineering geometric constraints and smooth transitions, making it difficult to directly import them into CNC machine tools (CAM software) for actual manufacturing. This means that the results of topology optimization often remain theoretical and are difficult to truly apply to mold production. Summary of the Invention

[0004] In view of this, in order to solve the problems mentioned in the background technology, a topology optimization method for stress distribution of flange-type sliding sleeve ejection structure is proposed.

[0005] The objective of this invention can be achieved through the following technical solution: This invention provides a topology optimization method for stress dispersion of flange-type sliding sleeve ejection structure, including: S1, constructing an initial three-dimensional finite element model of flange-type sliding sleeve and mold, obtaining physical parameters in the warm forging process and inputting them into the initial three-dimensional finite element model, performing transient thermo-mechanical coupling analysis, and extracting the spatial dynamic frictional resistance field and thermal shrinkage stress field.

[0006] S2. Map the spatial dynamic frictional resistance field and thermal shrinkage stress field to the time-varying boundary load conditions of the topology optimization design domain of the material ejection structure, and generate a dynamic load mapping model.

[0007] S3. Set the outer opening punch and the middle ejector as a dual-solid design domain for collaborative topology optimization. Define the contact surface as the sliding friction boundary. Based on the variable density method and dynamic load mapping model, establish a topology optimization model with the objective function of minimizing the maximum equivalent stress square difference between the two entities and the volume retention rate as the constraint condition.

[0008] S4. During the iterative optimization process of the topology optimization model, the material density field is adaptively redistributed between the intermediate ejector and the outer hole punch wall to calculate the optimal load sharing ratio and embed the processing and manufacturing geometric constraints to generate a three-dimensional mesh density field.

[0009] S5. Extract the isosurfaces of a given density threshold in the 3D mesh density field, and constrain the isosurfaces to geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm.

[0010] S6. Identify stress concentration regions in the isosurface, use B-spline curves to perform smooth chamfer transition calculations with continuous curvature, generate smooth reconstructed solid data and output it.

[0011] Compared with the prior art, the embodiments of the present invention have at least the following advantages or beneficial effects: (1) The present invention accurately extracts the spatial dynamic frictional resistance field and thermal shrinkage stress field generated by thermal expansion and contraction and sealing extrusion during warm forging through transient thermo-mechanical coupling analysis, and maps them as time-varying boundary loads. This technical means replaces the traditional static uniformly distributed thrust, accurately reproduces the complex time-varying load field of the entire demolding process, improves the accuracy of load input, and thus provides a real physical basis for topology optimization, ensuring the reliability of the final structural design.

[0012] (2) In this invention, the outer-perforated lower punch and the intermediate ejector are set as a co-optimized dual-solid design domain, and adaptive redistribution is performed with the goal of minimizing the maximum equivalent stress variance between the two entities. This method accurately calculates and solves the optimal load sharing ratio between the lower punch and the intermediate ejector, and reasonably disperses the local stress concentration on the slender ejector rod to the outer mold, realizing the balanced load bearing of each component of the ejector structure, effectively avoiding premature fatigue fracture of a single component, and improving the overall service life of the mold.

[0013] (3) This invention incorporates manufacturing geometric constraints during the iteration process and forcibly fits the isosurface into geometric primitives such as stepped cylindrical surfaces, conical surfaces, or planes through feature recognition. At the same time, it uses B-spline curves to perform smooth chamfer transition processing on stress concentration areas to ensure continuous curvature. This series of operations not only eliminates local high-stress points but also transforms the abstract mesh density field into standard solid data with smooth surfaces and clear geometric features, ensuring that the generated model can be directly imported into CNC machine tools for processing, thus breaking down the technical barriers from theoretical optimization to manufacturing implementation.

[0014] (4) This invention deeply integrates complex finite element physical simulation, time-varying load mapping, collaborative topology optimization, geometric feature fitting and solid smooth reconstruction to form a complete method and system closed loop. It breaks away from the traditional mode of relying heavily on the personal experience and trial and error of engineers, and realizes the full-process automation and intelligent design of flange-type sliding sleeve ejection structure from load analysis to the final machinable model output, shortening the mold development cycle and improving design efficiency. Attached Figure Description

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

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

[0017] Figure 2 This is a diagram illustrating the evolution of the transient thermo-mechanical coupling multiphysics field during demolding, as described in this invention.

[0018] Figure 3 This is the phase diagram of the iterative convergence of the collaborative topology optimization of the present invention.

[0019] Figure 4 This is a comparison diagram of the 2D thermodynamic gradient of the stress concentration region B spline smooth reconstruction according to the present invention. Detailed Implementation

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

[0021] Please see Figure 1 This invention provides a topology optimization method for stress distribution in a flange-type sliding sleeve ejection structure, including: (See reference) Figure 2 S1. Construct an initial three-dimensional finite element model of the flange-type sliding sleeve and the mold, obtain the physical parameters in the warm forging process and input them into the initial three-dimensional finite element model, perform transient thermo-mechanical coupling analysis, and extract the spatial dynamic frictional resistance field and thermal shrinkage stress field.

[0022] In a specific embodiment of the present invention, an initial three-dimensional finite element model of the flange-type sliding sleeve and the mold is constructed, the physical parameters in the warm forging process are obtained and input into the initial three-dimensional finite element model, transient thermo-mechanical coupling analysis is performed, and the spatial dynamic frictional resistance field and thermal shrinkage stress field are extracted, including: obtaining the geometric dimension data of the flange-type sliding sleeve and the structural parameters of the mold, and constructing the initial three-dimensional finite element model.

[0023] The heating temperature, upsetting and flange pressing tonnage, and temperature control cooling rate in the warm forging process are collected as physical parameters.

[0024] Input the physical parameters into the initial three-dimensional finite element model, run the transient thermo-mechanical coupling analysis algorithm, and simulate the instant of demolding of the sliding sleeve and the continuous ejection stage.

[0025] Extract the spatial dynamic frictional resistance field and thermal contraction stress field from the simulation results.

[0026] Specifically, this step aims to accurately reproduce the physical process of the flange-type sliding sleeve being ejected from the mold after warm forging through computer simulation, thereby obtaining the load data necessary for subsequent structural optimization.

[0027] The implementation process begins with building a digital simulation model, namely the initial three-dimensional finite element model. This model is a virtual representation of the precise geometric and physical properties of the flange sleeve and mold in the physical world. First, it is necessary to obtain the geometric dimensional data of the flange sleeve, including its outer diameter, inner diameter, height, flange diameter, and thickness, as well as the structural parameters of the mold, such as the contour dimensions of the mold cavity and the position and dimensions of the ejection mechanism. This data typically originates from the product's design drawings. Using this data, precise three-dimensional geometry of the flange sleeve and mold is created in 3D modeling software. This geometry is then imported into finite element analysis software, where it is discretized into a set of numerous tiny elements, each with independent physical properties, using mesh generation technology, thus completing the construction of the initial three-dimensional finite element model.

[0028] Next, to ensure the simulation closely resembles real-world conditions, key physical parameters from the warm forging process need to be collected. Warm forging is a plastic forming process performed below the metal recrystallization temperature. The physical parameters collected in this step comprise three core components. The first is the heating temperature, specifically the initial temperature to which the flange sleeve is heated before forging, which determines the initial thermal state and mechanical properties of the material. The second is the upsetting and flange clamping tonnage, which refers to the main forming force applied to the workpiece during forging, usually measured in tons, and needs to be converted to forces (Newtons) in the International System of Units (SI) for calculations. The conversion relationship is shown in the following formula: ; in, Represents the applied force, measured in Newtons (N). This represents the tonnage of upsetting and flange pressing, and the unit is kilogram (kg), which is converted from tonnage (1 ton = 1000 kg). This is the acceleration due to gravity, usually taken as 9.8 m / s².

[0029] The third is the temperature control cooling rate, which refers to how quickly the workpiece's temperature drops over time during demolding. This directly affects the workpiece's shrinkage behavior and the resulting internal stress. These physical parameters are set based on actual measurements taken on the production floor, for example, by statistical analysis of data collected from 200 sets of industrial sensors on an actual production line.

[0030] The collected physical parameters are input into the constructed initial three-dimensional finite element model as boundary and initial conditions for the simulation calculation. For example, the heating temperature is set as the initial temperature of all nodes of the flange-type sleeve in the model; the calculated forming force is applied to the corresponding action surface of the mold; and the temperature-controlled cooling rate is defined as the heat exchange condition between the sleeve and the surrounding environment. Subsequently, the transient thermo-mechanical coupling analysis algorithm is run. Transient thermo-mechanical coupling analysis is an advanced numerical simulation method that can simultaneously calculate the temperature and force fields of the model as they change over time, and consider the interaction between the two, i.e., thermal expansion and contraction generate stress, while plastic deformation and friction generate heat. This analysis focuses on simulating the instant of sleeve demolding and the continuous ejection stage, which are critical periods that determine the magnitude of the demolding force and the stress distribution.

[0031] After the simulation calculations were completed, two core physical fields were extracted from the massive dataset. The first is the thermal shrinkage stress field. Because the temperature of the flange sleeve is higher than that of the mold, during the cooling process, the shrinkage of the sleeve is greater than that of the mold. This mismatch in shrinkage causes the sleeve to tightly "hug" the mold core, thus generating a huge positive pressure at the contact interface, i.e., thermal shrinkage stress. Its magnitude can be explained in principle by the following formula: ; in, This represents thermal shrinkage stress, and the unit is Pascal (Pa). It is the elastic modulus of a material at the current temperature, and its unit is Pascal (Pa). It is the average linear expansion coefficient of the material in the current temperature range, and the unit is per degree Celsius (1 / °C). It is the effective temperature difference between the flange-type sliding sleeve and the mold, expressed in degrees Celsius (°C).

[0032] The second is the spatial dynamic frictional resistance field. This resistance arises from the relative sliding between the sleeve and the mold during the ejection process. Its magnitude mainly depends on the normal force generated by the thermal shrinkage stress field and the coefficient of friction between them. The frictional resistance can be expressed as: ; in, Represents frictional resistance (expressed in the form of stress), with the unit being Pascal (Pa). It is a dynamic friction coefficient that is related to temperature, contact pressure and relative sliding speed. It is a dimensionless parameter and is usually taken as 0.15. It is the normal stress component generated on the contact surface by the thermal shrinkage stress field, and its unit is Pascal (Pa).

[0033] Both of these fields are "spatially dynamic," meaning that the stress and friction values ​​at different locations on the sliding sleeve change continuously at different times during the demolding process. By using finite element post-processing tools to export the normal contact stress and friction shear stress data of all nodes on the contact interface at each time step, the required thermal shrinkage stress field and spatial dynamic friction resistance field are obtained.

[0034] For example, to more clearly illustrate the implementation process of this step, a specific example is given below. Assume the geometric dimensions of a flange-type sliding sleeve are: outer diameter 100 mm, height 80 mm, and flange diameter 150 mm. The structural parameters of the mold are matched accordingly. The flange-type sliding sleeve is made of 42CrMo steel, and the mold is made of H13 steel. Based on these data, an initial three-dimensional finite element model containing 50,000 tetrahedral elements is constructed.

[0035] The physical parameters collected and determined from the production line are as follows: heating temperature is 850°C; upsetting and flange pressing tonnage is 500 tons; temperature-controlled cooling rate is 10°C / s. First, convert the tonnage into the force required for calculation. 500 tons equals 500,000 kilograms. .

[0036] The initial condition of heating temperature 850°C was used as the initial condition for flange-type sliding sleeve, the force of 4,900,000 Newtons was used as the load boundary condition, and the cooling rate of 10°C / s was used as the heat flow boundary condition. These conditions were then input into the initial three-dimensional finite element model to run transient thermo-mechanical coupling analysis.

[0037] After the analysis was completed, the calculation results for a node A on the inner wall of the sliding sleeve were extracted at 0.5 seconds into the demolding process. The simulation results showed that the temperature of node A had dropped from 850°C to 750°C, while the surface temperature of the mold core it contacted was 250°C. The finite element software, taking into account material properties, temperature changes, and geometric constraints, calculated the component of the thermal shrinkage stress field of this node in the normal direction of the contact surface, i.e., the normal stress. It is 150 megapascals (MPa).

[0038] Meanwhile, according to experimental data, the dynamic friction coefficient between 42CrMo steel and H13 steel at 750°C... The value is 0.15. Using the normal stress at this node, the specific value of the spatial dynamic frictional resistance field at that point is calculated: .

[0039] This means that at 0.5 seconds into the demolding process, node A experiences a normal compressive stress of 150 MPa and a frictional resistance of 22.5 MPa in the ejection direction. By repeating this calculation and data extraction at every time step throughout the demolding process for all interface nodes in contact between the sleeve and the mold, a complete spatial dynamic frictional resistance field and thermal shrinkage stress field, varying with time and space, consisting of a series of data points, are finally obtained.

[0040] S2. Map the spatial dynamic frictional resistance field and thermal shrinkage stress field to the time-varying boundary load conditions of the topology optimization design domain of the material ejection structure, and generate a dynamic load mapping model.

[0041] In a specific embodiment of the present invention, the spatial dynamic frictional resistance field and thermal shrinkage stress field are mapped to the time-varying boundary load conditions of the topology optimization design domain of the material removal structure to generate a dynamic load mapping model, including: extracting the boundary node coordinates of the topology optimization design domain of the material removal structure.

[0042] Based on the boundary node coordinates, the spatial dynamic frictional resistance field and the thermal contraction stress field are spatially interpolated to obtain the node load vector.

[0043] The nodal load vectors are serialized according to the time step to construct time-varying boundary load conditions.

[0044] Time-varying boundary load conditions are applied to the topology optimization design domain of the material removal structure to generate a dynamic load mapping model.

[0045] Specifically, the core task of this step is to accurately convert and apply the complex load data obtained in step S1, which describes the physical process of the sleeve demolding, to the material removal structure model used for subsequent structural optimization. This conversion process ensures that the mechanical environment on which the optimization design is based can truly reflect the harsh working conditions.

[0046] The process begins by identifying the target object to which the load is applied, namely the topology optimization design domain of the ejection structure. This design domain is a three-dimensional digital model that will be used for structural morphology optimization in subsequent steps. It typically represents components directly involved in the ejection process, such as the intermediate ejector and the peripheral punch with openings in the mold. To accurately apply the load to its surface, the spatial coordinates of all nodes located on its outer boundary must first be extracted from the finite element mesh of the model. This set of coordinates constitutes the boundary node coordinates.

[0047] After obtaining the target locations, the next step is to calculate the load values ​​that these locations should bear. Since the spatial dynamic frictional resistance field and thermal shrinkage stress field obtained in step S1 are based on the mesh nodes of the flange-type sliding sleeve and mold contact surface, the spatial locations of these nodes usually do not completely coincide with the boundary node coordinates of the topology optimization design domain of the ejection structure. Therefore, spatial interpolation calculation is required. This calculation, for each boundary node, finds several nearest load data points from the spatial dynamic frictional resistance field and thermal shrinkage stress field in its surrounding space, and then calculates the load value that the boundary node should bear using a weighted average based on the location and load value of these data points. A commonly used interpolation method is the inverse distance weighted method, and its calculation formula is as follows: ; in, At a certain point in time At that time, a certain boundary node on the topology optimization design domain of the material return structure The load value to be determined is a nodal load vector containing normal and tangential components, with units of Pascal (Pa) or Megapascal (MPa). At the same time point The distance to the node obtained from the spatial dynamic frictional resistance field or thermal contraction stress field. The most recent The load values ​​of a known load data point are expressed in Pascals (Pa) or megapascals (MPa). It is a node With the The Euclidean distance between known load data points is expressed in meters (m) or millimeters (mm). This is the number of known load data points used in the interpolation calculation. The setting is based on ensuring calculation accuracy while taking into account calculation efficiency, and is usually set to 4 to 8 points based on experience. It is a power exponent used to adjust the influence of distance weights, typically taking a value of 2. By performing this calculation on all boundary nodes, the set of nodal load vectors of the entire design domain boundary at a specific time is obtained.

[0048] Since the demolding process is dynamic, the load varies with time. Therefore, the spatial interpolation calculation described above needs to be repeated at all discrete time steps simulated in step S1. A time step is the time interval for recording results in transient analysis. Arranging the set of nodal load vectors calculated at each time step in chronological order constructs the time-varying boundary load condition. This is essentially a four-dimensional dataset containing three-dimensional spatial coordinates and one-dimensional time-based load distribution information.

[0049] The final step involves applying the constructed time-varying boundary load conditions to the corresponding boundary nodes of the topology optimization design domain of the ejected structure using the load application function of the finite element analysis software. After loading, the topology optimization design domain of the ejected structure carries complete, time-varying dynamic load information, transforming into a dynamic load mapping model ready for subsequent analysis and optimization.

[0050] For example, continuing with the example from step S1, assume that spatial dynamic frictional resistance field and thermal contraction stress field data, recorded every 0.1 seconds from t=0 to t=1, are obtained from transient thermo-mechanical coupling analysis. Now, these loads need to be mapped onto the topology optimization design domain of the ejector structure. First, extract the coordinates of all boundary nodes on the surface of this design domain that contacts the sliding sleeve. Select one of the boundary nodes P, whose coordinates are (10.0, 20.0, 80.0), in millimeters.

[0051] Next, taking a time step of t=0.5 seconds as an example, spatial interpolation is performed on node P to obtain its nodal load vector. In the spatial dynamic frictional resistance field data, the three data points A, B, and C closest to node P are found, and their coordinates and frictional resistance at t=0.5 seconds are as follows: Data point A: coordinates (10.1, 20.0, 80.0), frictional resistance... Megapascals. Data point B: coordinates (9.9, 20.1, 80.0), frictional resistance. Megapascals. Data point C: coordinates (10.0, 19.9, 79.9), frictional resistance. Megapascal.

[0052] Calculate the distances from node P to A, B, and C: millimeters millimeters Millimeters.

[0053] Using the inverse distance weighting method, the power exponent is set. ,Pick ,and Calculate the frictional resistance at node P at t=0.5 seconds. : Megapascal.

[0054] The same method can be used to calculate the normal load generated by thermal shrinkage stress. Combining these two components yields the nodal load vector at node P at t=0.5 seconds. This interpolation process is applied to all boundary nodes and repeated for all 11 time steps of t=0, 0.1, 0.2, ..., 1.0 seconds to construct the complete time-varying boundary load conditions.

[0055] Finally, the time-varying boundary load condition, which contains 11 time-step load sequences, is applied to the topology optimization design domain of the unloaded structure, thus generating a dynamic load mapping model, which provides accurate input for subsequent topology optimization analysis.

[0056] See Figure 3 S3. Set the outer opening punch and the middle ejector as a dual-solid design domain for collaborative topology optimization. Define the contact surface as the sliding friction boundary. Based on the variable density method and dynamic load mapping model, establish a topology optimization model with the objective function of minimizing the maximum equivalent stress square difference between the two entities and the constraint condition of volume retention rate.

[0057] In a specific embodiment of the present invention, the outer opening punch and the intermediate ejector are set as a dual-solid design domain for collaborative topology optimization. The contact surface is defined as a sliding friction boundary. A topology optimization model is established based on the variable density method combined with the dynamic load mapping model, with the minimization of the maximum equivalent stress square difference between the two entities as the objective function and the volume retention rate as the constraint condition. This includes: dividing the outer opening punch and the intermediate ejector into mesh cells and setting them as dual-solid design domains.

[0058] Extract the interface between the two solid design domains and define it as a sliding friction boundary.

[0059] Based on the variable density method, the relative density of the mesh elements is used as the design variable, and the equivalent stress of each element is calculated by combining the dynamic load mapping model.

[0060] Calculate the maximum equivalent stress of the outer opening punch and the middle ejector, and construct an objective function to minimize the square difference of the maximum equivalent stress of the two entities.

[0061] Set a threshold for the ratio of total volume to initial volume, construct volume retention rate constraints, and generate a topology optimization model.

[0062] Specifically, the purpose of this step is to establish a mathematical optimization model to guide subsequent structural innovation design. This model will determine how the two working components, the peripheral punch and the intermediate ejector, distribute material optimally to jointly withstand the severe loads from the demolding process.

[0063] First, the two independent components—the peripheral punch and the intermediate ejector—are defined as a unified optimization object in the digital environment, i.e., a two-body design domain. This means that in subsequent optimization calculations, they will be treated as an interconnected system, rather than two isolated parts. The three-dimensional geometric models of these two components are meshed, discretized into thousands of tiny finite element elements, which forms the basis for numerical calculations.

[0064] During the unloading process, there is relative motion and force transmission between the outer punch with the outer opening and the middle ejector. Therefore, it is necessary to precisely define their interaction. The surfaces where these two entities geometrically contact each other, i.e., the interface, are extracted and defined as the sliding friction boundary. This boundary condition allows for frictional relative sliding between the two surfaces in the tangential direction and the transmission of pressure in the normal direction. Its frictional behavior follows Coulomb's law of friction, requiring the setting of a friction coefficient. This coefficient is based on typical frictional experimental data of die steel under high temperature and high pressure.

[0065] Next, the variable density method is adopted as the core algorithm for topology optimization. The variable density method is a clever material placement approach that assigns a variable property, namely relative density, to each mesh element in the two-solid design domain. This relative density is the design variable in the subsequent optimization process, and its value ranges from 0 to 1. When the relative density of an element is 1, it represents a solid material; when it is 0, it represents a hole that should be removed; values ​​between 0 and 1 represent a lower-density transitional material. To make the calculation results tend towards clear solids or holes, the mechanical properties of the elements, such as the elastic modulus, are nonlinearly related to their relative density, typically represented using the Solid Isotropic Material Penalty Model (SIMP). ; in, It is the first The actual elastic modulus of each unit; It is the elastic modulus of the original solid material; It is the first The relative density of each unit, i.e., the design variable; This is a penalty factor, typically set to 3, to ensure that the relative density in the optimization result is close to 0 or 1, thus obtaining a clear structural profile. After setting this relationship, combined with the time-varying boundary load conditions provided by the dynamic load mapping model generated in step S2, a finite element analysis is performed on the entire two-body design domain, and the equivalent stress borne by each element during the entire demolding process can be calculated.

[0066] The optimization goal is to ensure that the outer punch with the outer opening and the intermediate ejector can work together to bear the load, preventing one component from failing prematurely due to excessive load. To this end, an objective function is constructed that minimizes the square difference of the maximum equivalent stresses of the two components. First, during the entire dynamic load application time, the maximum equivalent stress in all elements of the outer punch with the outer opening and the maximum equivalent stress in all elements of the intermediate ejector are identified. Then, the square difference between these two maximum stress values ​​is calculated. The smaller the square difference, the closer the peak stress levels of the two components are, and the more balanced the load bearing. The objective function can be expressed as: ; in, It is the objective function value to be minimized, and its dimension is the square of the stress. It is the maximum equivalent stress value found globally in the punch with the outer opening; This is the maximum equivalent stress value found through a global search within the intermediate unloader. The optimization algorithm will adjust the relative density of all elements. To keep trying to reduce The value of .

[0067] Finally, to prevent the optimization algorithm from removing all material in pursuit of zero stress, and also to control the weight and cost of the final product, a volume retention constraint needs to be set. This constraint limits the total amount of material used in the optimized structure. A threshold ratio of total volume to initial volume is set, requiring that the sum of the volumes of all elements (considering relative density) must not exceed a specific percentage of the initial total volume. This constraint is expressed as: ; in, and They are the first The relative density and volume of each unit; It is the total number of elements in the dual-entity design domain; It is the total volume of the initial design domain; This is a set threshold for the volume retention ratio, for example, 0.5, meaning the final structure's volume does not exceed 50% of the initial volume. This threshold is typically set based on lightweight design specifications or cost budgets for the product. Combining the objective function with the constraints formally constructs the topology optimization model for subsequent solutions.

[0068] For example, continuing the calculation steps above, a topology optimization model is now established for the peripheral punch and the intermediate ejector. First, the geometric models of these two components are imported into the finite element software and divided into 100,000 mesh elements, which are collectively set as a two-solid design domain. The interface between the two is extracted and defined as the sliding friction boundary, and the coefficient of dynamic friction between them is set to 0.1 according to the material handbook.

[0069] The variable density method was used, and the material was H13 steel. Its elastic modulus was... The penalty factor is 190 gigapascals (GPa). Set it to 3. Apply the dynamic load mapping model generated in step S2 to this two-body design domain.

[0070] In one optimization iteration, the algorithm is assumed to provide a relative density distribution scheme. Based on this scheme, finite element analysis is performed to calculate the maximum equivalent stress of the lower punch with the outer opening during the entire demolding process. The maximum equivalent stress of the intermediate feeder is 450 megapascals (MPa). It is 550 megapascals (MPa).

[0071] Based on these values, calculate the objective function value for the current iteration: This non-zero objective function value indicates that the load on the two components is not yet balanced.

[0072] Simultaneously, volume retention rate constraints are set. Based on lightweight requirements, a volume retention rate ratio threshold is set. The value is 0.6. Calculate the relative density of all elements in the current iteration. and volume The obtained total volume is then divided by the initial total volume. The current volume retention rate is 0.58. Because... Therefore, the volume retention rate constraint is satisfied.

[0073] Through the above steps, a complete topology optimization model, including design variables, objective function, and constraints, has been constructed, preparing for the next step of iterative optimization. The optimization algorithm will then proceed in subsequent steps based on the current... Based on the objective function value and the constraint conditions, the material density distribution is adjusted in order to obtain a smaller objective function value in the next iteration.

[0074] S4. During the iterative optimization process of the topology optimization model, the material density field is adaptively redistributed between the intermediate ejector and the outer hole punch wall to calculate the optimal load sharing ratio and embed the processing and manufacturing geometric constraints to generate a three-dimensional mesh density field.

[0075] In a specific embodiment of the present invention, during the iterative optimization process of the topology optimization model, the material density field is adaptively redistributed between the intermediate ejector and the outer opening punch wall to calculate the optimal load sharing ratio and embed manufacturing geometric constraints to generate a three-dimensional mesh density field, including: calculating the sensitivity of the objective function of the topology optimization model with respect to the design variables.

[0076] The material density field is updated based on sensitivity to achieve adaptive redistribution.

[0077] Based on the updated material density field, the bearing capacity of the intermediate ejector and the outer opening punch wall is calculated, and the optimal load sharing ratio is determined.

[0078] During the density update process, the draft angle along the axial direction is forced to be greater than zero and the generation of internal holes is restricted as a geometric constraint for processing and manufacturing.

[0079] After multiple iterations until the objective function converges, a three-dimensional mesh density field is output.

[0080] Specifically, this step is the core of the topology optimization process. Through a series of iterative calculations, it automatically explores and determines the optimal distribution of materials within the design space, thereby achieving the preset optimization objective. This process is based on the topology optimization model established in step S3 and is a closed-loop dynamic optimization process.

[0081] The driving force behind the iterative optimization process comes from sensitivity analysis. First, it is necessary to calculate the sensitivity of the objective function of the topology optimization model with respect to each design variable. In this embodiment, the objective function is the maximum squared difference of equivalent stress between the two entities. The design variable is the relative density of each grid cell. Sensitivity is essentially the partial derivative of the objective function with respect to design variables. It measures how much the objective function changes when the relative density of a given element changes slightly. Its mathematical expression is: .

[0082] The sign and magnitude of this value reveal the adjustment unit. The density of a cell contributes to the optimization objective in terms of direction and extent. A large negative value means that increasing the density of the cell can significantly reduce the maximum stress difference between the two solids, thus making the cell a "favorable" location for adding material. Conversely, a positive value or a value close to zero means that increasing the density of the cell has little or no improvement on the objective, making it a "candidate" location for removing material.

[0083] Based on the calculated sensitivity, the optimization algorithm updates the entire material density field, achieving an adaptive redistribution of material between the intermediate ejector and the outer punch wall. This update process is not a simple linear increase or decrease, but is controlled by complex optimization algorithms, such as the Optimality-Criteria method or the Method of Moving Asymptotes. A simplified update logic can be understood as follows: ; in, It is the first Unit during the next iteration The new relative density; It is the first The relative density of the next iteration; This represents the update function of the optimization algorithm. Under the premise of satisfying the volume retention rate constraint, this function removes limited material resources from regions with low sensitivity (small contribution) and redistributes them to regions with high sensitivity (large contribution).

[0084] During the density update cycle, to ensure the final structure is practically manufacturable, manufacturing geometry constraints must be incorporated. For example, for parts requiring casting or forging, demolding feasibility must be considered. Therefore, a draft angle greater than zero along the axial direction is mandatory, meaning all vertical or near-vertical structural sidewalls must have a slight tilt angle to avoid interference or scratches during demolding. Simultaneously, to prevent the formation of closed cavities that cannot be removed by subsequent machining, the generation of internal cavities is limited. These constraints are achieved by applying specific density filtering or geometric projection techniques after density updates. Specifically, by introducing a geometric projection model based on the Heaviside step function, design variables are mapped to a physical density field with a specified tilt angle to achieve draft angle constraints. At the same time, the Virtual-Temperature-Method is adopted, which sets the low-density pore region to be self-heating and the high-density solid region to have high thermal conductivity. The excessively high virtual temperature generated in the closed pores is added as a penalty term to the optimization objective, thereby limiting and eliminating the low-density pore region completely covered by the solid material in the iteration, thus correcting the density distribution that does not meet manufacturing requirements.

[0085] The optimization process ends when the iterative process converges, meaning the objective function value no longer decreases significantly and the material density field stabilizes. The resulting material distribution is considered the optimal solution under the current constraints. Based on this final, stable material density field, the optimal load-sharing ratio can be calculated. By performing a static analysis on the final structure, the load-bearing capacity borne by the intermediate ejector and the outer opening punch wall under the total load is calculated, typically in the form of structural reactions. The optimal load-sharing ratio can be expressed as: ; in, It is the load-sharing ratio of the punch wall under the outer opening; and These figures represent the load-bearing capacities of the outer punch wall and the intermediate ejector, calculated using finite element analysis, in Newtons (N). This ratio visually reflects the synergistic load-bearing effect achieved by the optimized design. The iterative optimization process terminates after multiple iterations until the objective function converges. The final output is a dataset containing the relative density values ​​of all mesh elements in the two-entity design domain, i.e., a three-dimensional mesh density field. This field digitally depicts the macroscopic morphology and internal topology of the optimal structure.

[0086] For example, following the example from step S3, at the beginning of a certain iteration, the objective function... for First, calculate the objective function. The sensitivity is given for the relative density of all 100,000 elements in the two-body design domain. Assuming that for a specific element A in the low-stress region of the outer-aperture punch, the calculated sensitivity is... For a specific unit B in the high-stress region of the intermediate unloader, its sensitivity is... .

[0087] Based on the sensitivity, the optimization algorithm adaptively redistributes the material density field. Since element A has a positive sensitivity, increasing its density will increase the stress difference; therefore, the algorithm tends to decrease its relative density, for example, updating it from 0.8 to 0.7. Element B, on the other hand, has a large negative sensitivity, indicating that increasing its density effectively reduces the stress difference; the algorithm will increase its relative density, for example, updating it from 0.3 to 0.4. During the update process, all structural sidewalls are simultaneously checked and forced to meet the manufacturing geometry constraint that the draft angle is greater than 3 degrees.

[0088] This iterative process is repeated. Assume that after 150 iterations, the objective function converges. The objective function value at this point is... Down to The difference between the previous iteration and the current iteration is less than the set convergence threshold. The optimization process stops, and the final three-dimensional mesh density field is output.

[0089] Based on this converged three-dimensional mesh density field (i.e., the final structural morphology), a confirmatory static analysis is performed. The total reaction force borne by the punch wall under the outer opening is calculated. for Newton, the total reaction force borne by the intermediate ejector for Newton.

[0090] Therefore, the optimal load sharing ratio can be calculated: This indicates that the optimized structure allows the outer perforated punch wall to bear approximately 50.5% of the total load, achieving a nearly equal load distribution between the two components.

[0091] S5. Extract the isosurfaces of a given density threshold in the 3D mesh density field, and constrain the isosurfaces to geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm.

[0092] In a specific embodiment of the present invention, the isosurfaces of a given density threshold in the three-dimensional mesh density field are extracted, and the isosurfaces are constrained into geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes by a geometric feature fitting algorithm, including: traversing the three-dimensional mesh density field and extracting spatial point cloud data with unit density equal to the given density threshold.

[0093] The moving cube algorithm is used to triangulate spatial point cloud data and generate isosurfaces.

[0094] Calculate the local curvature and normal vector characteristics of the isosurface.

[0095] Based on the local curvature and normal vector characteristics, the isosurface is segmented and fitted into geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm.

[0096] Specifically, this step aims to transform the abstract three-dimensional mesh density field generated in step S4, consisting of discrete density values, into a solid model composed of smooth, regular geometric surfaces that meets engineering practice and manufacturing requirements. This process can be seen as a "translation" process from an optimized digital concept to a specific engineering design.

[0097] The process begins by extracting the structure's outline from a 3D mesh density field. This requires setting a given density threshold, which serves as the boundary distinguishing solid material from voids. Typically, this threshold is set to 0.5, meaning that cells with a relative density greater than 0.5 are considered solid material, while those with a density less than 0.5 are considered void areas. By traversing every cell in the 3D mesh density field, all spatial locations where the cell density exactly equals the given density threshold are identified. Since the cell density varies continuously, these locations effectively form a continuous surface. By performing linear interpolation within the cells, the 3D coordinates of these points can be accurately calculated, resulting in a vast spatial point cloud describing the structure's boundaries.

[0098] After obtaining the spatial point cloud data, the next step is to transform it from an unordered set of points into a continuous surface. The moving cube algorithm is used here. This algorithm examines each tiny cube in the 3D mesh density field one by one. Based on the comparison of the density values ​​of these 8 nodes with a given density threshold, there are 256 possible cases. The moving cube algorithm predetermines one or more combinations of small triangular facets for each case. These small triangular facets pass through the interior of the cube, connecting the edges where the density values ​​cross the threshold. By piecing together the small triangular facets generated in all the cubes, a triangular mesh surface that completely covers the structural boundary is constructed, i.e., an isosurface. This isosurface is a preliminary geometric representation of the topology optimization result, but its surface is usually uneven, composed of a large number of small triangular facets, making it inconvenient for direct manufacturing.

[0099] To normalize the model, feature identification and fitting are needed for this coarse isosurface. First, the local curvature and normal vector features of each vertex on the isosurface are calculated. The normal vector describes the orientation of the surface at that point, while the curvature describes the degree of bending of the surface at that point. Principal curvatures are two fundamental quantities describing the bending of the surface, and their calculation can be based on the geometric information of the triangular facets surrounding the vertex. For example, Gaussian curvature... and mean curvature It can be estimated using the following formula: ; ; in, and These are two principal curvatures. By analyzing the principal curvature values ​​at all points within a region, the geometric shape of that region can be determined. For example, if both principal curvatures of a region are close to zero, it is likely a plane; if one principal curvature is constant and the other is zero, it is likely a cylinder; if both principal curvatures are equal and not zero, it is likely a sphere.

[0100] Based on the calculated local curvature and normal vector characteristics, a geometric feature fitting algorithm is used to process the isosurface. This algorithm first uses region-growing as a growth criterion, dividing the entire isosurface into several regions with similar geometric features, where the difference in the angle between the normal vectors and the principal curvature between adjacent triangular facets is less than a set tolerance. For example, it identifies and clusters all triangular facet regions with near-zero curvature and consistent normal vector directions. Then, for each segmented region, the algorithm attempts to fit it with an ideal geometric primitive such as a stepped cylinder, cone, or plane. For example, for a region identified as a plane, the algorithm calculates an optimal fitting plane; for a region identified as a cylinder, it calculates its axis, radius, and other parameters. In this way, the original complex isosurface, composed of countless small triangular facets, is constrained and replaced with a clear model composed of a few standard engineering geometric primitives such as stepped cylinders, cones, or planes. These primitives are defined by precise mathematical equations, satisfying the geometric constraints for subsequent engineering analysis and CNC machining.

[0101] For example, assuming the 3D mesh density field generated in step S4 is ready, it is now geometrically reconstructed. First, a given density threshold is set to 0.5. The entire 3D mesh density field is traversed, and all locations with a density value equal to 0.5 are extracted, forming a spatial point cloud data containing millions of points.

[0102] Next, this spatial point cloud data is input into the moving cube algorithm. After processing, the algorithm generates an isosurface composed of 1.5 million small triangular facets. Although this isosurface shows an optimized structural profile, its surface is rough and cannot be directly used to generate CNC machining code.

[0103] Then, the local curvature and normal vector characteristics of each vertex on the isosurface are calculated. Taking a large region A on the isosurface as an example, calculations show that over 95% of the vertices within this region exhibit the same two principal curvatures. and The values ​​are all in the range of -0.001 to 0.001, and the normal vectors of these vertices are basically in the same direction, all pointing approximately to the (0,0,1) direction.

[0104] Based on these characteristics, the geometric feature fitting algorithm identifies region A as a quasi-plane. Subsequently, the algorithm uses the least squares method to fit the coordinates of all vertices within region A, calculating an optimally fitted plane equation: This means that region A can be replaced by an ideal plane located at a height of z = 50 mm.

[0105] Taking another slender, curved region B as an example, calculations revealed one of its principal curvatures. It is close to 0 at all points, while the other principal curvature This is approximated by a constant of 0.02 (unit: 1 / mm), and the normal vector always points to a common central axis. Based on this, the algorithm identifies region B as a quasi-cylindrical surface. By fitting points within region B, the corresponding cylindrical surface parameters are calculated as follows: the axis direction is the Y-axis, the axis equations are x=30, z=40, and the radius is... Millimeters.

[0106] By repeating this process across all regions of the isosurface, the complex triangular mesh model was successfully transformed into a model composed of several precisely defined geometric primitives representing engineering standards, such as planes, cylinders, and cones. This model possesses clear geometric features and a concise data structure, laying the foundation for subsequent manufacturing and analysis.

[0107] See Figure 4 S6. Identify stress concentration areas in the isosurface, use B-spline curves to perform smooth chamfer transition calculations with continuous curvature, generate smooth reconstructed solid data and output it.

[0108] In a specific embodiment of the present invention, stress concentration regions in the isosurface are identified, and B-spline curves are used to perform smooth chamfer transition calculations with continuous curvature, generating and outputting smooth reconstructed solid data, including: obtaining the stress tensor of each node of the isosurface.

[0109] The stress gradient is obtained by calculating the ratio of the stress difference between adjacent nodes to the distance.

[0110] Compare the stress gradient with a set threshold, and output the set of nodes in the stress concentration region.

[0111] Extract the boundary feature lines of the stress concentration region.

[0112] Obtain discrete data points of the boundary feature lines.

[0113] Construct B-spline basis functions to establish a system of linear equations relating discrete data points and control point coordinates.

[0114] By introducing curvature continuity constraints, the linear equation system is solved to obtain the updated control point coordinates.

[0115] Generate a smooth transition surface based on the updated control point coordinates.

[0116] Perform Boolean operations on the smoothed surface and the aforementioned geometric primitives to generate smoothed reconstructed entity data.

[0117] Transform smoothly reconstructed entity data into a standard entity data exchange format and output it.

[0118] Specifically, this step is the final refinement stage of the entire optimization design process. Its core objective is to eliminate local high-stress points that may be generated during the geometric reconstruction process, ensuring that the final solid model not only has the optimal shape, but also has the expected mechanical properties and durability.

[0119] The process begins with stress checking the isosurface model generated in step S5, composed of the aforementioned geometric primitives, to identify potential stress concentration regions. These regions typically appear at sharp corners where different geometric primitives intersect. To accurately locate these regions, the stress gradient at each node of the isosurface needs to be calculated. First, a final finite element static analysis is performed on the isosurface model to obtain the stress tensor at each node of the model under the peak load condition determined in step S2. The stress tensor is a mathematical object describing the stress state in all directions at a point, from which the equivalent stress, which comprehensively measures the stress level, can be calculated. Subsequently, the stress gradient between adjacent nodes is calculated, which characterizes the rate of stress change in space. The calculation formula is: ; in, It is a node With nodes The stress gradient between them is expressed in Pascals per meter (Pa / m) or megapascals per millimeter (MPa / mm). and Adjacent nodes and The equivalent stress value, in Pascals (Pa) or megapascals (MPa). This is the spatial distance between the two nodes, measured in meters (m) or millimeters (mm). A high stress gradient value means that the stress increases sharply over a very short distance, which is a typical characteristic of stress concentration. By comparing the calculated stress gradient with a pre-set threshold, all areas where the stress gradient exceeds the limit can be identified, i.e., stress concentration areas. This typical value is usually 1000 megapascals / mm² (1000 MPa / mm). This threshold is determined based on the material's fatigue limit and fracture toughness data, ensuring that the identified areas are the weakest points in the structure most likely to develop cracks during long-term service.

[0120] After identifying stress concentration areas, the sharp geometric features of these areas need to be smoothed. Boundary feature lines formed by the intersection of different geometric elements in these areas are extracted; these lines are typically sharp edges.

[0121] Next, B-spline curves are used to replace these boundary feature lines to achieve a smooth chamfer transition with continuous curvature. A B-spline curve is a smooth parametric curve determined by a series of control points, possessing excellent local controllability and smoothness. To generate this transition surface, discrete data points are first collected from the boundary feature lines. Then, B-spline basis functions are constructed, and a system of linear equations connecting these discrete data points with the coordinates of the unknown control points is established. The mathematical expression for the B-spline surface is: ; in, It is a generated smooth transition surface; and These are the parametric coordinate variables of the surface; These are the coordinates of the control point grid to be solved; and Is and The order of the B-spline basis functions in the direction and Usually, a value of 3 or 4 is used to ensure curvature continuity; and This refers to the number of control points. To solve for the control points... In addition to requiring the surface to pass through or approximate discrete data points, a curvature continuity constraint must also be introduced. Mathematically, this constraint requires the transition surface and the original geometric primitives to satisfy the following conditions at their interface: Geometric continuity, that is, the coincidence of the boundary nodes at the intersection. Continuous, the tangent planes determined by the first derivative are coplanar ( (continuous), and the magnitude and direction of the principal curvature determined by the second derivative are exactly the same ( (Continuous). This constraint requires that at the connection points between the transition surface and the original geometric primitives, not only should the tangent directions be consistent, but the curvature changes should also be smoothly transitioned, thus mathematically guaranteeing that there are no geometrical abrupt changes at the connection points. Solving this constrained system of linear equations yields the updated coordinates of the control points, thereby generating an ideal, smooth transition surface.

[0122] After generating a smooth transition surface to replace the sharp corners, it is fused with the engineering-standard geometric primitive model generated in step S5 using Boolean operations in 3D modeling. The Boolean operations precisely remove the sharp corners in the original model and geometrically stitch the newly generated smooth surface with the original model, ultimately generating a smooth reconstructed solid data that maintains the optimized topology while eliminating local stress concentrations.

[0123] Finally, this final, smoothly reconstructed entity data is converted into a standard entity data exchange format commonly used in the industry, such as STEP (Standard for the Exchange of Product Model Data) or IGES (Initial Graphics Exchange Specification), and then exported. This output file can be directly read by various computer-aided manufacturing (CAM) software to generate machine tool machining paths or for further detailed engineering analysis.

[0124] For example, following the example from step S5, the generated model contains a plane. It intersects with a cylindrical surface with a radius of 50 mm, forming a sharp 90-degree internal turn.

[0125] First, stress analysis is performed on the model to identify stress concentration areas. Near the corner, two adjacent nodes, M and N, are selected. Node M is located on a plane, 0.1 mm from the corner line, and its equivalent stress... The equivalent stress is 150 megapascals (MPa). Node N is located on the cylindrical surface, also 0.1 mm from the corner line, and its equivalent stress is... The distance between two points is 480 megapascals. It is approximately 0.2 millimeters.

[0126] Calculate the stress gradient between two points: Assuming a stress gradient threshold of 1000 MPa / mm based on the material properties of H13 mold steel, the corner region is marked as a stress concentration area because the calculated value of 1650 is greater than 1000.

[0127] Next, the boundary feature line of this region is extracted, namely the circular edge where the plane and the cylindrical surface intersect. Then, a smooth chamfer transition calculation is performed on this circular edge using a B-spline curve. By collecting a series of discrete data points on this edge and setting constraints to achieve curvature continuity at the connection between the transition surface and the original plane and cylindrical surfaces, the coordinates of the control points of the B-spline surface are solved. The calculation result generates a chamfered surface with continuously changing curvature, replacing the original 90-degree sharp corner.

[0128] Next, Boolean operations are performed. This newly generated chamfered surface is merged with the original planar and cylindrical surface models. The operation automatically removes sharp corners and seamlessly integrates the chamfered surface, generating smooth reconstructed solid data. Finally, this final model containing the smooth chamfer is output as a standard solid data exchange format file named "Optimized_Retainer_Final.stp". This step is now complete.

[0129] The above content is merely an example and illustration of the concept of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described, or use similar methods to replace them, as long as they do not deviate from the concept of the invention or exceed the scope defined by the present invention, and all such modifications and additions should fall within the protection scope of the present invention.

Claims

1. A method for topology optimization of stress distribution in a flange-type sliding sleeve ejection structure, characterized in that, include: S1. Construct an initial three-dimensional finite element model of the flange-type sliding sleeve and the mold, obtain the physical parameters in the warm forging process and input them into the initial three-dimensional finite element model, perform transient thermo-mechanical coupling analysis, and extract the spatial dynamic frictional resistance field and thermal shrinkage stress field. S2. Map the spatial dynamic frictional resistance field and thermal shrinkage stress field to the time-varying boundary load conditions of the topology optimization design domain of the material ejection structure, and generate a dynamic load mapping model. S3. Set the outer opening punch and the middle ejector as a dual-solid design domain for collaborative topology optimization. Define the contact surface as the sliding friction boundary. Based on the variable density method and dynamic load mapping model, establish a topology optimization model with the objective function of minimizing the maximum equivalent stress square difference between the two entities and the constraint condition of volume retention rate. S4. During the iterative optimization process of the topology optimization model, the material density field is adaptively redistributed between the intermediate ejector and the outer hole punch wall to calculate the optimal load sharing ratio and embed the processing and manufacturing geometric constraints to generate a three-dimensional mesh density field. S5. Extract the isosurfaces from the three-dimensional mesh density field with a given density threshold, and constrain the isosurfaces to geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm. S6. Identify stress concentration regions in the isosurface, use B-spline curves to perform smooth chamfer transition calculations with continuous curvature, generate smooth reconstructed solid data and output it.

2. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: The initial three-dimensional finite element model of the flange-type sliding sleeve and the mold is constructed. The physical parameters in the warm forging process are obtained and input into the initial three-dimensional finite element model. Transient thermo-mechanical coupling analysis is performed to extract the spatial dynamic frictional resistance field and thermal shrinkage stress field, including: Obtain the geometric dimensions of the flange-type sliding sleeve and the structural parameters of the mold, and construct an initial three-dimensional finite element model; The heating temperature, upsetting and flange pressing tonnage, and temperature control cooling rate in the warm forging process are collected as physical parameters. Input the physical parameters into the initial three-dimensional finite element model, run the transient thermo-mechanical coupling analysis algorithm, and simulate the instantaneous demolding of the sliding sleeve and the continuous ejection stage. Extract the spatial dynamic frictional resistance field and thermal contraction stress field from the simulation results.

3. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: The process of mapping the spatial dynamic frictional resistance field and thermal contraction stress field to the time-varying boundary load conditions of the topology optimization design domain of the material removal structure, and generating a dynamic load mapping model, includes: Extract the boundary node coordinates of the topology optimization design domain of the material return structure; Based on the boundary node coordinates, the spatial dynamic frictional resistance field and the thermal contraction stress field are spatially interpolated to obtain the node load vector; Serialize the nodal load vectors according to the time step to construct time-varying boundary load conditions; Time-varying boundary load conditions are applied to the topology optimization design domain of the material removal structure to generate a dynamic load mapping model.

4. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: The outer-perforated punch and the intermediate ejector are set as a two-solid design domain for collaborative topology optimization. The contact surface is defined as a sliding friction boundary. Based on the variable density method and dynamic load mapping model, a topology optimization model is established with the minimization of the maximum equivalent stress square difference between the two solids as the objective function and the volume retention rate as the constraint condition. This includes: The grid cells of the outer opening punch and the middle ejector are divided and set as a dual-solid design domain; Extract the interface between the two solid design domains and define it as the sliding friction boundary; Based on the variable density method, the relative density of the mesh elements is used as the design variable, and the equivalent stress of each element is calculated by combining the dynamic load mapping model. Calculate the maximum equivalent stress between the outer opening punch and the middle ejector, and construct an objective function to minimize the square difference of the maximum equivalent stress between the two entities; Set a threshold for the ratio of total volume to initial volume, construct volume retention rate constraints, and generate a topology optimization model.

5. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: In the iterative optimization process of the topology optimization model, the material density field is adaptively redistributed between the intermediate ejector and the outer opening punch wall, the optimal load sharing ratio is calculated, and manufacturing geometric constraints are embedded to generate a three-dimensional mesh density field, including: Calculate the sensitivity of the objective function of the topology optimization model with respect to the design variables; The material density field is updated based on sensitivity to achieve adaptive redistribution; Based on the updated material density field, the bearing capacity of the intermediate ejector and the outer opening punch wall is calculated, and the optimal load sharing ratio is determined. During the density update process, the draft angle along the axial direction is forced to be greater than zero and the generation of internal holes is restricted as a geometric constraint for processing and manufacturing. After multiple iterations until the objective function converges, a three-dimensional mesh density field is output.

6. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: The extraction of isosurfaces from a given density threshold in a three-dimensional mesh density field, followed by constraining these isosurfaces into geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm, includes: Traverse the 3D grid density field and extract spatial point cloud data with a unit density equal to a given density threshold. The moving cube algorithm is used to triangulate spatial point cloud data and generate isosurfaces; Calculate the local curvature and normal vector characteristics of the isosurface; Based on the local curvature and normal vector characteristics, the isosurface is segmented and fitted into geometric primitives selected from stepped cylindrical surfaces, conical surfaces, or planes using a geometric feature fitting algorithm.

7. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 1, characterized in that: The process involves identifying stress concentration regions in the isosurface, performing smooth chamfer transition calculations using B-spline curves to generate and output smooth reconstructed solid data, including: Calculate the stress gradient at each node of the isosurface and mark the region where the stress gradient exceeds a set threshold as a stress concentration region; Extract the boundary feature lines of the stress concentration region; B-spline curves are used to fit the boundary feature lines, and the coordinates of control points are calculated to achieve a smooth chamfer transition with continuous curvature. Perform Boolean operations on the smoothed surface and the aforementioned geometric primitives to generate smoothed reconstructed entity data; Transform smoothly reconstructed entity data into a standard entity data exchange format and output it.

8. The method for topology optimization of stress distribution in a flange-type sliding sleeve ejection structure according to claim 7, characterized in that: The calculation of the stress gradient at each node of the isosurface, and marking regions where the stress gradient exceeds a set threshold as stress concentration regions, includes: Obtain the stress tensor at each node of the isosurface; The stress gradient is obtained by calculating the ratio of the stress difference between adjacent nodes to their distance. Compare the stress gradient with a set threshold, and output the set of nodes in the stress concentration region.

9. The method for topology optimization of force distribution in a flange-type sliding sleeve ejection structure according to claim 7, characterized in that: The method of fitting boundary feature lines with B-spline curves and calculating control point coordinates to achieve a smooth chamfer transition with continuous curvature includes: Obtain discrete data points of the boundary feature lines; Construct B-spline basis functions and establish a system of linear equations relating discrete data points and control point coordinates; By introducing curvature continuity constraints, the linear equation system is solved to obtain the updated control point coordinates; Generate a smooth transition surface based on the updated control point coordinates.

Images