An accuracy control system and method for a large curvature and shallow drawing formed part

Abstract

The present invention belongs to the technical field of workpiece precision control. The present invention discloses a precision control system and method for large-curvature and shallow-drawing formed parts, including: 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; discretizing the finite element analysis model of the formed part into n nodes, and calculating the stress distribution at the n nodes during the unloading process of the formed part, the strain energy generated by springback, and the springback bending moment formed thereby; based on the stress characteristics at the n nodes, clustering the n nodes by characteristics to obtain m node regions; using the eagle-eye algorithm to solve the overall springback deformation amount and springback strain energy of the m node regions; presetting a tolerance range, and for the node regions where the springback deformation amount exceeds the tolerance range, increasing the springback bending stiffness of the corresponding regions based on the characteristic geometric shape and offsetting the deformation amount generated by the springback bending moment, effectively improving the dimensional accuracy and shape quality of large-curvature and shallow-drawing formed parts.

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-drawing formed parts. Background Art

[0002] The patent with the application publication number CN109513931A discloses a control method for residual thermal stress and its induced deformation in additive manufacturing, including the following steps: obtaining the residual stress distribution and stress values of the formed part by experimental or numerical simulation methods; designing a continuous or discontinuous flexible base structure according to the specific structure form of the formed part, 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 base material. By designing multiple groups of flexible base structures, the zonal release and active control of residual stress based on the flexible structure are realized. During the additive forming process, gradually change the characteristic size of the flexible structure to reduce the structure flexibility coefficient and obtain a gradient flexible structure to meet the service load requirements of the component. Through structural deformation, spatially discrete and release the residual stress, and realize the discrete control of macroscopic thermal stress in space within the allowable deformation range.

[0003] During the metal sheet forming process, the formed part will produce reverse bending deformation during unloading, and 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 difficult control of the dimensional accuracy and shape quality of the formed part; different regions of the formed part may generate different degrees of strain energy and residual stress during the forming process, and existing methods often treat the entire formed part as a whole and cannot optimize different regions specifically, 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 extremely prone to large reverse bending deformation after stamping, resulting in problems such as interference with other components and excessive gaps during assembly; if the reverse bending deformation cannot be effectively controlled, it will seriously affect the overall quality and performance of the automobile.

[0004] In view of this, the present invention proposes a precision control system and method for large-curvature and shallow-drawing formed 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 object, the present invention provides the following technical solution: A precision control system for large-curvature and shallow-drawing formed parts, comprising: a model construction module, configured to obtain the initial geometric parameters and material mechanical property parameters of the formed part; and establish a finite element analysis model of the formed part;

[0006] An energy analysis model is 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 during the unloading process of the formed part;

[0007] A deformation amount fitting module clusters the n nodes based on the strain energy and residual stress at the n nodes to obtain m node regions; and uses the eagle eye algorithm to solve the reverse bending deformation amount of the m node regions;

[0008] A detection and adjustment module is used to preset a tolerance range, and add a stiffness enhancement structure to the node regions where the reverse bending deformation amount exceeds the tolerance range. Each module is connected by wired 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 plastic parameters.

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

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

[0012] For the preliminary finite element model, extract its outer surface as the outer boundary surface, and use geometric Boolean operations to extract the inner boundary surface from the preliminary finite element model; based on the outer boundary surface and the inner boundary surface, extract the intermediate surface;

[0013] Predefine the parameter domain, project the extracted intermediate surface onto the parameter domain by equiangular projection, and perform mesh division in the parameter domain to generate the mesh topology relationship of rectangular elements;

[0014] The method of performing mesh division includes:

[0015] Uniformly insert nodes on the inside and boundary of the parameter domain and number the nodes. Each adjacent four nodes form a rectangular element. For each rectangular element, determine the numbers of the nodes at its four vertices, and determine the topological relationship of the element according to the numbers of the nodes;

[0016] Define the physical domain. According to the mesh topology relationship of the rectangular elements, construct shell elements in the physical domain; for the nodes of each 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 numbers of the nodes in the parameter domain as the numbers of the nodes of the shell elements 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] The 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 quantity on the parameter domain δ; ∫()dδ represents the parameter domain integration over the parameter domain δ;

[0021] The load matrix F_i = ∫((NF(i)) T ·b)dδ; where b is the body force, and (NF(i)) T is the transpose of NF(i);

[0022] Assemble the stiffness matrices and load matrices of all shell elements respectively to obtain the global stiffness matrix K and global load matrix F, and process the continuity between shell elements to obtain the finite element analysis model of the formed part.

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

[0024] Define the signed distance factor d(x, y, z), which represents 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 abscissa, ordinate, and vertical coordinate of any point L in the space coordinate system where the preliminary finite element model is located, respectively;

[0025] Define the intermediate surface deviation where 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;

[0026] Sample and discretize the signed distance factor into a three-dimensional voxel grid as the 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;

[0027] For each voxel, check whether there is an intersection point where d(x, y, z) = t with adjacent voxels. If there is an intersection point, calculate the coordinates of the intersection point using linear interpolation; smooth-connect the calculated coordinates of the intersection points 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] Perform an initial mesh division on the finite element analysis model of the formed part to obtain an initial mesh;

[0030] Define a finite element function, and the formula of the finite element function is: K·uh = f; solve the finite element function on the initial mesh to obtain an initial solution uh; where f is the load vector;

[0031] Based on the initial solution uh, calculate the corresponding error estimator for each shell element i

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

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

[0034] The method of performing encryption and refinement includes:

[0035] Represent the entire computational domain with a hexahedral element as the root node, perform equal division and refinement on the root node to generate 8 sub - hexahedral elements as the child nodes of the root node, and recursively perform equal division and refinement on each child node until the preset maximum refinement level is reached; that is, construct a complete hexahedral tree hierarchical structure, and each node represents a hexahedral element;

[0036] Correspond the marked shell element i to the nodes in the hexahedral tree hierarchical structure, that is, the nodes are marked. For each marked node, generate 8 new child nodes, representing 8 refined sub - hexahedral elements, delete all the child nodes of the original node, and replace them with the generated 8 new child nodes, and recurse until the preset maximum refinement level is reached;

[0037] Preset an error threshold And re - solve the finite element function on the new mesh to obtain a new solution uh′. If the corresponding error estimator calculated for each shell element i based on the new solution uh′ is greater than or equal to the error threshold Then repeat the encryption refinement until the calculated error estimator is less than the error threshold. Stop the encryption refinement to obtain the final mesh, which contains n nodes.

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

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

[0040] Based on the strain field {∈} and stress field {ρ} of each shell element, calculate the strain energy U 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 non-linear constitutive matrix, μp is the plastic strain, R is the temperature, θ is the strain rate; JY is the geometric mapping Jacobian determinant 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 elastic matrix, calculate the residual stress ER of each shell element;

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

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

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

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

[0048] Define the number of clusters \(K\) and the clustering objective function \(f_s\);

[0049] where \(C_j\) is the center of the \(j\)-th cluster, \(X_{p1}\) is the position of the \(p1\)-th firefly, and \(PUX\) is the inter-cluster penalty function;

[0050] where \(q\) is the clustering index, \(w_{jq}\) is the weight coefficient between the \(j\)-th cluster and the \(q\)-th cluster; \(Dd_{jq}\) is the Euclidean distance between the \(j\)-th cluster and the \(q\)-th cluster, \(\beta1\) is the distance power; \(EX\) is the energy difference penalty function; \(\mu1\) is the penalty coefficient, \(N(j)\) is the number of nodes in the \(j\)-th cluster, and \(N(q)\) is the number of nodes in the \(q\)-th cluster;

[0051] The energy difference penalty function \(EX=(1 + \mu2\times\Delta E_{jq}+(1 - \mu2)\times\Delta Q_{jq})\); where \(\Delta E_{jq}\) is the difference in average strain energy between the \(j\)-th cluster and the \(q\)-th cluster, \(\Delta Q_{jq}\) is the difference in average residual stress between the \(j\)-th cluster and the \(q\)-th cluster; \(\mu2\) is the weight adjustment parameter;

[0052] Define that the luminous intensity of the firefly is inversely proportional to the value of the objective function; calculate the attractiveness \(\gamma_{(p1,p2)}\) between the \(p1\)-th firefly and the \(p2\)-th firefly;

[0053] \(\gamma_{(p1,p2)} = a\times\exp(-c'\times(RQ(p1,p2)) 2 )\times(1 + \tau2\times\cos(V(p1,p2)-\varepsilon2))\); where \(RQ(p1,p2)\) is the Euclidean distance between the \(p1\)-th firefly and the \(p2\)-th firefly, \(a\) is the preset maximum attractiveness, \(c'\) is the preset light absorption coefficient; \(\tau2\) is the adjustment parameter, \(V(p1,p2)\) is the direction angle between the \(p1\)-th firefly and the \(p2\)-th firefly, and \(\varepsilon2\) is the preset desired direction angle;

[0054] Move each firefly according to the attractiveness of all other fireflies, and the movement formula is: where \(\alpha3\) is the step size factor, \(rd\) is a random number in the range \([0, 1]\), \(X'_{p1}\) is the new position of the \(p1\)-th firefly, and \(X_{p2}\) is the position of the \(p2\)-th firefly before movement;

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

[0056] Iterate until the preset number of iterations is satisfied to obtain the final positions of the fireflies. According to the final positions of the fireflies, group adjacent fireflies into the same cluster, and each cluster corresponds to a node region, thereby obtaining m node regions.

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

[0058] For each node region, extract the coordinates of all nodes within the node region as the eagles; initialize the parameters of the eagle eye algorithm, where the parameters include the scale of the eagle group, i.e., the number of eagles M4, the maximum number of iterations, and the discovery rate.

[0059] Encode each eagle, and the encoding method adopts real number encoding or binary encoding, and the encoding length is equal to the number of nodes within the node region multiplied by 3.

[0060] Define the competition function BN; calculate the value of the competition function for each eagle, denoted as the competition value; sort the eagles within 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, calculate the competition value of each eagle within the eagle group again. According to the preset territory occupancy rate ra, retain the first ra×M4 eagles, and the remaining eagles are eliminated.

[0062] For the eagles retained in each iteration, mutate with a preset probability rp to obtain the new positions of the eagles.

[0063] The eagles at this time form a new generation of eagle group. Calculate the competition value of the eagles within the new generation population, and retain the eagle with the highest competition value as the seed eagle for the next generation of eagle group.

[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 within the final eagle group, take the encoding corresponding to the eagle with the highest competition value as the optimal value, decode the optimal value, and obtain the optimal plane fitting coordinates of the nodes within the node region; take the difference between the coordinates of each node within the node region and the optimal plane fitting coordinates as the reverse bending deformation amount of the node; take the average value of the reverse bending deformation amounts of each node within the node region as the reverse bending deformation amount of the node region.

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

[0066] where w′, w0

[0067] w1 is the discovered balance weight parameter; En(U_g, ρ′_g) is the deformation state function, where U_g is the strain energy of node g in the node area and ρ′_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, and P′_g is the fitting plane for node g in the node area.

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

[0069] Cu(g, P_g) = x1 × (∑ g ||g - P′_g||) + (1 - x1) × max(||g - P′_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 eagle X at the (it + 1)-th iteration, X(it) is the position of eagle X at the it-th iteration, π is the pi, and AK is the saccade radius factor.

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

[0074] The 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 randomly updated numbers in the interval [0, 1].

[0075] The formula for mutation is:

[0076] X_new = X_old + SF × (X_old - 2 × r_6 × (ub - lb)); where SF is the scaling factor, r_6 is a 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 mutated eagle, and X_old is the position of the eagle before mutation.

[0077] A precision control method for large curvature and shallow drawing formed parts, which is realized based on the described precision control system for large curvature and shallow drawing formed parts, includes: S1. Obtain the initial geometric parameters and material mechanical property parameters of the formed part; and establish a finite element analysis model of the formed part;

[0078] S2. 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 during the unloading process of the formed part;

[0079] S3. Based on the strain energy and residual stress at the n nodes, cluster the n nodes to obtain m node regions; use the eagle eye algorithm to solve the reverse bending deformation amount of the m node regions;

[0080] S4. Preset a tolerance range, and add a stiffness enhancement structure to the node regions where the reverse bending deformation amount exceeds the tolerance range.

[0081] Technical effects and advantages of the precision control system and method for large curvature and shallow drawing formed parts of the present invention:

[0082] The present invention can improve the forming precision of large curvature and shallow drawing formed parts, and effectively reduce the occurrence of deformation defects; through modeling and 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 realizes the partition processing of the deformation in different regions, improving the flexibility and pertinence of deformation control; secondly, the efficient solution algorithm for the reverse bending deformation amount improves the calculation efficiency and solution accuracy, can quickly and accurately solve complex reverse bending deformation amounts, greatly improves the analysis efficiency, and shortens the product design cycle; by adding a stiffness enhancement structure, the stiffness of specific regions is adjusted specifically to further optimize the deformation control effect. In short, it can not only significantly improve the product quality and production efficiency, but also has many advantages such as strong flexibility and excellent calculation performance. Brief Description of the Drawings

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

[0084] Figure 2 It is a schematic diagram of the precision control method for large curvature and shallow drawing formed parts of the present invention. Detailed Embodiments

[0085] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0086] Embodiment 1

[0087] Please refer to Figure 1 As shown, a precision control system for a large curvature and shallow drawing forming part in this embodiment includes:

[0088] A model construction module, which is used to obtain the initial geometric parameters and material mechanical property parameters of the forming part; and establish a finite element analysis model of the forming part;

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

[0090] A deformation amount fitting module, based on the strain energy and residual stress at the n nodes, clusters the n nodes to obtain m node regions; uses the eagle eye algorithm to solve the reverse bending deformation amount of the m node regions;

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

[0092] Further, 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 plastic parameters.

[0093] The initial geometric parameters are directly measured and obtained; the material mechanical property parameters are given in relevant material manuals and can be directly consulted and obtained; if the required parameters cannot be directly consulted, they are indirectly obtained through reverse analysis; for example, given the stress-strain response of a material under certain specific working conditions, the combination of material parameters that satisfies this response is deduced.

[0094] Further, the establishment method of the finite element analysis model of the forming part includes:

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

[0096] Extract the middle surface from the preliminary finite element model and generate a triangular mesh surface. Specifically, for the preliminary finite element model, extract its outer surface as the outer boundary surface, and use geometric Boolean operations to extract the inner boundary surface from the preliminary finite element model.

[0097] Define the signed distance factor d(x, y, z), which represents the shortest distance from any point L inside 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. Here, x, y, and z are the abscissa, ordinate, and vertical coordinate of any point L in the space coordinate system where the preliminary finite element model is located, respectively.

[0098] Define the middle surface deviation where 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. Take the middle surface deviation as the objective function, and use the least squares method to solve the value of the weighting coefficient that can minimize the objective function to obtain the optimal weighting coefficient.

[0099] Sample and discretize the signed distance factor into a three-dimensional voxel grid as the 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 where d(x, y, z) = t between it and its adjacent voxels. 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 points to obtain the middle surface.

[0101] Pre-define a parameter domain (usually choose a rectangular domain or a triangular domain), project the extracted middle surface isometrically into the parameter domain, and perform mesh division in the parameter domain to generate the mesh topology relationship of rectangular elements. The isometric projection aims to keep the angles of the surface unchanged in the parameter domain.

[0102] The methods of performing mesh division include:

[0103] Uniformly insert nodes inside and on the boundary of the parameter domain and number the nodes. Each adjacent four nodes form a rectangular element. For each rectangular element, determine the numbers of the nodes at its four vertices, and determine the topological relationship of the element according to the numbers of the nodes.

[0104] That is, (m + 1)×(n + 1) nodes are generated within the parameter domain, as well as m×n rectangular elements; each element is determined by the numbers of four nodes, thereby constructing the topological relationship of the elements.

[0105] Define the physical domain (the actual three-dimensional space domain, which corresponds one-to-one with the parameter domain). According to the grid topological relationship of the rectangular elements, shell elements are constructed in the physical domain; specifically, for the nodes of each rectangular element within the parameter domain, the corresponding coordinates in the physical domain are calculated through the surface equation, and these coordinates in the physical domain are the control point coordinates of the shell elements; the numbers of the nodes within the parameter domain are used as the numbers of the nodes of the shell elements within 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 α, 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 to enable the shape function to 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 β, 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, usually taken as orthogonal. By taking different values of α and β within the domain of definition, a series of parameter points can be generated within the parameter domain; the shape function NF(i) has stronger geometric expression ability and at the same time retains piecewise continuity and affine invariance.

[0108] Obtain the strain-displacement matrix B and stress-strain matrix D of the shell element. The strain-displacement matrix and stress-strain matrix are derived through the shape function of the shell element and the material constitutive relationship; based on 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, and calculate the load matrix F_i 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 quantity on the parameter domain δ; ∫()dδ represents the parameter domain integration over the parameter domain δ.

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

[0111] Assemble the stiffness matrices and load matrices of all shell elements respectively to obtain the global stiffness matrix K and the global load matrix F, and handle the continuity between shell elements (two adjacent shell elements share a boundary, ensuring that the displacements on this boundary are 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] Perform an initial mesh division on the finite element analysis model of the formed part to obtain an initial mesh; the method of initial mesh division can adopt the mapping method, the feedforward method, or the feedforward mapping method; a set of nodes generated during the initial mesh division is relatively evenly distributed.

[0114] Define the finite element function, and the formula of the finite element function is: K·uh = f; solve the finite element function on the initial mesh to obtain the initial solution uh (the nodal displacement field); where f is the load vector (the external load acting on the structure).

[0115] Based on the initial solution uh, calculate the corresponding error estimator for each shell element i

[0116] where 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 approximation is constructed using the recovery technique or the projection technique.

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

[0118] The method of performing encryption and refinement includes:

[0119] Represent the entire computational domain (the boundary of the finite element analysis model of the formed part describes the geometric shape and range of the entire computational domain) with a hexahedron element as the root node, perform equal division and refinement on the root node to generate 8 sub - hexahedron elements as the child nodes of the root node, and recursively perform equal division and refinement on each child node until the preset maximum refinement level is reached; that is, construct a complete hexahedron tree hierarchical structure, and each node (the root node and the child nodes) represents a hexahedron element.

[0120] Correspond the shell element \(i\) marked as needing encryption to nodes on the hexahedron tree hierarchy, that is, the nodes are marked. For each marked node, generate 8 new child nodes, representing 8 refined sub - hexahedron elements. Delete all the child nodes of the original node and replace them with the 8 newly generated child nodes. Recursively repeat this process until the preset maximum refinement level is reached.

[0121] Preset error threshold And re - solve the finite - element function on the new mesh to obtain a new solution \(u_h'\). If the corresponding error estimator calculated for each shell element \(i\) based on the new solution \(u_h'\) is greater than or equal to the error threshold Then repeat the encryption and refinement process; until the calculated error estimator is less than the error threshold Stop the encryption and refinement process to obtain the final mesh, and the final mesh contains \(n\) nodes.

[0122] Furthermore, the calculation methods of the strain energy and residual stress at the \(n\) nodes include:

[0123] Solve the finite - element function on the final mesh to obtain the final solution \(u_h''\) (the final nodal displacement field); the final solution \(u_h''\) contains the nodal displacement vectors of all shell elements. Based on the final solution \(u_h''\), calculate the strain field \(\{\epsilon\}=u_h''\cdot B\) and stress field \(\{\sigma\}=u_h''\cdot D\) for each shell element.

[0124] Based on the strain field \(\{\epsilon\}\) and stress field \(\{\sigma\}\) of each shell element, calculate the strain energy \(U\) of each shell element;

[0125] where, \(\int()dV\) represents the volume integral over the volume \(V\) of the shell element; \(\{\epsilon\}^T\) T is the transpose of the strain field \(\{\epsilon\}\), \([D_L(\mu_p,R,\theta)]\) is the non - linear constitutive matrix, \(\mu_p\) is the plastic strain, \(R\) is the temperature, \(\theta\) is the strain rate; \(J_Y\) 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 parametric space to the physical space. When the shell element undergoes large deformations and large rotations, the value of \(J_Y\) will change accordingly, so as to be able to consider the influence of geometric non - linear effects on the strain energy of the shell element; it is derived from the shape functions 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, describing the elastic behavior of the material, related to the strain field {∈} and temperature R, and describing the nonlinear elastic properties of the material; [Dp(μp, R, θ)] is the plastic matrix, describing the plastic behavior of the material, determined by the yield criterion (such as the von Mises criterion) and the associated flow rule; [Dvp({∈}, μp, R, θ)] is the viscoplastic matrix, describing the viscoplastic behavior of the material, and solved using a 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] Among them, {τp} is the plastic strain vector, is the tensor product operator. For matrices and vectors, the tensor product is equivalent to the matrix multiplication operation.

[0129] The plastic strain vector is solved by the return mapping algorithm, such as the steepest descent method, the tangent method, etc.; during the calculation process, it is judged whether the current stress state reaches the yield condition. If it does not reach the yield, all strains are elastic strains; if it reaches the yield, 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, based on the calculated strain energy and residual stress, the average strain energy and average residual stress of each shell element are calculated. Based on the average strain energy and average residual stress of each shell element, interpolation operations are performed at n nodes using an interpolation function to obtain the strain energy and residual stress of each node; the interpolation function adopts the same functional form as the shape function. For example, for a 20-node hexahedron element, the interpolation function selects a cubic complete polynomial function of 20 nodes.

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

[0132] Furthermore, the ways to cluster the n nodes include:

[0133] The strain energy and residual stress of each node are respectively formed into a w-dimensional feature vector; the feature vectors of the 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 \(j\)-th cluster, \(X_{p1}\) is the position of the \(p1\)-th firefly (feature vector), and PUX is the inter-cluster penalty function.

[0136] Among them, \(q\) is the index of the cluster (not indexing the same cluster as \(j\) at the same time), \(w_{jq}\) is the weight coefficient between the \(j\)-th cluster and the \(q\)-th cluster; \(Dd_{jq}\) is the Euclidean distance or other \(Lp\) distance between the \(j\)-th cluster and the \(q\)-th cluster, \(\beta1\) is the distance exponent, which is adjusted according to the actual situation to make the distance metric more flexible; EX is the energy difference penalty function; \(\lambda\) is the penalty coefficient, which can adjust the penalty intensity, \(N(j)\) is the number of nodes in the \(j\)-th cluster, and \(N(q)\) is the number of nodes in the \(q\)-th cluster.

[0137] The energy difference penalty function \(EX=(1 + \mu2\times\Delta E_{jq}+(1 - \mu2)\times\Delta Q_{jq})\); where, \(\Delta E_{jq}\) is the difference in the average strain energy between the \(j\)-th cluster and the \(q\)-th cluster, \(\Delta Q_{jq}\) is the difference in the average residual stress between the \(j\)-th cluster and the \(q\)-th cluster; \(\mu2\) is the weight adjustment parameter, which can adjust the penalty intensity.

[0138] It is defined that the luminous brightness of the firefly is inversely proportional to the value of the objective function, that is, the firefly with high brightness corresponds to a better clustering effect.

[0139] Calculate the attractiveness \(\gamma_{(p1,p2)}\) between the \(p1\)-th firefly and the \(p2\)-th firefly;

[0140] \(\gamma_{(p1,p2)} = a\times\exp(-c'\times(RQ(p1,p2)) 2 )\times(1 + \tau2\times\cos(V(p1,p2)-\varepsilon2))\); where, \(RQ(p1,p2)\) is the Euclidean distance between the \(p1\)-th firefly and the \(p2\)-th firefly, \(a\) is the preset maximum attractiveness, \(c'\) is the preset light absorption coefficient; \(\tau2\) is the adjustment parameter, which controls the influence degree on the attractiveness, \(V(p1,p2)\) is the direction angle between the \(p1\)-th firefly and the \(p2\)-th firefly, and \(\varepsilon2\) is the preset expected direction angle, so as to tend to cluster nodes with similar directions together, thereby obtaining a more compact cluster.

[0141] Move each firefly according to the attractiveness of all other fireflies, and the movement formula is: Among them, \(\alpha3\) is the step size factor, \(rd\) is a random number in the interval \([0, 1]\) used to increase diversity, \(X'_{p1}\) is the new position of the \(p1\)-th firefly, and \(X_{p2}\) is the position of the \(p2\)-th firefly before moving.

[0142] Calculate the value of the clustering objective function for the new position of the firefly. If the value of the clustering objective function for the new position is higher than that of the original position, update the position of the firefly to the new position; and update the luminous intensity of each firefly according to the new position of the firefly.

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

[0144] Further, the solution method of the reverse bending deformation amount includes:

[0145] For each node area, extract the coordinates of all nodes in the node area as eagles; initialize the parameters of the eagle eye algorithm, and the parameters include the scale of the eagle group (the number of eagles M4), the maximum number of iterations, and the discovery rate.

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

[0147] Define the competition function Among them, w′, w0, and w1 are discovery balance weight parameters, which balance the fitting accuracy and smoothness to obtain a more reasonable solution for the reverse bending deformation amount of the node; En(U_g,ρ′_g) is the deformation state function, U_g is the strain energy of node g in the node area, and ρ′_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. For all nodes in the given node area, an optimal plane needs to be found to fit these nodes, and the fitting plane is the coordinate of the projection point of node g on this plane.

[0148] The competition function can better balance the fitting accuracy and smoothness to obtain a more reasonable solution for the reverse bending deformation amount of the node.

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

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

[0151] x1 is a planar balance parameter used to balance the weights of two items, and the planar deviation function is used to prevent the situation where individual nodes deviate too much.

[0152] Calculate the value of the competition function for each eagle, denoted 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 eagles.

[0153] The scanning behavior of the eagles simulates the eagles looking around during the predation process to find a better position; 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; where X(it + 1) is the position of eagle X at the (it + 1)-th iteration, X(it) is the position of eagle X at the it-th 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 eagles, DPW is the global factor, and X_r is the position of the eagle with the highest competition value in the current eagle group.

[0156] The 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 randomly updated numbers within the interval [0, 1].

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

[0158] For the eagles retained in each iteration, mutate with a preset probability rp to obtain the new position of the eagles. The formula for mutation 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 a random mutation number within the interval [0, 1], ub is the upper boundary of the node coordinates, lb is the lower boundary of the node coordinates; X_new is the position of the mutated eagle, X_old is the position of the eagle before mutation. If a certain coordinate component in the new position of an eagle exceeds the preset boundary range, boundary processing needs to be performed on this component.

[0160] At this time, the eagles only form a new generation of eagle group. Calculate the competition values of the eagles within the new generation of population, and retain the eagle with the highest competition value as the seed eagle of the next generation of eagle group.

[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 values of each eagle within the final eagle group. Take the encoding corresponding to the eagle with the highest competition value as the optimal value, and decode the optimal value to obtain the optimal plane fitting coordinates of the nodes within the node area; take the difference between the coordinates of each node within the node area and the optimal plane fitting coordinates as the reverse bending deformation amount of the node; take the average value of the reverse bending deformation amounts of each node within the node area as the reverse bending deformation amount of the node area.

[0162] Furthermore, according to actual requirements, the tolerance range is preset for the reverse bending deformation amount of the formed part, such as ±θ3; where θ3 is the maximum allowable reverse bending deformation amount; determine whether the reverse bending deformation amount of each node area exceeds the tolerance range.

[0163] For the node areas that need to increase stiffness, add stiffness enhancement structures in the finite element model, such as stiffening, thickening, etc. The form, size, position, etc. of the stiffness enhancement structures need to be designed according to specific circumstances. Usually, add stiffness enhancement structures within or near the node areas to improve the overall stiffness of the node areas.

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

[0165] Recalculate the reverse bending deformation amounts of each node area within the finite element analysis model of the formed part, and check whether the reverse bending deformation amounts of the original node areas that need to increase stiffness have been controlled within the tolerance range. If there are still areas exceeding the tolerance, repeat the steps to continue optimizing the design of the stiffness enhancement structures; when the reverse bending deformation amounts of all node areas are controlled within the tolerance range, output the final finite element analysis model of the formed part and the design scheme of its stiffness enhancement structures.

[0166] This embodiment can improve the forming accuracy of large-curvature and shallow-drawing formed parts, effectively reducing 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 realizes the zonal treatment of deformations in different regions, improving the flexibility and pertinence of deformation control; secondly, the efficient solution algorithm for the reverse bending deformation amount improves the calculation efficiency and solution accuracy, can quickly and accurately solve the complex reverse bending deformation amount, greatly improves the analysis efficiency, and shortens the product design cycle; by adding a stiffness enhancement structure to specifically adjust the stiffness of a specific region, the deformation control effect is further optimized. In short, it can not only significantly improve the product quality and production efficiency, but also has many advantages such as strong flexibility and excellent computing performance.

[0167] Embodiment 2

[0168] Please refer to Figure 2 As shown, for the parts not described in detail in this embodiment, refer to the description content of Embodiment 1. A method for controlling the accuracy of large-curvature and shallow-drawing formed parts is provided, including:

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

[0170] S2. 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 during the unloading process of the formed part;

[0171] S3. Based on the strain energy and residual stress at the n nodes, cluster the n nodes to obtain m node regions; use the eagle-eye algorithm to solve the reverse bending deformation amount of the m node regions;

[0172] S4. Preset a tolerance range, and add a stiffness enhancement structure to the node regions where the reverse bending deformation amount exceeds the tolerance range to increase the stiffness of the corresponding node regions.

[0173] Embodiment 3

[0174] This embodiment publicly provides an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the computer program, it implements the operation mode of the above-provided method for controlling the accuracy of large-curvature and shallow-drawing formed parts.

[0175] Since the electronic device introduced in this embodiment is the electronic device used to implement a method for controlling the accuracy of a large-curvature and shallow-drawing formed part in the embodiments of the present application, based on the method for controlling the accuracy of a large-curvature and shallow-drawing formed part introduced in the embodiments of the present application, those skilled in the art can understand the specific implementation manners and various variations of the electronic device in this embodiment. Therefore, the specific implementation of how this electronic device implements the method in the embodiments of the present application will not be described in detail here. As long as those skilled in the art implement the electronic device used in the method for controlling the accuracy of a large-curvature and shallow-drawing formed part in the embodiments of the present application, it falls within the scope of protection of the present application.

[0176] The above formulas are all dimensionless and take their numerical values for calculation. The formulas are obtained by collecting a large amount of data for software simulation to obtain a formula that is closest to the actual situation. The preset parameters and threshold selection in the formulas are set by those skilled in the art according to the actual situation.

[0177] The above are only the preferred implementation manners of the present invention. The protection scope of the present invention is not limited to the above embodiments. All technical solutions within the concept of the present invention belong to the protection scope of the present invention. It should be noted that for ordinary technical users in the technical field, several improvements and refinements made without departing from the principle of the present invention should also be regarded as within the protection scope of the present invention.

Claims

[[ID=--0]]1. A precision control system for a large curvature and shallow drawing formed part, characterized in that, Including: A model construction module, which is used to obtain the initial geometric parameters and material mechanical property parameters of the formed part; And establish a finite element analysis model of the formed part; An energy analysis module, which is 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 during the unloading process of the formed part; A deformation amount fitting module, which clusters the n nodes based on the strain energy and residual stress at the n nodes to obtain m node regions; uses the eagle eye algorithm to solve the reverse bending deformation amount of the m node regions; A detection and adjustment module, which is used to preset a tolerance range, and add a stiffness enhancement structure to the node region where the reverse bending deformation amount exceeds the tolerance range; Each module is connected in a wired and / or wireless manner; The solution method of the reverse bending deformation amount includes: For each node area, extract the coordinates of the nodes within the node area and use them as the positions of the eagles; initialize the parameters of the eagle-eye algorithm, where the parameters include the number of eagles , the maximum number of iterations, and the discovery rate; Encoding each eagle, and the encoding method adopts real number encoding or binary encoding, and the encoding length is equal to the number of nodes in the node region multiplied by 3; Define the competition function ; Calculate the value of the competition function for each eagle, denoted as the competition value; Sort the eagles in the eagle group according to the competition value, and update the position of each eagle based on the position of the eagle with the highest competition value in the current eagle group and its own current position; After each position update, recalculate the competition value of each eagle in the eagle group, and based on the preset territory occupancy rate , retain the top eagles, and eliminate the remaining eagles; For the eagles retained in each iteration, with a preset probability perform mutation to obtain the new positions of the eagles; At this time, the eagles form a new generation of eagle group, calculate the competition value of the eagles in the new generation population, and retain the eagle with the highest competition value as the seed eagle of the next generation of eagle group; 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 encoding corresponding to the eagle with the highest competition value as the optimal value, decode the optimal value, and obtain the optimal plane fitting coordinates of the nodes in the node region; take the difference between the coordinates of each node in the node region and the optimal plane fitting coordinates as the reverse bending deformation amount of the node; take the average value of the reverse bending deformation amounts of each node in the node region as the reverse bending deformation amount of the node region; The competition function has the following formula: ; wherein, , and are for finding balance weight parameters; is a deformation state function, is the strain energy of node within the node region, is the residual stress of node within the node region, is the curvature of node within the node region, is a plane deviation function, is the fitting plane for node within the node region; Deformation state function ; wherein is the deformation balance parameter; ; wherein, is a planar balance parameter; The formula for updating the position is: ; where, is the position of the eagle at the -th iteration, and is the position of the eagle at the -th iteration, is the circumference ratio, is the saccade radius factor, is the global factor; is the position of the eagle with the highest competition value in the current eagle group; is the position of the eagle with the highest competition value in the current eagle group; Sweeping radius factor ; wherein is a constant Global factor ; wherein , , , and are all randomly updated numbers within the range of [0, 1]; The formula for performing mutation is: ; where, is the scaling factor, is a mutation random number within the range of [0, 1], is the upper boundary of the coordinates of the node, is the lower boundary of the coordinates of the node; is the position of the mutated eagle, is the position of the eagle before mutation.

2. The precision control system for large curvature and shallow drawing formed parts according to claim 1, wherein The initial geometric parameters include length, width, thickness, and curvature radius; the material mechanical property parameters include elastic modulus, Poisson's ratio, yield strength, and plastic parameters.

3. The precision control system for large curvature and shallow drawing formed parts according to claim 2, wherein The establishment method of the finite element analysis model of the formed part includes: According to the initial geometric parameters of the formed part, establish a three-dimensional solid model of the formed part; define material properties on the three-dimensional solid model based on the material mechanical property parameters to obtain a preliminary finite element model; For the preliminary finite element model, extract its outer surface as the outer boundary surface, and use geometric Boolean operations to extract the inner boundary surface from the preliminary finite element model; based on the outer boundary surface and the inner boundary surface, extract the intermediate surface; Pre-define a parameter domain, project the extracted intermediate surface isogonally into the parameter domain, and perform mesh division in the parameter domain to generate a mesh topology relationship of rectangular elements; The method of performing mesh division includes: Uniformly insert nodes on the inside and boundary of the parameter domain and number the nodes. Each adjacent four nodes form a rectangular element. For each rectangular element, determine the numbers of the nodes at its four vertices, and determine the element topology relationship according to the numbers of the nodes; Define the physical domain. According to the mesh topology relationship of rectangular elements, shell elements are constructed in the physical domain. For each node in the parametric domain of a rectangular element, the coordinates in the corresponding physical domain are calculated through a surface equation, and the coordinates in the physical domain are the control point coordinates of the shell element. The numbering of the nodes in the parametric domain is used as the numbering of the nodes of the shell element in the physical domain. For each shell element , construct its shape functions ; ; among them, is the parameter direction is the number of the upper node, is the one-dimensional B-spline basis function in the parameter direction ; is the weight of the node numbered in the parameter direction ; is the parameter direction is the number of the upper node, is the one-dimensional B-spline basis function in the parameter direction ; is the weight of the node numbered in the parameter direction ; Obtain the strain-displacement matrix of the shell element and the stress-strain matrix ; Based on the strain-displacement matrix of the shell element and the stress-strain matrix calculate the stiffness matrix of the shell element ; Calculate the load matrix based on the shape function of the shell element ; Assemble the stiffness matrices and load matrices of all shell elements separately to obtain the global stiffness matrix and the global load matrix , and handle the continuity between shell elements to obtain the finite element analysis model of the formed part.

4. The precision control system for large curvature and shallow drawing formed parts according to claim 3, characterized in that, The method for extracting the intermediate surface includes: Define the signed distance factor , representing the shortest distance from an arbitrary point within the preliminary finite element model to the inner boundary surface and the outer boundary surface; where , and are the abscissa, ordinate, and vertical coordinate of the arbitrary point in the space coordinate system where the preliminary finite element model is located, respectively; Define the deviation of the intermediate surface ; where is the shortest distance value from any point to the inner boundary surface, is the shortest distance value from any point to the outer boundary surface, is the weighting coefficient of the inner boundary surface, is the weighting coefficient of the outer boundary surface; Discretize the signed distance factor into a three-dimensional voxel grid, which serves 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. For each voxel, check whether there is an intersection with its adjacent voxels If there is an intersection, calculate the coordinates of the intersection point using linear interpolation; smoothly connect the calculated coordinates of the intersection points to obtain the intermediate surface.

5. The precision control system for large curvature and shallow drawing formed parts according to claim 4, characterized in that The method for discretizing the finite element analysis model of the formed part into n nodes includes: Perform an initial mesh division on the finite element analysis model of the formed part to obtain an initial mesh. Define a finite element function, and the formula of the finite element function is: ; Obtain the initial solution by solving the finite element function on the initial mesh ; Among them, is the load vector; Based on the initial solution , for each shell element calculate the corresponding error estimator ; ; wherein, is the true solution, is the elasticity matrix, is the strain tensor of Preset relative error , for shell elements marked as needing encryption, the shell elements marked as needing encryption are encrypted and refined to obtain a new mesh; The method for encryption and refinement includes: Represent the entire computational domain with a hexahedral element as the root node. Perform equal division and refinement on the root node to generate 8 sub-hexahedral elements as the child nodes of the root node. Recursively perform equal division and refinement on each child node until the preset maximum refinement level is reached, that is, construct a complete hexahedral tree hierarchical structure, and each node represents a hexahedral element. For the shell elements marked for encryption , map them to the corresponding nodes in the hexahedron tree hierarchy, i.e., the nodes are marked. For the marked nodes, generate 8 new child nodes, representing the 8 refined sub-hexahedron elements. Delete all the child nodes of the original node and replace them with the 8 generated new child nodes. Recursively repeat this process until the preset maximum refinement level is reached; Preset error threshold and re-solve the finite element function on the new mesh to obtain a new solution If, based on the new solution For each shell element calculate that the corresponding error estimator is greater than or equal to the error threshold , then repeat the encryption refinement; until the calculated error estimator is less than the error threshold Stop the encryption refinement to obtain the final mesh, and the final mesh contains n nodes.

6. The precision control system for large curvature and shallow drawing formed parts according to claim 5, characterized in that The calculation method for the strain energy and residual stress at the n nodes includes: Solve the finite element function on the final mesh to obtain the final solution ; the final solution contains the nodal displacement vectors of all shell elements; based on the final solution , calculate the strain field and stress field ; Based on the strain field of each shell element and the stress field , calculate the strain energy of each shell element ; ; wherein, represents the volume integration of the shell element over its volume; is the transpose of the strain field ; is the non-linear constitutive matrix, is the plastic strain, is the temperature, is the strain rate; is the geometric mapping Jacobian determinant of the shell element; ; wherein is the elastic matrix; is the plastic matrix; is the viscoplastic matrix; Based on the plastic matrix and the elastic matrix, the residual stress of each shell element is calculated ; ; wherein, is the plastic strain vector, is the tensor product operator; For each shell element, calculate the average strain energy and average residual stress of each shell element based on the calculated strain energy and residual stress. Based on the average strain energy and average residual stress of each shell element, use an interpolation function to perform interpolation operations at the n nodes to obtain the strain energy and residual stress of each node.

7. The precision control system for large curvature and shallow drawing formed parts according to claim 6, characterized in that The method for clustering the n nodes includes: Construct a w-dimensional feature vector from the strain energy and residual stress of each node respectively. Define the feature vectors of the n nodes as the initial positions of n fireflies in the w-dimensional space. Define the number of clusters and the clustering objective function ; Define that the luminous intensity of a firefly is inversely proportional to the value of the objective function; calculate the attractiveness between the th firefly and the th firefly; ; wherein, is the Euclidean distance between the th firefly and the th firefly, is the preset maximum attractiveness, is the preset light absorption coefficient; is the adjustment parameter, is the th firefly and the th firefly between the direction angle, is the preset desired direction angle; Move each firefly according to the attractiveness of all other fireflies. The movement formula is as follows: ; where is the step factor, is a random number within the range of [0, 1], is the new position of the -th firefly, is the position of the -th firefly before movement, is the position of the -th firefly before movement; Calculate the value of the clustering objective function for the new position of the firefly. If the value of the clustering objective function for the new position is higher than that of the original position, update the firefly to the new position, and update the luminous intensity of each firefly according to the new position of the firefly. Iterate until the preset number of iterations is satisfied to obtain the final positions of the fireflies. According to the final positions of the fireflies, group adjacent fireflies into the same cluster, and each cluster corresponds to a node region, thereby obtaining m node regions.

8. A method for controlling the accuracy of a large-curvature and shallow-drawing formed part, which is implemented based on the large-curvature and shallow-drawing formed part accuracy control system described in any one of claims 1 to 7, characterized in that It includes: S1. Obtain the initial geometric parameters and material mechanical property parameters of the formed part; And establish a finite element analysis model of the formed part; S2. 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 during the unloading process of the formed part; S3. Based on the strain energy and residual stress at the n nodes, cluster the n nodes to obtain m node regions; use the eagle-eye algorithm to solve the reverse bending deformation of the m node regions; S4. Preset a tolerance range, and add a stiffness enhancement structure to the node regions where the reverse bending deformation exceeds the tolerance range.

Images