Precision control system and method for large-curvature and shallow-drawing formed part

Abstract

The invention belongs to the technical field of workpiece precision control, and discloses a precision control system and method for a large-curvature and shallow-drawing formed part, and the method comprises the steps: obtaining initial geometric parameters and material mechanical property parameters of the formed part; establishing a finite element analysis model of the formed part; discretizing the finite element analysis model of the formed part into n nodes, and calculating stress distribution at the n nodes of the formed part in the unloading process, strain energy generated by rebound and rebound bending moment formed by the strain energy; based on the stress characteristics of the n nodes, carrying out characteristic clustering on the n nodes to obtain m node areas; solving the overall rebound deformation and rebound strain energy of the m node areas by using an eagle-eye algorithm; the tolerance range is preset, the rebound bending rigidity of the corresponding area of the node area with the rebound deformation exceeding the tolerance range is increased on the basis of the characteristic geometrical shape design, the deformation generated by rebound bending moment is offset, and the size precision and shape quality of the large-curvature and shallow-drawing forming part are effectively improved.

Description

Technical Field

[0001] The present invention relates to the technical field of workpiece precision control, and more specifically, to a precision control system and method for large curvature and shallow drawn formed parts. Background Art

[0002] The patent application publication number CN109513931A discloses a control method for residual thermal stress and induced deformation in additive manufacturing, including the following steps: using experimental or numerical simulation methods to obtain residual stress distribution and stress value of formed parts; designing a continuous or discontinuous flexible substrate structure according to the specific structural form of the formed parts, so that the equivalent strain of the designed flexible structure is lower than the allowable deformation range, and has an equivalent elastic modulus lower than that of the main material, and by designing multiple groups of flexible substrate structures, the residual stress partition release and active control based on the flexible structure are realized. In the additive forming process, the characteristic size of the flexible structure is gradually changed, the structural flexibility coefficient is reduced, and a gradient flexible structure is obtained to meet the service load requirements of the component. Through structural deformation, spatial discretization and release of residual stress, within the allowable deformation range, discrete control of spatial macro thermal stress is realized.

[0003] During the metal sheet forming process, the formed part will produce reverse bending deformation when unloading. This deformation is closely related to the strain energy and residual stress inside the formed part. However, existing control methods often ignore this point and cannot accurately predict and evaluate the trend of reverse bending deformation, resulting in difficulty in controlling the dimensional accuracy and shape quality of the formed part. Different areas of the formed part may produce different degrees of strain energy and residual stress during the forming process. Existing methods often treat the entire formed part as a whole and cannot optimize different areas in a targeted manner, thus affecting the overall accuracy and quality of the formed part. Taking automobile body panels as an example, due to their large curvature radius and drawing ratio, they are very likely to produce large reverse bending deformation after stamping, resulting in interference with other components during assembly, excessive gaps and other problems. If the reverse bending deformation cannot be effectively controlled, it will seriously affect the overall quality and performance of the car.

[0004] In view of this, the present invention proposes a precision control system and method for large curvature and shallow drawn parts to solve the above problems. Summary of the invention

[0005] In order to overcome the above-mentioned defects of the prior art and to achieve the above-mentioned purpose, the present invention provides the following technical solutions: A precision control system for large curvature and shallow drawing formed parts, comprising: a model building module for obtaining initial geometric parameters and material mechanical property parameters of the formed parts; and establishing a finite element analysis model of the formed parts;

[0006] Energy analysis model, used to discretize the finite element analysis model of the formed part into n nodes, and calculate the strain energy and residual stress at the n nodes of the formed part during the unloading process;

[0007] The deformation fitting module clusters n nodes based on the strain energy and residual stress at n nodes to obtain m node areas; the Hawkeye algorithm is used to solve the reverse bending deformation of the m node areas;

[0008] The detection and adjustment module is used to preset the tolerance range and add a stiffness reinforcement structure to the node area where the reverse bending deformation exceeds the tolerance range. The modules are connected by wire and / or wireless means.

[0009] Furthermore, the initial geometric parameters include length, width, thickness and radius of curvature; the material mechanical property parameters include elastic modulus, Poisson's ratio, yield strength and plasticity parameters.

[0010] Furthermore, the method of establishing the finite element analysis model of the formed part includes:

[0011] According to the initial geometric parameters of the formed part, a three-dimensional solid model of the formed part is established; based on the mechanical property parameters of the material, material properties are defined on the three-dimensional solid model to obtain a preliminary finite element model;

[0012] For the preliminary finite element model, its outer surface is extracted as the outer boundary surface, and the inner boundary surface is extracted from the preliminary finite element model by using geometric Boolean operation; based on the outer boundary surface and the inner boundary surface, the intermediate surface is extracted;

[0013] Predefine the parameter domain, project the extracted intermediate surface into the parameter domain in an equiangular manner, and perform meshing in the parameter domain to generate a mesh topological relationship of rectangular units;

[0014] The method of performing grid division includes:

[0015] Nodes are inserted evenly inside and on the boundary of the parameter domain and numbered. Every four adjacent nodes form a rectangular unit. For each rectangular unit, the numbers of the nodes at its four vertices are determined, and the topological relationship of the unit is determined according to the node numbers.

[0016] Define the physical domain, and construct shell elements in the physical domain according to the mesh topology of the rectangular elements. For each node of the rectangular element in the parameter domain, calculate the corresponding coordinates of the physical domain through the surface equation, and the coordinates of the physical domain are the control point coordinates of the shell element. Use the node numbers in the parameter domain as the node numbers of the shell element in the physical domain.

[0017] For each shell element i, construct its shape function NF(i);

[0018] NF(i)=∑I w_I(α)×N_I(α)+∑ J w_J(β)×N_J(β); where I is the number of the node in the parameter direction α, N_I(α) is the one-dimensional B-spline basis function in the parameter direction α, and w_I(α) is the weight of the node numbered I in the parameter direction α; J is the number of the node in the parameter direction β, N_J(β) is the one-dimensional B-spline basis function in the parameter direction β, and w_J(β) is the weight of the node numbered J in the parameter direction β;

[0019] Obtain the strain-displacement matrix B and stress-strain matrix D of the shell element; calculate the stiffness matrix K_i of the shell element i based on the strain-displacement matrix B and stress-strain matrix D of the shell element, and calculate the load matrix F_i based on the shape function of the shell element;

[0020] Stiffness matrix K_i=∫(B T ·D·B)dδ; where B T is the transpose of the strain-displacement matrix B; dδ represents the differential integral on the parameter domain δ; ∫()dδ represents the parameter domain integration on the parameter domain δ;

[0021] Load matrix F_i=∫((NF(i)) T b)dδ; where b is physical strength, (NF(i)) T is the transpose of NF(i);

[0022] The stiffness matrices and load matrices of all shell elements are assembled separately to obtain the global stiffness matrix K and the global load matrix F, and the continuity between shell elements is processed to obtain the finite element analysis model of the formed part.

[0023] Furthermore, the method of extracting the intermediate surface includes:

[0024] Define a signed distance factor d(x,y,z) to represent the shortest distance from any point L in the preliminary finite element model to the inner boundary surface and the outer boundary surface; where x, y and z are the horizontal coordinate, vertical coordinate and vertical coordinate of any point L in the spatial coordinate system of the preliminary finite element model, respectively;

[0025] Define midsurface deviation Wherein, d_in is the shortest distance from any point L to the inner boundary surface, d_out is the shortest distance from any point L to the outer boundary surface, w_in is the weighting coefficient of the inner boundary surface, and w_out is the weighting coefficient of the outer boundary surface;

[0026] Discretize the signed distance factor sampling into a three-dimensional voxel grid as the distance field, the three-dimensional voxel grid contains N voxels, each voxel corresponds to a scalar value, that is, the distance field value at the center point of the voxel;

[0027] For each voxel, check whether there is an intersection point with d(x, y, z) = t between it and the adjacent voxel. If there is an intersection point, use linear interpolation to calculate the coordinates of the intersection point; smoothly connect the calculated coordinates of the intersection point to obtain the intermediate surface.

[0028] Furthermore, the method of discretizing the finite element analysis model of the formed part into n nodes includes:

[0029] Performing initial mesh division on the finite element analysis model of the formed part to obtain an initial mesh;

[0030] Define the finite element function, the formula of the finite element function is: K·uh=f; by solving the finite element function on the initial grid, the initial solution uh is obtained; where f is the load vector;

[0031] Based on the initial solution uh, for each shell element i, the corresponding error estimate is calculated

[0032] Among them, u is the true solution, C is the elastic matrix, ε(u-uh) is the strain tensor of u-uh;

[0033] Preset relative error kl, for The shell element i is marked as needing to be encrypted, and the shell element i marked as needing to be encrypted is encrypted and refined to obtain a new mesh;

[0034] The encryption refinement method includes:

[0035] The entire computational domain is represented by a hexahedral unit as the root node, and the root node is equally divided and refined to generate 8 child hexahedral units as the child nodes of the root node. Each child node is recursively equally divided and refined until the preset maximum refinement level is reached; that is, a complete hexahedral tree hierarchy is constructed, and each node represents a hexahedral unit;

[0036] The shell element i marked as needing to be encrypted is mapped to a node in the hexahedral tree hierarchy, that is, the node is marked. For each marked node, 8 new child nodes are generated, representing the 8 refined child hexahedral elements. All the child nodes of the original node are deleted and replaced with the generated 8 new child nodes, and the process is repeated until the preset maximum refinement level is reached.

[0037] Preset error threshold The finite element function is re-solved on the new grid to obtain a new solution uh′. If the corresponding error estimate calculated for each shell element i based on the new solution uh′ is greater than or equal to the error threshold Repeat the encryption and refinement until the calculated error estimate is less than the error threshold Stop the encryption and refinement to obtain the final mesh, which contains n nodes.

[0038] Furthermore, the calculation method of the strain energy and residual stress at the n nodes includes:

[0039] The finite element function is solved on the final mesh to obtain the final solution uh″; the final solution uh″ contains the node displacement vectors of all shell elements; based on the final solution uh″, the strain field {∈}=uh″·B and the stress field {ρ}=uh″·D of each shell element are calculated;

[0040] Calculate the strain energy U of each shell element based on the strain field {∈} and stress field {ρ} of each shell element;

[0041] Where ∫()dV represents the volume integral of the volume V of the shell element; {∈} T is the transpose of the strain field {∈}, [DL(μp,R,θ)] is the nonlinear constitutive matrix, μp is the plastic strain, R is the temperature, and θ is the strain rate; JY is the geometric mapping Jacobian of the shell element;

[0042] [DL(μp,R,θ)]=[De({∈},R)]+[Dp(μp,R,θ)]+[Dvp({∈},μp,R,θ)]; where [De({∈},R)] is the elastic matrix; [Dp(μp,R,θ)] is the plastic matrix; [Dvp({∈},μp,R,θ)] is the viscoplastic matrix;

[0043] Based on the plastic matrix and the elastic matrix, the residual stress ER of each shell element is calculated;

[0044] Where {τp} is the plastic strain vector, is the tensor product operator;

[0045] For each shell element, the average strain energy and average residual stress of each shell element are calculated based on the calculated strain energy and residual stress. Based on the average strain energy and average residual stress of each shell element, the interpolation function is used at n nodes to perform interpolation operations to obtain the strain energy and residual stress of each node.

[0046] Furthermore, the method of clustering the n nodes includes:

[0047] The strain energy and residual stress of each node constitute a w-dimensional feature vector respectively; the feature vectors of n nodes are defined as the initial positions of n fireflies in the w-dimensional space;

[0048] Define the number of clusters K and the clustering objective function fs;

[0049] Among them, C_j is the center of the jth cluster, X_p1 is the position of the p1th firefly, and PUX is the inter-cluster penalty function;

[0050] Where q is the index of the cluster, w_jq is the weight coefficient between the jth cluster and the qth cluster; Dd_jq is the Euclidean distance between the jth cluster and the qth cluster, β1 is the distance order; EX is the energy difference penalty function; is the penalty coefficient, N(j) is the number of nodes in the jth cluster, and N(q) is the number of nodes in the qth cluster;

[0051] Energy difference penalty function EX = (1 + μ2 × ΔE_jq + (1-μ2) × ΔQ_jq); where ΔE_jq is the difference between the average strain energy of the jth cluster and the qth cluster, ΔQ_jq is the difference between the average residual stress of the jth cluster and the qth cluster; μ2 is the weight adjustment parameter;

[0052] Define that the brightness of a firefly is inversely proportional to the value of the objective function; calculate the attraction γ_(p1,p2) between the p1th firefly and the p2th firefly;

[0053] γ_(p1,p2)=a×exp(-c′×(RQ(p1,p2)) 2 )×(1+τ2×cos(V(p1,p2)-ε2)); where RQ(p1,p2) is the Euclidean distance between the p1th firefly and the p2th firefly, a is the preset maximum attraction, c′ is the preset light absorption coefficient; τ2 is the adjustment parameter, V(p1,p2) is the direction angle between the p1th firefly and the p2th firefly, and ε2 is the preset expected direction angle;

[0054] Move each firefly according to the attractiveness of all other fireflies. The formula for moving is: Among them, α3 is the step size factor, rd is a random number in the interval [0, 1], X′_p1 is the new position of the p1th firefly, and X_p2 is the position of the p2th firefly before moving;

[0055] Calculate the value of the clustering objective function for the new position of the firefly. If the new position is higher than the value of the clustering objective function at the original position, update the firefly to the new position; and update the brightness of each firefly according to the new position of the firefly;

[0056] Iterate until the preset number of iterations is met to obtain the final position of the firefly. According to the final position of the firefly, the adjacent fireflies are grouped into the same cluster. Each cluster corresponds to a node area, thus obtaining m node areas.

[0057] Furthermore, the method for calculating the reverse bending deformation includes:

[0058] For each node area, the coordinates of all nodes in the node area are extracted as eagles; the parameters of the eagle-eye algorithm are initialized, including the scale of the eagle group, that is, the number of eagles M4, the maximum number of iterations and the discovery rate;

[0059] Encode each eagle using real number coding or binary coding, and the length of the code is equal to the number of nodes in the node area multiplied by 3;

[0060] Define the competition function BN; calculate the value of the competition function of each eagle, recorded as the competition value; sort the eagles in the eagle group according to the competition value, and update the position of each eagle according to the position of the eagle with the highest competition value in the current eagle group and its own current position;

[0061] After each position update, the competition value of each eagle in the eagle group is calculated again. According to the preset territory occupancy rate ra, the first ra×M4 eagles are retained, and the rest are eliminated;

[0062] For each eagle retained in each iteration, mutation is performed with a preset probability rp to obtain the new position of the eagle;

[0063] At this time, the eagles form a new generation of eagles. The competition value of the eagles in the new generation of population is calculated, and the eagle with the highest competition value is retained as the seed eagle of the next generation of eagles.

[0064] Repeat the iteration until the number of iterations reaches the defined maximum number of iterations to obtain the final eagle group, calculate the competition value of each eagle in the final eagle group, take the code corresponding to the eagle with the highest competition value as the optimal value, decode the optimal value to obtain the optimal plane fitting coordinates of the nodes in the node area; take the difference between the coordinates of each node in the node area and the optimal plane fitting coordinates as the node's inverse bending deformation; take the average of the inverse bending deformation of each node in the node area as the inverse bending deformation of the node area.

[0065] Furthermore, the formula of the competition function BN is:

[0066] Among them, w′, w0

[0067] and w1 are the parameters for finding equilibrium weights; En(U_g,ρ′_g) is the deformation state function, U_g is the strain energy of node g in the node region, ρ′_g is the residual stress of node g in the node region, is the curvature of node g in the node region, Cu(g,P′_g) is the plane deviation function, and P′_g is the fitting plane for node g in the node region;

[0068] Deformation state function En(U_g,ρ′_g)=z1×∑ g U_g+(1-z1)×∑ g ρ′_g; where z1 is the deformation equilibrium parameter;

[0069] Cu(g,P_g)=x1×(∑ g ||gP′_g||)+(1-x1)×max(||gP′_g||); where,

[0070] x1 is the plane balance parameter;

[0071] The formula for updating the position is:

[0072] X(it+1)=X(it)+AK×cos(2π×r_2)+(2×r_3)×(X_r-X(it))+DPW; where X(it+1) is the position of the eagle X at the it+1th iteration, X(it) is the position of the eagle X at the itth iteration, π is pi, and AK is the sweep radius factor;

[0073] The scanning radius factor AK = 2a1 × r_1-a1; where a1 is a constant, DPW is a global factor; X_r is the position of the eagle with the highest competition value in the current eagle group;

[0074] Global factor DPW = (2×r_4)×(X_l-X(it))+(2×r_5)×(X_g-X(it)); where r_1, r_2, r_3, r_4 and r_5 are all random update numbers in the interval [0, 1];

[0075] The formula for performing mutation is:

[0076] X_new=X_old+SF×(X_old-2×r_6×(ub-lb)); where SF is the scaling factor, r_6 is the mutation random number in the interval [0, 1], ub is the upper boundary of the node's coordinates, and lb is the lower boundary of the node's coordinates; X_new is the position of the eagle after the mutation, and X_old is the position of the eagle before the mutation.

[0077] A method for controlling the precision of a large curvature and shallow drawing formed part is implemented based on the large curvature and shallow drawing formed part precision control system, comprising: S1, obtaining initial geometric parameters and material mechanical property parameters of the formed part; and establishing a finite element analysis model of the formed part;

[0078] S2, discretizing the finite element analysis model of the formed part into n nodes, and calculating the strain energy and residual stress at the n nodes of the formed part during the unloading process;

[0079] S3. Based on the strain energy and residual stress at the n nodes, the n nodes are clustered to obtain m node regions; the reverse bending deformation of the m node regions is solved by the Hawkeye algorithm;

[0080] S4. Preset the tolerance range and increase the stiffness of the node area where the reverse bending deformation exceeds the tolerance range to strengthen the structure.

[0081] The technical effects and advantages of the precision control system and method of large curvature and shallow drawing formed parts of the present invention are as follows:

[0082] The present invention can improve the forming accuracy of large curvature and shallow drawn parts, and effectively reduce the occurrence of deformation defects; through modeling analysis, the finite element analysis model of the formed part is discretized into multiple nodes, and the strain energy and residual stress at each node are analyzed, which can more accurately predict and evaluate the reverse bending deformation trend of the formed part, and realize the zoning processing of deformation in different areas, thereby improving the flexibility and pertinence of deformation control; secondly, the efficient solution algorithm of the reverse bending deformation improves the calculation efficiency and solution accuracy, and can quickly and accurately solve complex reverse bending deformation, greatly improving the analysis efficiency and shortening the product design cycle; by increasing the stiffness enhancement structure, the stiffness of the specific area is adjusted in a targeted manner, and the deformation control effect is further optimized. In short, it can not only significantly improve product quality and production efficiency, but also has many advantages such as strong flexibility and excellent computing performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 It is a schematic diagram of a precision control system for large curvature and shallow drawn formed parts of the present invention;

[0084] Figure 2 It is a schematic diagram of a method for controlling the precision of large curvature and shallow drawn parts according to the present invention. DETAILED DESCRIPTION

[0085] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0086] Example 1

[0087] See also Figure 1 As shown, the precision control system for large curvature and shallow drawn parts described in this embodiment includes:

[0088] Model building module, used to obtain the initial geometric parameters and material mechanical performance parameters of the formed part; and to establish a finite element analysis model of the formed part;

[0089] Energy analysis model, used to discretize the finite element analysis model of the formed part into n nodes, and calculate the strain energy and residual stress at the n nodes of the formed part during the unloading process;

[0090] The deformation fitting module clusters n nodes based on the strain energy and residual stress at n nodes to obtain m node areas; the Hawkeye algorithm is used to solve the reverse bending deformation of the m node areas;

[0091] The detection and adjustment module is used to preset the tolerance range, add a stiffness reinforcement structure to the node area where the reverse bending deformation exceeds the tolerance range, and improve the stiffness of the corresponding node area; each module is connected by wired and / or wireless means to realize data transmission between modules.

[0092] Furthermore, the initial geometric parameters include length, width, thickness and radius of curvature; and the material mechanical property parameters include elastic modulus, Poisson's ratio, yield strength and plasticity parameters.

[0093] The initial geometric parameters are directly measured and obtained; the material mechanical properties parameters are given in the relevant material manual and can be directly obtained by reference; if the required parameters cannot be directly found, they can be obtained indirectly through inverse analysis; for example, if the stress-strain response of a material under certain specific working conditions is known, the material parameter combination that satisfies this response can be inferred.

[0094] Furthermore, the method of establishing the finite element analysis model of the formed part includes:

[0095] According to the initial geometric parameters of the formed part, a three-dimensional solid model of the formed part is established; based on the material mechanical property parameters, material properties are defined on the three-dimensional solid model to obtain a preliminary finite element model (with the properties of the three-dimensional solid model).

[0096] An intermediate surface is extracted from the preliminary finite element model, and a triangular mesh surface is generated; specifically, for the preliminary finite element model, its outer surface is extracted as an outer boundary surface, and an inner boundary surface is extracted from the preliminary finite element model using geometric Boolean operations.

[0097] Define the signed distance factor d(x,y,z) to represent the shortest distance from any point L in the preliminary finite element model to the inner boundary surface and the outer boundary surface; for points located inside the preliminary finite element model, the signed distance factor d(x,y,z) takes a positive value; for points located outside, the signed distance factor d(x,y,z) takes a negative value, and on the boundary of the preliminary finite element model, the signed distance factor d(x,y,z)=0; where x, y and z are the horizontal, vertical and vertical coordinates of the arbitrary point L in the spatial coordinate system of the preliminary finite element model, respectively.

[0098] Define midsurface deviation Among them, d_in is the shortest distance value from any point L to the inner boundary surface, d_out is the shortest distance value from any point L to the outer boundary surface, w_in is the weighting coefficient of the inner boundary surface, and w_out is the weighting coefficient of the outer boundary surface; the deviation of the intermediate surface is used as the objective function, and the least squares method is used to solve the value of the weighting coefficient that can minimize the objective function to obtain the optimal weighting coefficient.

[0099] The signed distance factor sampling is discretized into a three-dimensional voxel grid as a distance field. The three-dimensional voxel grid contains N voxels, and each voxel corresponds to a scalar value, that is, the distance field value at the center point of the voxel.

[0100] For each voxel, check whether there is an intersection point with d(x, y, z) = t between it and the adjacent voxel. If there is an intersection point, use linear interpolation to calculate the coordinates of the intersection point; smoothly connect the calculated coordinates of the intersection point to obtain the intermediate surface.

[0101] A parameter domain is predefined (usually a rectangular domain or a triangular domain is selected), the extracted intermediate surface is conformally projected into the parameter domain, and meshing is performed in the parameter domain to generate a mesh topological relationship of rectangular units; the conformal projection is intended to keep the angle of the surface in the parameter domain unchanged.

[0102] The methods of meshing include:

[0103] Nodes are uniformly inserted inside and on the boundary of the parameter domain and are numbered. Every four adjacent nodes form a rectangular unit. For each rectangular unit, the numbers of the nodes at its four vertices are determined, and the topological relationship of the unit is determined according to the node numbers.

[0104] That is, (m+1)×(n+1) nodes and m×n rectangular units are generated in the parameter domain; each unit is determined by the numbers of four nodes, thereby constructing the topological relationship of the unit.

[0105] Define the physical domain (the actual three-dimensional space domain, which corresponds one-to-one to the parameter domain), and construct shell elements in the physical domain according to the mesh topological relationship of the rectangular elements. Specifically, for each node of the rectangular element in the parameter domain, calculate the coordinates of the corresponding physical domain through the surface equation, and these coordinates of the physical domain are the control point coordinates of the shell element. The node numbers in the parameter domain are used as the node numbers of the shell elements in the physical domain.

[0106] For each shell element i, construct its shape function NF(i);

[0107] NF(i)=∑ I w_I(α)×N_I(α)+∑ J w_J(β)×N_J(β); where I is the number of the node in the parameter direction α, N_I(α) is the one-dimensional B-spline basis function in the parameter direction α, and w_I(α) is the weight of the node numbered I in the parameter direction α. ​​During the calculation process, the value of the weight is dynamically adjusted so that the shape function can more accurately express the deformed geometry and improve the calculation accuracy. J is the number of the node in the parameter direction β, N_J(β) is the one-dimensional B-spline basis function in the parameter direction β, and w_J(β) is the weight of the node numbered J in the parameter direction β. It should be noted that the parameter directions α and β are two coordinate directions in the parameter domain, which are usually taken as orthogonal. By taking different values ​​of α and β in the definition domain, a series of parameter points can be generated in the parameter domain. The shape function NF(i) has stronger geometric expression ability while retaining piecewise continuity and affine invariance.

[0108] The strain-displacement matrix B and stress-strain matrix D of the shell element are obtained. The strain-displacement matrix and stress-strain matrix are derived through the shape function of the shell element and the material constitutive relationship. The stiffness matrix K_i of the shell element i is calculated based on the strain-displacement matrix B and stress-strain matrix D of the shell element, and the load matrix F_i is calculated based on the shape function of the shell element.

[0109] Stiffness matrix K_i=∫(B T ·D·B)dδ; where B T is the transpose of the strain-displacement matrix B; dδ represents the differential integral over the parameter domain δ; ∫()dδ represents the parameter domain integration over the parameter domain δ.

[0110] Load matrix F_i=∫((NF(i)) Tb)dδ; where b is the body force (the volume force acting on the volume of the shell element, which is a vector field, including gravity, inertia and electromagnetic force), (NF(i)) T is the transpose of NF(i).

[0111] The stiffness matrices and load matrices of all shell elements are assembled separately to obtain the global stiffness matrix K and the global load matrix F, and the continuity between shell elements is processed (two adjacent shell elements share a boundary to ensure that the displacement on this boundary is continuous) to obtain the finite element analysis model of the formed part.

[0112] Furthermore, the method of discretizing the finite element analysis model of the formed part into n nodes includes:

[0113] The finite element analysis model of the formed part is initially meshed to obtain an initial mesh; the initial meshing method can adopt a mapping method, a feedforward method or a feedforward mapping method; a group of nodes generated during the initial meshing are distributed relatively evenly.

[0114] Define the finite element function, the formula of the finite element function is: K·uh=f; by solving the finite element function on the initial grid, the initial solution uh (node ​​displacement field) is obtained; where f is the load vector (external load acting on the structure).

[0115] Based on the initial solution uh, for each shell element i, the corresponding error estimate is calculated

[0116] Among them, u is the true solution, C is the elastic matrix, and ε(u-uh) is the strain tensor of u-uh. Since the true solution u is unknown, an approximate quantity is constructed and obtained using recovery technology or projection technology.

[0117] Preset relative error kl (based on accuracy requirements), for The shell element i marked as needing to be encrypted is encrypted and refined to obtain a new mesh.

[0118] The encryption refinement method includes:

[0119] The entire calculation domain (the boundary of the finite element analysis model of the formed part describes the geometric shape and scope of the entire calculation domain) is represented by a hexahedral unit as the root node. The root node is equally divided and refined to generate 8 sub-hexahedral units as the child nodes of the root node. Each sub-node is recursively equally divided and refined until the preset maximum refinement level is reached; that is, a complete hexahedral tree hierarchy is constructed, and each node (root node and child node) represents a hexahedral unit.

[0120] The shell element i marked as needing encryption is corresponded to the node in the hexahedral tree hierarchy, that is, the node is marked. For each marked node, 8 new child nodes are generated, representing the 8 refined sub-hexahedral elements. All the child nodes of the original node are deleted and replaced with the generated 8 new child nodes, and the process is recursive until the preset maximum refinement level is reached.

[0121] Preset error threshold The finite element function is re-solved on the new grid to obtain a new solution uh′. If the corresponding error estimate calculated for each shell element i based on the new solution uh′ is greater than or equal to the error threshold Repeat the encryption and refinement until the calculated error estimate is less than the error threshold Stop the encryption and refinement to obtain the final mesh, which contains n nodes.

[0122] Furthermore, the calculation method of strain energy and residual stress at n nodes includes:

[0123] The finite element function is solved on the final grid to obtain the final solution uh′′ (final nodal displacement field); the final solution uh″ contains the nodal displacement vectors of all shell elements; based on the final solution uh″, the strain field {∈}=uh″·B and stress field {ρ}=uh″·D of each shell element are calculated.

[0124] Calculate the strain energy U of each shell element based on the strain field {∈} and stress field {ρ} of each shell element;

[0125] Where ∫()dV represents the volume integral of the volume V of the shell element; {∈} T is the transpose of the strain field {∈}, [DL(μp,R,θ)] is the nonlinear constitutive matrix, μp is the plastic strain, R is the temperature, and θ is the strain rate; JY is the geometric mapping Jacobian determinant of the shell element; the geometric mapping Jacobian determinant is the determinant value of the Jacobian matrix, which represents the transformation ratio from the parameter space to the physical space. When the shell element undergoes large deformation or rotation, the value of JY will change accordingly, so as to consider the influence of geometric nonlinear effects on the strain energy of the shell element; it is derived from the shape function of the shell element and the coordinates of the nodes.

[0126] [DL(μp,R,θ)]=[De({∈},R)]+[Dp(μp,R,θ)]+[Dvp({∈},μp,R,θ)]; where [De({∈},R)] is the elastic matrix, which describes the elastic behavior of the material, is related to the strain field {∈} and temperature R, and describes the nonlinear elastic properties of the material; [Dp(μp,R,θ)] is the plasticity matrix, which describes the plastic behavior of the material, which is determined by the yield criterion (such as the vonMises criterion) and the associated flow rules; [Dvp({∈},μp,R,θ)] is the viscoplastic matrix, which describes the viscoplastic behavior of the material and is solved using the viscoplastic constitutive model (such as the Perzyna model).

[0127] Based on the plastic matrix and the elastic matrix, the residual stress ER of each shell element is calculated;

[0128] Where {τp} is the plastic strain vector, It is the tensor product operator. For matrices and vectors, the tensor product is equivalent to the matrix product operation.

[0129] The plastic strain vector is solved by a return mapping algorithm, such as the steepest descent method, the tangent method, etc. During the calculation process, it is determined whether the current stress state reaches the yield condition. If the yield condition is not reached, all strains are elastic strains. If the yield condition is reached, the strain needs to be decomposed into an elastic part and a plastic part, and the plastic strain vector is updated.

[0130] For each shell element, the average strain energy and average residual stress of each shell element are calculated based on the calculated strain energy and residual stress. Based on the average strain energy and average residual stress of each shell element, the interpolation function is used at n nodes to perform interpolation operations to obtain the strain energy and residual stress of each node; the interpolation function adopts the same function form as the shape function. For example, for a hexahedral element with 20 nodes, the interpolation function selects a cubic complete polynomial function of 20 nodes.

[0131] It should be noted that for nodes on the structure boundary, since only some units contribute, only the weighted average of the contributing units is calculated during interpolation.

[0132] Furthermore, the method of clustering n nodes includes:

[0133] The strain energy and residual stress of each node constitute a w-dimensional eigenvector respectively; the eigenvectors of n nodes are defined as the initial positions of n fireflies in the w-dimensional space.

[0134] Define the number of clusters K and the clustering objective function fs;

[0135] Among them, C_j is the center of the jth cluster, X_p1 is the position of the p1th firefly (feature vector), and PUX is the inter-cluster penalty function.

[0136] Where q is the index of the cluster (it does not index the same cluster as j), w_jq is the weight coefficient between the jth cluster and the qth cluster; Dd_jq is the Euclidean distance or other Lp distance between the jth cluster and the qth cluster, β1 is the distance order, which is adjusted according to the actual situation to make the distance measurement more flexible; EX is the energy difference penalty function; is the penalty coefficient, the penalty strength can be adjusted, N(j) is the number of nodes in the jth cluster, and N(q) is the number of nodes in the qth cluster.

[0137] Energy difference penalty function EX = (1 + μ2 × ΔE_jq + (1-μ2) × ΔQ_jq); where ΔE_jq is the difference between the average strain energy of the j-th cluster and the q-th cluster, ΔQ_jq is the difference between the average residual stress of the j-th cluster and the q-th cluster; μ2 is a weight adjustment parameter that can adjust the penalty intensity.

[0138] The brightness of fireflies is defined to be inversely proportional to the value of the objective function, that is, fireflies with higher brightness correspond to better clustering effects.

[0139] Calculate the attraction γ_(p1,p2) between the p1th firefly and the p2th firefly;

[0140] γ_(p1,p2)=a×exp(-c′×(RQ(p1,p2)) 2 )×(1+τ2×cos(V(p1,p2)-ε2)); where RQ(p1,p2) is the Euclidean distance between the p1th firefly and the p2th firefly, a is the preset maximum attraction, c′ is the preset light absorption coefficient; τ2 is an adjustment parameter that controls the degree of influence on the attraction, V(p1,p2) is the direction angle between the p1th firefly and the p2th firefly, and ε2 is the preset expected direction angle, which tends to cluster nodes with similar directions together to obtain more compact clusters.

[0141] Move each firefly according to the attractiveness of all other fireflies. The formula for moving is: Among them, α3 is the step size factor, rd is a random number in the interval [0, 1], which is used to increase diversity, X′_p1 is the new position of the p1th firefly, and X_p2 is the position of the p2th firefly before moving.

[0142] The value of the clustering objective function is calculated for the new position of the firefly. If the new position is higher than the value of the clustering objective function at the original position, the firefly is updated to the new position; and the luminous brightness of each firefly is updated according to the new position of the firefly.

[0143] Iterate until the preset number of iterations is met or the clustering objective function converges to obtain the final position of the firefly. According to the final position of the firefly, the adjacent fireflies are grouped into the same cluster. Each cluster corresponds to a node area, thus obtaining m node areas.

[0144] Furthermore, the solution to the reverse bending deformation includes:

[0145] For each node area, the coordinates of all nodes in the node area are extracted as eagles; the parameters of the Hawkeye algorithm are initialized, including the size of the eagle group (the number of eagles M4), the maximum number of iterations and the discovery rate.

[0146] Each eagle is encoded using real number encoding or binary encoding, and the length of the encoding is equal to the number of nodes in the node area multiplied by 3 (x, y, z three coordinates).

[0147] Defining the competition function Among them, w′, w0 and w1 are found to balance the weight parameters, balance the fitting accuracy and smoothness, and obtain a more reasonable node bending deformation solution; En(U_g,ρ′_g) is the deformation state function, U_g is the strain energy of node g in the node area, ρ′_g is the residual stress of node g in the node area, is the curvature of node g in the node area, Cu(g,P′_g) is the plane deviation function, P′_g is the fitting plane for node g in the node area, and for all nodes in a given node area, it is necessary to find an optimal plane to fit these nodes. The fitting plane is the coordinate of the projection point of node g on the plane.

[0148] The competitive function can better balance the fitting accuracy and smoothness and obtain a more reasonable solution for the node inverse bending deformation.

[0149] Deformation state function En(U_g,ρ′_g)=z1×∑ g U_g+(1-z1)×∑ g ρ′_g; where z1 is the deformation equilibrium parameter, which is used to balance the weights of residual stress and strain energy for the function.

[0150] Cu(g,P_g)=x1×(∑ g ||gP′_g||)+(1-x1)×max(||gP′_g||); where,

[0151] x1 is the plane balance parameter, which is used to balance the weights of the two items. The plane deviation function is used to prevent individual nodes from deviating too much.

[0152] Calculate the value of the competition function of each eagle, recorded as the competition value; sort the eagles in the eagle group according to the competition value, and update the position of each eagle by simulating the scanning behavior of the eagle.

[0153] The scanning behavior of the eagle simulates the eagle looking around to find a better position during the hunting process; each eagle will update its position according to the position of the eagle with the highest competition value in the current eagle group and its own current position. The formula for updating the position is:

[0154] X(it+1)=X(it)+AK×cos(2π×r_2)+(2×r_3)×(X_r-X(it))+DPW; wherein X(it+1) is the position of the eagle X at the it+1th iteration, X(it) is the position of the eagle X at the itth iteration, π is the pi, and AK is the scanning radius factor.

[0155] The scanning radius factor AK = 2a1×r_1-a1; where a1 is a constant that controls the scanning radius of the eagle; DPW is a global factor; X_r is the position of the eagle with the highest competition value in the current eagle group.

[0156] Global factor DPW = (2×r_4)×(X_l-X(it))+(2×r_5)×(X_g-X(it)); where r_1, r_2, r_3, r_4 and r_5 are all random update numbers in the interval [0, 1].

[0157] After each position update, the competition value of each eagle in the eagle group is calculated again. According to the preset territory occupancy rate ra, the first ra×M4 eagles are retained and the remaining eagles are eliminated.

[0158] For each eagle retained in each iteration, mutation is performed with a preset probability rp to obtain the new position of the eagle; the mutation formula is:

[0159] X_new=X_old+SF×(X_old-2×r_6×(ub-lb)); where SF is the scaling factor that controls the mutation step size, r_6 is the mutation random number in the interval [0, 1], ub is the upper boundary of the node coordinates, and lb is the lower boundary of the node coordinates; X_new is the position of the eagle after the mutation, and X_old is the position of the eagle before the mutation. If a coordinate component in the new position of a eagle exceeds the preset boundary range, boundary processing of the component is required.

[0160] The eagles at this time constitute a new generation of eagle population. The competition value of the eagles in the new generation population is calculated, and the eagle with the highest competition value is retained as the seed eagle of the next generation of eagle population.

[0161] Repeat the iteration until the number of iterations reaches the defined maximum number of iterations to obtain the final eagle group, calculate the competition value of each eagle in the final eagle group, take the code corresponding to the eagle with the highest competition value as the optimal value, decode the optimal value to obtain the optimal plane fitting coordinates of the nodes in the node area; take the difference between the coordinates of each node in the node area and the optimal plane fitting coordinates as the node's inverse bending deformation; take the average of the inverse bending deformation of each node in the node area as the inverse bending deformation of the node area.

[0162] Furthermore, the tolerance range of the bending deformation of the formed part is pre-set according to actual requirements, for example, ±θ3; wherein θ3 is the maximum allowable bending deformation; and it is determined whether the bending deformation of each node area exceeds the tolerance range.

[0163] For node areas that need to increase stiffness, stiffness enhancement structures such as reinforcement and thickening are added to the finite element model. The form, size, and position of the stiffness enhancement structure need to be designed according to the specific situation. Usually, stiffness enhancement structures are added in or near the node area to improve the overall stiffness of the node area.

[0164] Add the designed stiffness enhancement structure to the finite element analysis model of the formed part, and update the mesh, unit topology, material properties, etc. of the finite element analysis model of the formed part; ensure that the connection and transition between the stiffness enhancement structure and the original model are reasonable and smooth.

[0165] Recalculate the reverse bending deformation of each node area in the finite element analysis model of the formed part, check whether the reverse bending deformation of the node area that originally needs to increase the stiffness has been controlled within the tolerance range. If there are still areas exceeding the tolerance, repeat the steps to continue to optimize the design of the stiffness enhancement structure; when the reverse bending deformation of all node areas is controlled within the tolerance range, output the final finite element analysis model of the formed part and the design scheme of its stiffness enhancement structure.

[0166] This embodiment can improve the forming accuracy of large curvature and shallow drawn parts, and effectively reduce the occurrence of deformation defects; through modeling analysis, the finite element analysis model of the formed part is discretized into multiple nodes, and the strain energy and residual stress at each node are analyzed, which can more accurately predict and evaluate the reverse bending deformation trend of the formed part, and realize the zoning processing of deformation in different areas, thereby improving the flexibility and pertinence of deformation control; secondly, the efficient solution algorithm of the reverse bending deformation improves the calculation efficiency and solution accuracy, and can quickly and accurately solve complex reverse bending deformation, greatly improving the analysis efficiency and shortening the product design cycle; by increasing the stiffness enhancement structure, the stiffness of specific areas is adjusted in a targeted manner to further optimize the deformation control effect. In short, it can not only significantly improve product quality and production efficiency, but also has many advantages such as strong flexibility and excellent computing performance.

[0167] Example 2

[0168] See also Figure 2 As shown, the part not described in detail in this embodiment is described in Example 1, which provides a method for controlling the precision of large curvature and shallow drawn parts, including:

[0169] S1. Obtaining the initial geometric parameters and material mechanical property parameters of the formed part; and establishing a finite element analysis model of the formed part;

[0170] S2, discretizing the finite element analysis model of the formed part into n nodes, and calculating the strain energy and residual stress at the n nodes of the formed part during the unloading process;

[0171] S3. Based on the strain energy and residual stress at the n nodes, the n nodes are clustered to obtain m node regions; the reverse bending deformation of the m node regions is solved by the Hawkeye algorithm;

[0172] S4. Preset the tolerance range, add a stiffness reinforcement structure to the node area where the reverse bending deformation exceeds the tolerance range, and improve the stiffness of the corresponding node area.

[0173] Example 3

[0174] This embodiment discloses an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the operation mode of the method for controlling the precision of large curvature and shallow drawn parts provided above is implemented.

[0175] Since the electronic device introduced in this embodiment is an electronic device used to implement a method for controlling the precision of a large curvature and shallow drawn formed part in the embodiment of this application, based on the method for controlling the precision of a large curvature and shallow drawn formed part introduced in the embodiment of this application, a person skilled in the art can understand the specific implementation of the electronic device of this embodiment and its various variations, so how the electronic device implements the method in the embodiment of this application is not described in detail here. As long as a person skilled in the art implements the electronic device used in the method for controlling the precision of a large curvature and shallow drawn formed part in the embodiment of this application, it belongs to the scope of protection of this application.

[0176] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters and thresholds in the formula are set by technicians in this field according to actual conditions.

[0177] The above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technical users in this technical field, some improvements and modifications without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A precision control system for large curvature and shallow drawn parts, characterized in that: include: Model building module, used to obtain the initial geometric parameters and material mechanical properties parameters of the formed part; And establish the finite element analysis model of the formed parts; Energy analysis model, used to discretize the finite element analysis model of the formed part into n nodes, and calculate the strain energy and residual stress at the n nodes of the formed part during the unloading process; The deformation fitting module clusters n nodes based on the strain energy and residual stress at n nodes to obtain m node areas; the Hawkeye algorithm is used to solve the reverse bending deformation of the m node areas; The detection and adjustment module is used to preset the tolerance range and add a stiffness reinforcement structure to the node area where the reverse bending deformation exceeds the tolerance range. The modules are connected by wire and / or wireless means.

2. The precision control system for large curvature and shallow drawn parts according to claim 1, characterized in that: The initial geometric parameters include length, width, thickness and radius of curvature; the material mechanical property parameters include elastic modulus, Poisson's ratio, yield strength and plasticity parameters.

3. The precision control system for large curvature and shallow drawn parts according to claim 2, characterized in that: The method for establishing the finite element analysis model of the formed part includes: According to the initial geometric parameters of the formed part, a three-dimensional solid model of the formed part is established; based on the mechanical property parameters of the material, material properties are defined on the three-dimensional solid model to obtain a preliminary finite element model; For the preliminary finite element model, its outer surface is extracted as the outer boundary surface, and the inner boundary surface is extracted from the preliminary finite element model by using geometric Boolean operation; based on the outer boundary surface and the inner boundary surface, the intermediate surface is extracted; Predefine the parameter domain, project the extracted intermediate surface into the parameter domain in an equiangular manner, and perform meshing in the parameter domain to generate a mesh topological relationship of rectangular units; The method of performing grid division includes: Nodes are inserted evenly inside and on the boundary of the parameter domain and numbered. Every four adjacent nodes form a rectangular unit. For each rectangular unit, the numbers of the nodes at its four vertices are determined, and the topological relationship of the unit is determined according to the node numbers. Define the physical domain, and construct shell elements in the physical domain according to the mesh topology of the rectangular elements. For each node of the rectangular element in the parameter domain, calculate the corresponding coordinates of the physical domain through the surface equation, and the coordinates of the physical domain are the control point coordinates of the shell element. Use the node numbers in the parameter domain as the node numbers of the shell element in the physical domain. For each shell element i, construct its shape function NF(i); NF(i)=∑ I w_I(α)×N_I(α)+∑ J w_J(β)×N_J(β); where I is the number of the node in the parameter direction α, N_I(α) is the one-dimensional B-spline basis function in the parameter direction α, and w_I(α) is the weight of the node numbered I in the parameter direction α; J is the number of the node in the parameter direction β, N_J(β) is the one-dimensional B-spline basis function in the parameter direction β, and w_J(β) is the weight of the node numbered J in the parameter direction β; Obtain the strain-displacement matrix B and stress-strain matrix D of the shell element; calculate the stiffness matrix K_i of the shell element i based on the strain-displacement matrix B and stress-strain matrix D of the shell element, and calculate the load matrix F_i based on the shape function of the shell element; Stiffness matrix K_i=∫(B T ·D·B)dδ; where B T is the transpose of the strain-displacement matrix B; dδ represents the differential integral on the parameter domain δ; ∫()dδ represents the parameter domain integration on the parameter domain δ; Load matrix F_i=∫((NF(i)) T b)dδ; where b is physical strength, (NF(i)) T is the transpose of NF(i); The stiffness matrices and load matrices of all shell elements are assembled separately to obtain the global stiffness matrix K and the global load matrix F, and the continuity between shell elements is processed to obtain the finite element analysis model of the formed part.

4. The precision control system for large curvature and shallow drawn parts according to claim 3, characterized in that: The method of extracting the intermediate surface includes: Define a signed distance factor d(x,y,z) to represent the shortest distance from any point L in the preliminary finite element model to the inner boundary surface and the outer boundary surface; where x, y and z are the horizontal coordinate, vertical coordinate and vertical coordinate of any point L in the spatial coordinate system of the preliminary finite element model, respectively; Define midsurface deviation Wherein, d_in is the shortest distance from any point L to the inner boundary surface, d_out is the shortest distance from any point L to the outer boundary surface, w_in is the weighting coefficient of the inner boundary surface, and w_out is the weighting coefficient of the outer boundary surface; Discretize the signed distance factor sampling into a three-dimensional voxel grid as the distance field, the three-dimensional voxel grid contains N voxels, each voxel corresponds to a scalar value, that is, the distance field value at the center point of the voxel; For each voxel, check whether there is an intersection point with d(x, y, z) = t between it and the adjacent voxel. If there is an intersection point, use linear interpolation to calculate the coordinates of the intersection point; smoothly connect the calculated coordinates of the intersection point to obtain the intermediate surface.

5. The precision control system for large curvature and shallow drawn parts according to claim 4, characterized in that: The method of discretizing the finite element analysis model of the formed part into n nodes includes: Performing initial mesh division on the finite element analysis model of the formed part to obtain an initial mesh; Define the finite element function, the formula of the finite element function is: K·uh=f; by solving the finite element function on the initial grid, the initial solution uh is obtained; where f is the load vector; Based on the initial solution uh, for each shell element i, the corresponding error estimate is calculated Among them, u is the true solution, C is the elastic matrix, ε(u-uh) is the strain tensor of u-uh; Preset relative error kl, for The shell element i is marked as needing to be encrypted, and the shell element i marked as needing to be encrypted is encrypted and refined to obtain a new mesh; The encryption refinement method includes: The entire computational domain is represented by a hexahedral unit as the root node, and the root node is equally divided and refined to generate 8 child hexahedral units as the child nodes of the root node. Each child node is recursively equally divided and refined until the preset maximum refinement level is reached; that is, a complete hexahedral tree hierarchy is constructed, and each node represents a hexahedral unit; The shell element i marked as needing to be encrypted is mapped to a node in the hexahedral tree hierarchy, that is, the node is marked. For each marked node, 8 new child nodes are generated, representing the 8 refined child hexahedral elements. All the child nodes of the original node are deleted and replaced with the generated 8 new child nodes, and the process is repeated until the preset maximum refinement level is reached. Preset error threshold The finite element function is re-solved on the new mesh to obtain a new solution uh′. If the error estimate corresponding to each shell element i calculated based on the new solution uh′ is greater than or equal to the error threshold Repeat the encryption and refinement until the calculated error estimate is less than the error threshold Stop the encryption and refinement to obtain the final mesh, which contains n nodes.

6. The precision control system for large curvature and shallow drawn parts according to claim 5, characterized in that: The calculation method of the strain energy and residual stress at the n nodes includes: The finite element function is solved on the final mesh to obtain the final solution uh″; the final solution uh″ contains the node displacement vectors of all shell elements; based on the final solution uh″, the strain field {∈}=uh″·B and the stress field {ρ}=uh″·D of each shell element are calculated; Calculate the strain energy U of each shell element based on the strain field {∈} and stress field {ρ} of each shell element; Where ∫()dV represents the volume integral of the volume V of the shell element; {∈} T is the transpose of the strain field {∈}, [DL(μp,R,θ)] is the nonlinear constitutive matrix, μp is the plastic strain, R is the temperature, and θ is the strain rate; JY is the geometric mapping Jacobian of the shell element; [DL(μp,R,θ)]=[De({∈},R)]+[Dp(μp,R,θ)]+[Dvp({∈},μp,R,θ)]; where [De({∈},R)] is the elastic matrix; [Dp(μp,R,θ)] is the plastic matrix; [Dvp({∈},μp,R,θ)] is the viscoplastic matrix; Based on the plastic matrix and the elastic matrix, the residual stress ER of each shell element is calculated; Where {τp} is the plastic strain vector, is the tensor product operator; For each shell element, the average strain energy and average residual stress of each shell element are calculated based on the calculated strain energy and residual stress. Based on the average strain energy and average residual stress of each shell element, the interpolation function is used at n nodes to perform interpolation operations to obtain the strain energy and residual stress of each node.

7. The precision control system for large curvature and shallow drawn parts according to claim 6, characterized in that: The method of clustering n nodes includes: The strain energy and residual stress of each node constitute a w-dimensional feature vector respectively; the feature vectors of n nodes are defined as the initial positions of n fireflies in the w-dimensional space; Define the number of clusters K and the clustering objective function fs; Among them, C_j is the center of the jth cluster, X_p1 is the position of the p1th firefly, and PUX is the inter-cluster penalty function; Where q is the index of the cluster, w_jq is the weight coefficient between the jth cluster and the qth cluster; Dd_jq is the Euclidean distance between the jth cluster and the qth cluster, β1 is the distance order; EX is the energy difference penalty function; θ1 is the penalty coefficient, N(j) is the number of nodes in the jth cluster, and N(q) is the number of nodes in the qth cluster; Energy difference penalty function EX = (1 + μ2 × ΔE_jq + (1-μ2) × ΔQ_jq); where ΔE_jq is the difference between the average strain energy of the jth cluster and the qth cluster, ΔQ_jq is the difference between the average residual stress of the jth cluster and the qth cluster; μ2 is the weight adjustment parameter; Define that the brightness of a firefly is inversely proportional to the value of the objective function; calculate the attraction γ_(p1,p2) between the p1th firefly and the p2th firefly; γ_(p1,p2)=a×exp(-c′×(RQ(p1,p2)) 2 )×(1+τ2×cos(V(p1,p2)-ε2)); where RQ(p1,p2) is the Euclidean distance between the p1th firefly and the p2th firefly, a is the preset maximum attraction, c′ is the preset light absorption coefficient; τ2 is the adjustment parameter, V(p1,p2) is the direction angle between the p1th firefly and the p2th firefly, and ε2 is the preset expected direction angle; Move each firefly according to the attractiveness of all other fireflies. The formula for moving is: Among them, α3 is the step size factor, rd is a random number in the interval [0, 1], X′_p1 is the new position of the p1th firefly, and X_p2 is the position of the p2th firefly before moving; Calculate the value of the clustering objective function for the new position of the firefly. If the new position is higher than the value of the clustering objective function at the original position, update the firefly to the new position; and update the brightness of each firefly according to the new position of the firefly; Iterate until the preset number of iterations is met to obtain the final position of the firefly. According to the final position of the firefly, the adjacent fireflies are grouped into the same cluster. Each cluster corresponds to a node area, thus obtaining m node areas.

8. The precision control system for large curvature and shallow drawn parts according to claim 7, characterized in that: The method for calculating the reverse bending deformation includes: For each node area, the coordinates of all nodes in the node area are extracted as eagles; the parameters of the eagle-eye algorithm are initialized, including the scale of the eagle group, that is, the number of eagles M4, the maximum number of iterations and the discovery rate; Encode each eagle using real number coding or binary coding, and the length of the code is equal to the number of nodes in the node area multiplied by 3; Define the competition function BN; calculate the value of the competition function of each eagle, recorded as the competition value; sort the eagles in the eagle group according to the competition value, and update the position of each eagle according to the position of the eagle with the highest competition value in the current eagle group and its own current position; After each position update, the competition value of each eagle in the eagle group is calculated again. According to the preset territory occupancy rate ra, the first ra×M4 eagles are retained, and the rest are eliminated; For each eagle retained in each iteration, mutation is performed with a preset probability rp to obtain the new position of the eagle; At this time, the eagles form a new generation of eagles. The competition value of the eagles in the new generation of population is calculated, and the eagle with the highest competition value is retained as the seed eagle of the next generation of eagles. Repeat the iteration until the number of iterations reaches the defined maximum number of iterations to obtain the final eagle group, calculate the competition value of each eagle in the final eagle group, take the code corresponding to the eagle with the highest competition value as the optimal value, decode the optimal value to obtain the optimal plane fitting coordinates of the nodes in the node area; take the difference between the coordinates of each node in the node area and the optimal plane fitting coordinates as the node's inverse bending deformation; take the average of the inverse bending deformation of each node in the node area as the inverse bending deformation of the node area.

9. The precision control system for large curvature and shallow drawn parts according to claim 8, characterized in that: The formula of the competition function BN is: BN=w′×En(U_g,ρ′_g)+w0×∑ g θ′_g+w1×Cu(g,P_g); where w′, w0 and w1 are found equilibrium weight parameters; En(U_g,ρ′_g) is the deformation state function, U_g is the strain energy of node g in the node region, ρ′_g is the residual stress of node g in the node region, θ′_g is the curvature of node g in the node region, Cu(g,P′_g) is the plane deviation function, and P′_g is the fitting plane for node g in the node region; Deformation state function En(U_g,ρ′_g)=z1×∑ g U_g+(1-z1)×∑ g ρ′_g; where z1 is the deformation equilibrium parameter; Cu(g,P_g)=x1×(∑ g ||gP′_g||)+(1-x1)×max(||gP′_g||); where x1 is the plane balance parameter; The formula for updating the position is: X(it+1)=X(it)+AK×cos(2π×r_2)+(2×r_3)×(X_r-X(it))+DPW; where X(it+1) is the position of the eagle X at the it+1th iteration, X(it) is the position of the eagle X at the itth iteration, π is pi, and AK is the sweep radius factor; The scanning radius factor AK = 2a1 × r_1-a1; where a1 is a constant, DPW is a global factor; X_r is the position of the eagle with the highest competition value in the current eagle group; Global factor DPW = (2×r_4)×(X_l-X(it))+(2×r_5)×(X_g-X(it)); where r_1, r_2, r_3, r_4 and r_5 are all random update numbers in the interval [0, 1]; The formula for performing mutation is: X_new=X_old+SF×(X_old-2×r_6×(ub-lb)); where SF is the scaling factor, r_6 is the mutation random number in the interval [0, 1], ub is the upper boundary of the node's coordinates, and lb is the lower boundary of the node's coordinates; X_new is the position of the eagle after the mutation, and X_old is the position of the eagle before the mutation.

10. A method for controlling the precision of large curvature and shallow drawn parts, which is implemented based on the precision control system for large curvature and shallow drawn parts according to any one of claims 1 to 9, characterized in that: include: S1. Obtaining initial geometric parameters and material mechanical property parameters of the formed part; And establish the finite element analysis model of the formed parts; S2, discretizing the finite element analysis model of the formed part into n nodes, and calculating the strain energy and residual stress at the n nodes of the formed part during the unloading process; S3. Based on the strain energy and residual stress at the n nodes, the n nodes are clustered to obtain m node regions; the reverse bending deformation of the m node regions is solved by the Hawkeye algorithm; S4. Preset the tolerance range and increase the stiffness of the node area where the reverse bending deformation exceeds the tolerance range to strengthen the structure.

Images