Optimization Method for Converting Movable Deformable Components Oriented to Pixel-based Topology Optimization Results

Abstract

This paper proposes a method for converting movable and deformable components based on pixel-based topology optimization results. This method can transform structural optimization results represented by a two-dimensional finite element mesh into a component layout represented by an explicit function within a given error, establishing a geometric bridge between mesh-based and component-based structural optimization methods. Furthermore, a shape optimization method for minimizing maximum stress is proposed, which can significantly optimize the shapes generated by the conversion method, alleviate stress concentration in the structure, and effectively improve mechanical properties.

Description

Technical Field

[0001] The present invention belongs to the field of computer-aided design, and in particular relates to a shape approximation method for converting pixel-based topology optimization results into deformable movable components, and a corresponding shape optimization algorithm for the deformable movable components. Background Art

[0002] In the past decade, two important structural optimization methods, topology optimization and shape optimization, have garnered widespread attention in industry. Their specific objectives are to optimize material distribution within a given design area based on given load conditions, constraints, and performance metrics. Topology optimization can modify the overall structural defect, while shape optimization only adjusts the overall structural boundary shape. These methods, which combine geometric optimization techniques with finite element and equivalent geometric analysis calculations, effectively guide geometric modelers in improving the actual physical performance of designed structures and reduce the time spent on physical testing, feedback, and redesign during the product design process.

[0003] The research on topology optimization and shape optimization methods is a challenging task with great practical application value. Since the introduction of penalized solid isotropic material topology optimization methods, this field has seen explosive growth. Based on the geometric tools employed, the more mature methods currently include penalized solid isotropic material optimization and progressive structural optimization based on finite element mesh descriptions, and level set and image field methods based on implicit boundary descriptions. Both of these methods produce outputs with implicit descriptions of geometric boundaries and topology, hindering their direct application in secondary modeling and decomposition. To address this issue, Fuchs and Shemesh et al. proposed a shape optimization method using spline curves / surfaces. This method produces an accurate spline representation of the output geometric boundaries, allowing direct input into CAD software for processing. However, the topology description remains implicit. To address this issue, Guo Xu et al. proposed a topology optimization method based on movable deformable components in 2014. This method combines the advantages of both methods and innovatively introduces the movable deformable component as a geometric tool, significantly facilitating the control of topology and geometry during the optimization process. In the framework of the method based on deformable mobile components, the energy function converges quickly, but the optimality of the solution is difficult to theoretically guarantee due to the non-convex component function construction, and the accuracy is slightly lower than that of the method based on finite element meshes. Therefore, the present invention proposes a componentized shape approximation method for pixel-based topology optimization results. The topology optimization pixel-based results generated by the finite element mesh method can be used as input to quickly generate a mobile componentized approximation under the user-given error. In addition, based on the above-mentioned conversion method, the present invention proposes an optimization method for minimizing the maximum von Mises stress of the component, realizing mechanical optimization of the shape approximation results, so that the output model achieves better results in both shape error and mechanical properties. Summary of the Invention

[0004] The present invention addresses the shortcomings of existing technologies by proposing a method for optimizing the conversion of movable and deformable components based on pixel-based topology optimization results. This method analyzes local symmetry based on the mesh's inherent skeletal structure and establishes a transformation from the pixelated quadrilateral mesh generated by the topology optimization method to a movable and deformable component layout. Shape optimization is then performed within this component layout, aiming to minimize the maximum von Mises stress. This ensures consistent topological connectivity between components during the shape optimization process.

[0005] The specific solution of the present invention is: a method for optimizing the conversion of movable deformable components based on pixel-based topology optimization results, the method comprising the following steps:

[0006] Step 1: Input the 2D regular quadrilateral mesh generated by the structural optimization method, define the mesh unit set contained in the two, and retain the area as X in , the invalid area is X bdy, retain the cell set on the region boundary as X, and calculate X by the 8SSEDT method in and X bdy The signed distance field on the region, in addition to the maximum allowed approximation error D max ;

[0007] Step 2: A flux-based skeletonization algorithm is used to accurately extract the skeleton structure of the retained area. The skeleton line is composed of grid cells with a single unit width. A continuous parameterized relationship between the boundary line and the skeleton line is established based on the parameter relationship between the boundary points and the skeleton points. Then, the involved quantities are calculated by difference in the discrete case.

[0008] Step 3: Based on the parameterized representation in step 2, the distribution of bilateral symmetry lines can be further obtained. Then, based on the local symmetry analysis, a sparse optimization formula and a maximum approximation error constraint are established. The solution is optimized using the augmented Lagrangian optimizer, and the segmentation points on the bilateral symmetry lines are extracted using the threshold to achieve decomposition of the bilateral symmetry lines.

[0009] Step 4: Extract the component symmetry axis corresponding to each segment on the bilateral symmetry through principal component analysis, and calculate its position parameters;

[0010] Step 5: Under the premise of fixing the component symmetry axis position parameters obtained in step 4, the component symmetry axis position and distribution are calculated by using the variational shape approximation method;

[0011] Step 6: Using component parameters and mechanical boundary conditions as input, establish a model that minimizes the maximum von Mises stress. Use the augmented Lagrangian optimizer combined with finite element analysis to iteratively optimize the shape. When the convergence condition is reached, obtain a component layout that eliminates stress concentration.

[0012] As a preference, in step 1, the 8SSEDT distance field calculation method is used to calculate the signed distance on the two-dimensional area. In step 2, a flux-based skeletonization algorithm is used to generate a skeleton, and the structure is used as a reference to extract the component symmetry axis distribution. In step 3, the theory of bilateral symmetry lines and bilateral symmetry is established, and based on this, a sparse optimization model for minimizing bilateral symmetry line segmentation is designed. Then, the segmentation points are obtained by solving the augmented Lagrangian optimizer. In step 3, the theory of bilateral symmetry lines and bilateral symmetry is established, and the sparse optimization model for minimizing bilateral symmetry line segmentation is designed.

[0013] In step 4, principal component analysis is used to estimate the position of the component's symmetry axis, and a variational shape approximation optimization model is given in step 5.

[0014] In step 6, a shape optimization framework based on component description is established, and a specific optimization model is established to minimize the maximum von Mises stress problem. Finally, the optimal parameters of the component are obtained through an iterative optimization process. More specifically:

[0015] Step 1: First, input the 2D regular quadrilateral mesh generated by the topology optimization method. The unit side lengths in the mesh are equal and there are no singular points. It is a rectangular range, which contains two types of areas: the reserved area and the invalid area. The unit is accompanied by a label to define the set of mesh units contained in the two. The reserved area is X in , the invalid area is X bdy , the set of cells on the boundary of the reserved area is X. The 8SSEDT (8-point signed sequential euclidean distance transform) method is used to calculate the signed distance field from the cells in the grid area to the boundary of the reserved area. The detailed steps of this method are given in Section 2 of the "Note on "FastRaster Scan Distance Propagation on the Discrete Rectangular Lattice". This distance field records the distance from the center point P of each cell to the center point of the nearest cell on the boundary of the reserved area, which is recorded as the minimum offset distance D x (P), also records the nearest cell center point P' on the boundary of the reserved area. The unit vector from P to P' is called the minimum offset distance vector. The above distance measurement is positive when the cell is inside the reserved area and negative when it is outside. In addition, the maximum allowable approximation error D max , this error value is used to limit the maximum deviation distance between the converted component boundary and the effective area grid boundary.

[0016] Step 2: Based on the distance field, the flux-based skeletonization algorithm is used to calculate the skeleton structure of the retained area. The specific calculation steps are given in Section 4 of Flux invariants for shape. The skeleton structure uses the grid unit as the basic unit, and the skeleton points are as follows: Figure 1 As shown, there are four types of classifications: regular skeleton points, intersection points, pseudo intersection points and end points. The entire skeleton structure can be decomposed into m branch lines using intersection points and pseudo intersection points as segmentation points. The skeleton obtained by the above operation can be expressed as follows:

[0017] SK=∪S i (t)(t∈[0,1],i=1,2,…,m)

[0018] Where SK is the skeleton, S i(t) is the i-th parameterizable branch line. The specific parameterization is established using the chord length parameterization method. The specific process is calculated according to the paper "Curves with chord length parameterization". The skeleton point P on the above branch line (P is the center of the grid unit) and the closest point Q on both sides of the boundary are i (i=1,2) the following relationship exists:

[0019] Q i =P+D X (P)R((-1) (i+1) α)n(P)(i=1,2)

[0020] Where α is the angle between the skeleton line tangent and the minimum offset distance vector of the boundary of P, which is called the projection angle, R((-1) (i+1) α) is the rotation matrix for clockwise and counterclockwise rotations of angle α when 1 and 2 are respectively, and n(P) is the unit tangent vector along the branch line at point P. The following relationship can also be established on the branch line S(t) with continuous parameterization:

[0021]

[0022] Among them C i (t)(i=1,2) is the parameterized representation of the boundary curves on both sides, r(t) is the minimum offset distance of the skeleton point at the parameter t position on the branch line, S′(t) is the unit tangent vector at the parameter t position, S′ ⊥ (t) is the unit vector obtained by rotating the tangent vector at the parameter t position counterclockwise by π / 2. The above mapping relationship can be accurately described in the continuous curve boundary. In the discrete case with pixel grid as the basic unit, there will be a loss of accuracy. If the point S(t) falls in the unit i, then the center point P of the unit i is used as the i The approximate value S′(P i ),α(P i ), D X (P i ) to replace S′(t), α(t), and r(t). Furthermore, since cosα(t) = |r(t)′|, where α(t) is the projection angle at the position of parameter t, the above equations support the approximate calculation of several members as follows:

[0023]

[0024]

[0025] Where H is P iThe number of points on the branch line that are included in the difference calculation in the forward and backward directions is mainly used to reduce the numerical error under discrete sampling. In this paper, H is set to 2, which can effectively reduce the impact of numerical errors.

[0026] Step 3: The skeleton structure is obtained in step 2. i (t) can continue to establish the bilateral symmetry line M i (t) are as follows:

[0027] M i (t) = S i (t)-r i (t)S′ i (t)r′ i (t)

[0028] C1(t) and C2(t) in M i (t) about S′ i (t) Directional symmetry. If there exists t∈(a,b), the symmetry directions of the points on the bilateral symmetry line are consistent, that is, S″ i (t) = 0, then it is said that there is consistent bilateral symmetry on the interval. When the interval satisfies consistent bilateral symmetry, the branch line can be used as the symmetry axis of a component. Therefore, the target can be established on each branch line. S′ i (t) is the model of the optimized variable. By optimizing this formula, the shape corresponding to each branch line can be constructed with the least components. In the discrete case, the optimized variable S in the above formula is i '(t) can be transformed into the skeleton point P j i (i=1,…,m,j=1,…,k i )(k i is the number of skeleton points on the i-th branch line) divided by the original S′ i (0) is the counterclockwise rotation angle of the reference vector, β * (P j i ). Based on the above description, a LASSO (Least absolute shrinkage and selection operator) regression optimization model can be established in the discrete case:

[0029]

[0030]

[0031]

[0032]

[0033]

[0034] In the above formula, β * (P) and β(P) are the optimized rotation angle and initial rotation angle at point P, respectively. i To improve the sparsity coefficient of the solution on the smooth skeleton branch, the purpose of the first set of constraints is to control the approximation error, and the purpose of the second and third sets of equality constraints is to fix the rotation angles of the first and last points on the branch line to prevent the solution from degenerating into the same constant. The specific optimization process of the above optimization formula can be described as follows:

[0035] (1) Calculate all skeleton points P on the branch line j i The rotation angle β(P j i ), the initial value of the optimization variable is set to β * (P j i )=β(P j i ), set equality constraints and inequality constraints;

[0036] (2) Use the augmented Lagrangian method to solve the above-mentioned weighted LASSO regression problem according to the initial values and constraints in (1). The specific steps of the augmented Lagrangian method are described in "A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds";

[0037] (3) Output the optimization results. If the absolute value of the rotation angle difference between nodes is Then point P j i Extract as segmentation points;

[0038] (4) On each bilateral symmetry line M i The segmentation point is used as the boundary, and it can be decomposed into several segments. For the bilateral symmetry line M i For example, if there are c segment points on it, it can be decomposed into c+1 segments, where the set of bilateral symmetric points corresponding to each segment is C ij (j=1,…,c+1), where is the coordinate vector of the kth point in the set, C ij There are n ij bilateral symmetry points.

[0039] Step 4: Based on the bilateral symmetric point set obtained in step 3, the component symmetry axis information is robustly extracted. First, the definition and schematic diagram of the movable deformable component are given, as shown in Figure 2 The figure shows a schematic diagram of a movable deformable component. For a two-dimensional direct axis movable deformable component i, the implicit function definition is as follows:

[0040]

[0041]

[0042] Among them L i is the half length of the component, and the center point of the component is (x oi ,y oi ), (x′, y′) is the counterclockwise rotation of θ around the center point i The coordinates obtained after the degree of measurement, P is the power value of the control boundary accuracy, which is set to 20 in the present invention. i (x′,T i ) is the boundary curve located on both sides of the symmetry axis, which can be a Lagrange interpolation curve or a spline curve, which takes the local coordinates and the half-width parameter set as parameters, where is the half-width parameter set. When the boundary curve is represented by a Lagrange interpolation curve, it can be listed as follows:

[0043]

[0044] When using B-spline curve, it can be expressed as follows:

[0045]

[0046] Among them B j,p (ζ) is the j-th p-order B-spline basis function, and its analytical expression is defined as follows:

[0047]

[0048]

[0049]

[0050] For component i, there are the following symmetry axes:

[0051]

[0052] Component i is symmetrical about the above-mentioned symmetry axis. When the position parameters of the symmetry axis are solved, the preliminary layout of the component is determined. The bilateral symmetry point set C obtained in step 3 is ij Convert to matrix form:

[0053]

[0054] For the matrix C ij Subtracting the mean center coordinate vector from each column vector in yields the following matrix:

[0055]

[0056] Among them are

[0057]

[0058] From this, the covariance matrix can be calculated as follows:

[0059]

[0060] The Jacobi iteration method is used to solve the eigenvectors and corresponding eigenvalues of the above covariance matrix. The largest eigenvalue is The corresponding eigenvector is The minimum eigenvalue is The corresponding eigenvector is Through the above principal component analysis, each set of bilateral symmetric points can be established into a component symmetry axis, where if the bilateral symmetric point set C ij The corresponding component symmetry axis is D d , then its specific parameters can be calculated as follows:

[0061]

[0062]

[0063]

[0064]

[0065] where v0 = [1,0] T , from which the symmetry axis parameters of the corresponding component d can be calculated;

[0066] Step 5: Sampling is established based on the input grid. The sampling points are divided into two types: reserved area sampling points and invalid area sampling points. The reserved area sampling points are established at the center points of the 4-layer grid cells inward from the boundary of the reserved area and the 3-layer grid cells near the intersection. The invalid area sampling points are established at the center points of the 4-layer grid cells outward from the boundary of the reserved area. The set of reserved area sampling points obtained is P in , the invalid area sampling point set is P bdy Based on the above sampling point set, the following discrete shape approximation energy function can be established, which is the error value (in units) between the component and the retained mesh area during the approximation process:

[0067]

[0068]

[0069]

[0070]

[0071] Where Ω0(P) is the piecewise characteristic function, which is 1 when the point is inside the reserved region, 0 when it is on the boundary, and -1 when it is outside the reserved region. Φ(P) is the global implicit function representation that uses an approximate max function to aggregate the components, and H(x) is the approximate Heaviside function as a filter, with a threshold of ε = 10 -4 The above shape approximation energy function is optimized using the Subplex method. The specific steps of Subplex are given in "Functional Stability Analysis of Numerical Algorithms". The entire optimization process can be described as follows:

[0072] Initialize all component half-width parameters as follows:

[0073]

[0074] (1) Using the Subplex method to iteratively optimize E Appr Minimization problem of

[0075] (2) Output the component width parameter. The other parameters of the component have been obtained in step 4.

[0076] The shape optimization method for the movable deformable component representation can be described as follows:

[0077] Step 1: Use the component parameters obtained from the above shape approximation method as input and expand the degrees of freedom. Convert the original straight-axis symmetric component into a straight-axis asymmetric component and set the endpoints of the symmetry axis on both sides to be variable in the parameterized direction. The definition of the straight-axis asymmetric component is as follows:

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084]

[0085]

[0086] In the above definition, the parameter θ of the input component i is i , x oi ,y oi , L i Just bring in the original value. For the width parameter set on the component boundary, it can be initialized to T l i =T i ,

[0087] Step 2: Set the mechanical boundary conditions on the background grid, i.e., external force vectors, displacement constraints, etc. For the problem of minimizing the maximum stress with volume constraints, the following optimization formula can be established:

[0088] Minimizeσ max (Φ)

[0089] st

[0090] K(Φ)U=F

[0091]

[0092] Where Φ is the global implicit function representation of the asymmetric component, f v To preserve the maximum area fraction of the material, V0 is the total area of the design domain, that is, the total area of the original background grid, V(Φ) is the total area of the region satisfying Φ(x, y)>0, K(Φ) is the global stiffness matrix under the current component distribution, U is the displacement vector, and F is the force vector. The maximum von Mises stress σ max The P-norm global stress assessment strategy can be approximated as follows:

[0093]

[0094]

[0095]

[0096] σ j =DB e U e

[0097] The unit geometry matrix B in the above column e , D is the elastic coefficient matrix, U e The components of the displacement vector U at the nodes of element e obtained by finite element calculation, the initial element stiffness matrix Where t is the unit thickness, A eThe initial unit stiffness matrix combined with the implicit function projection can be used to calculate the unit stiffness matrix K in the current state. e (Φ). The mapping of implicit functions to grid space is specifically in the following form:

[0098]

[0099]

[0100]

[0101] Where V e is the volume fraction of unit e, is the initial element stiffness matrix of element e, K e (Φ) is the element stiffness matrix obtained after component projection. After the above calculation, the stiffness matrix K(Φ) can be obtained by the stiffness matrix K of all elements. e (Φ). The above optimization model is optimized and solved using the augmented Lagrangian method. The optimization process can be described as follows:

[0102] (1) Input n parameters of the direct axis symmetric component and convert them into the direct axis asymmetric component according to the conversion method given in step 1. Given the displacement and force boundary conditions, if the number of grid points is N, the displacement boundary condition is recorded as the vector

[0103] The force boundary conditions are recorded as vectors Iteration number i = 1;

[0104] (2) Calculate the mapping value V from the current component implicit function Φ to the grid space e , and then calculate the area mapping value and unit stiffness matrix mapping value V(Φ) and K e (Φ);

[0105] (3) Assemble the unit stiffness matrix in (2) into a global stiffness matrix, calculate the finite element analysis, solve K(Φ)U=F, and further calculate the energy function value σ max ;

[0106] (4) Using the augmented Lagrangian method to perform the i-th optimization;

[0107] (5) If the stress function value of the i-th optimization output is and the stress function value outputted at the i-1th time satisfy Then stop the iteration and execute step (6), otherwise i increases by 1 and returns to step (2);

[0108] (6) Output the parameters of n direct-axis asymmetric components.

[0109] The beneficial effects of the present invention are as follows: the present invention combines skeleton extraction, sparse optimization methods and principal component analysis to establish a conversion method from pixel-based topology optimization results to movable deformable components, thereby realizing explicit topological construction of discrete topology optimization results. The output result of this method inherits the component result of the skeleton structure. Compared with other implicit conversion methods, this structure has greater application value in engineering design. The topology can be globally or locally maintained on the basis of this structure to adjust the shape and perform redesign calculations. In addition, the present invention also proposes a shape optimization method for minimizing the maximum von Mises stress for the above-mentioned component structure. Combined with the above-mentioned conversion method, the component parameters output by the conversion method are used as input, and the optimal shape design for stress can be achieved under given boundary conditions, which has the advantage of maintaining the skeleton topological connection during the optimization process. BRIEF DESCRIPTION OF THE DRAWINGS

[0110] Figure 1 This is a schematic diagram of skeleton point classification;

[0111] Figure 2 Schematic diagram of the movable deformable component;

[0112] Figure 3 Bridge topology optimization mesh;

[0113] Figure 4 Signed distance field to the boundaries of the preserved regions in the bridge topology optimization mesh;

[0114] Figure 5 Flux field in the mesh-preserved area of bridge topology optimization;

[0115] Figure 6 The skeleton branch lines of the bridge topology optimization mesh retaining area;

[0116] Figure 7 The bridge topology optimization mesh retains the bilateral symmetry lines of the region;

[0117] Figure 8 The skeleton segmentation points extracted from the above skeleton;

[0118] Figure 9 right Figure 1 The component symmetry axis extracted from the example;

[0119] Figure 10 right Figure 1 The component layout extracted from the example;

[0120] Figure 11 right Figure 1 The boundary condition settings used in the example shape optimization;

[0121] Figure 12 exist Figure 11 Under the boundary condition setting Figure 1 The example computes the minimum maximum von Mises stress,

[0122] Optimized component layout;

[0123] Figure 13 exist Figure 11 Under the boundary condition setting Figure 10 Calculate von Mises stresses in the structure composed of component layout;

[0124] Figure 14 exist Figure 11 Calculate von Mises stress on the structure formed by shape optimization results under boundary condition settings; DETAILED DESCRIPTION

[0125] The present invention will be further described below with reference to the accompanying drawings.

[0126] Optimization method for transforming movable deformable components based on pixel-based topology optimization results: Step 1: First, input the 2D regular quadrilateral mesh generated by the topology optimization method, where the reserved area is X in , is a white unit, and the invalid area is X bdy , are black cells, and the cell set on the boundary of the reserved region is X, such as Figure 3 The following is an example of an input bridge mesh. The 8SSEDT method is used to calculate the signed distance field from the cells in the mesh area to the boundary X of the reserved area. The distance field is as follows: Figure 4 In addition, the maximum allowable approximation error D must be entered. max ;

[0127] Step 2: Based on the distance field, a flux-based skeletonization algorithm is used to calculate the skeleton structure of the retained area. The skeleton structure uses grid cells as the basic unit. Figure 5 The flux field results are shown in Figure 6 The detailed skeleton branching results are shown in Figure 1, which has four categories: regular skeleton points, intersection points, pseudo-intersection points, and end points. Using intersection points and pseudo-intersection points as segmentation points, the entire skeleton structure can be decomposed into 34 branch lines, that is, m = 38. The skeleton obtained through the above operation can be expressed as follows:

[0128] SK=∪S i (t)(t∈[0,1],i=1,2,…,m)

[0129] Where SK is the skeleton, S i (t) is the i-th parameterizable branch line, and the specific parameterization is established using the chord length parameterization method. The skeleton point P on the above branch line (P is the center of the grid unit) and the closest point Q on both sides of the boundary i (i=1,2) the following relationship exists:

[0130] Q i =P+D X (P)R((-1) (i+1) α)n(P)(i=1,2)

[0131] Where α is the angle between the skeleton line tangent and the minimum offset distance vector of the boundary of P, which is called the projection angle, R((-1) (i+1) α) is the rotation matrix for clockwise and counterclockwise rotations of angle α when 1 and 2 are respectively, and n(P) is the unit tangent vector along the branch line at point P. The following relationship can also be established on the branch line S(t) with continuous parameterization:

[0132]

[0133] Among them C i (t)(i=1,2) is the parameterized representation of the boundary curves on both sides, r(t) is the minimum offset distance of the skeleton point at the parameter t position on the branch line, S′(t) is the unit tangent vector at the parameter t position, S′ ⊥ (t) is the unit vector obtained by rotating the tangent vector at the parameter t position counterclockwise by π / 2. The above mapping relationship can be accurately described in the continuous curve boundary. In the discrete case with pixel grid as the basic unit, there will be a loss of accuracy. If the point S(t) falls in the unit i, then the center point P of the unit i is used as the i The approximate value S′(P i ),α(P i ), D X (P i ) to replace S′(t), α(t), and r(t). Furthermore, since cosα(t) = |r(t)′|, where α(t) is the projection angle at the position of parameter t, the above equations support the approximate calculation of several members as follows:

[0134]

[0135]

[0136] Where H is P i The number of points on the branch line that are included in the differential calculation in the forward and backward directions is important for reducing the numerical error under discrete sampling. In this paper, H is taken as 2;

[0137] Step 3: The skeleton structure is obtained in step 2. i (t) can continue to establish the bilateral symmetry line M i (t) are as follows:

[0138] M i (t) = S i (t)-r i (t)S′i (t)r′ i (t)

[0139] C1(t) and C2(t) in M i (t) about S′ i (t) Direction symmetry, such as Figure 7 The results of the bilateral symmetry line calculation are shown. Calculate the LASSO regression optimization model in the discrete case:

[0140]

[0141]

[0142]

[0143]

[0144]

[0145] In the above formula, β * (P) and β(P) are the optimized rotation angle and initial rotation angle at point P, respectively. i To improve the sparsity coefficient of the solution on the smooth skeleton branch, the purpose of the first set of constraints is to control the approximation error, and the purpose of the second and third sets of equality constraints is to fix the rotation angles of the first and last points on the branch line to prevent the solution from degenerating into the same constant. The specific optimization process of the above optimization formula can be described as follows:

[0146] (5) Calculate all skeleton points P on the branch line j i The rotation angle β(P j i ), the initial value of the optimization variable is set to β * (P j i )=β(P j i ), set equality constraints and inequality constraints;

[0147] (6) Using the augmented Lagrangian method, solve the above weighted LASSO regression problem according to the initial values and constraints in (1);

[0148] (7) Output the optimization results. If the absolute value of the rotation angle difference between nodes is Then point P j i Extract as segmentation points;

[0149] (8) On each bilateral symmetry line M i The segmentation point is used as the boundary, and it can be decomposed into several segments. For the bilateral symmetry line M iFor example, if there are c segment points on it, it can be decomposed into c+1 segments, where the set of bilateral symmetric points corresponding to each segment is C ij (j=1,…,c+1), where is the coordinate vector of the kth point in the set, C ij There are n ij bilateral symmetry points.

[0150] Through the above calculations, we can extract Figure 8 The segmentation point distribution shown in the figure, where the white points on the light white skeleton line are segmentation points;

[0151] Step 4: Based on the bilateral symmetric point set obtained in step 3, the component symmetry axis information is robustly extracted. First, the definition and schematic diagram of the movable deformable component are given, as shown in Figure 2 The figure shows a schematic diagram of a movable deformable component. For a two-dimensional direct axis movable deformable component i, the implicit function definition is as follows:

[0152]

[0153]

[0154] Among them L i is the half length of the component, and the center point of the component is (x oi ,y oi ), (x′, y′) is the counterclockwise rotation of θ around the center point i The coordinates obtained after the degree of measurement, P is the power value of the control boundary accuracy, which is set to 20 in the present invention. i (x′,T i ) is the boundary curve located on both sides of the symmetry axis, which can be a Lagrange interpolation curve or a spline curve, which takes the local coordinates and the half-width parameter set as parameters, where is the half-width parameter set. When the boundary curve is represented by a Lagrange interpolation curve, it can be listed as follows:

[0155]

[0156] When using B-spline curve, it can be expressed as follows:

[0157]

[0158] Among them B j,p (ζ) is the j-th p-order B-spline basis function, and its analytical expression is defined as follows:

[0159]

[0160]

[0161]

[0162] For component i, there are the following symmetry axes:

[0163]

[0164] Component i is symmetrical about the above-mentioned symmetry axis. When the position parameters of the symmetry axis are solved, the preliminary layout of the component is determined. The bilateral symmetry point set C obtained in step 3 is ij Convert to matrix form:

[0165]

[0166] For the matrix C ij Subtracting the mean center coordinate vector from each column vector in yields the following matrix:

[0167]

[0168] Among them are

[0169]

[0170] From this, the covariance matrix can be calculated as follows:

[0171]

[0172] The Jacobi iteration method is used to solve the eigenvectors and corresponding eigenvalues of the above covariance matrix. The largest eigenvalue is The corresponding eigenvector is The minimum eigenvalue is The corresponding eigenvector is Through the above principal component analysis, each set of bilateral symmetric points can be established into a component symmetry axis, where if the bilateral symmetric point set C ij The corresponding component symmetry axis is D d , then its specific parameters can be calculated as follows:

[0173]

[0174]

[0175]

[0176]

[0177]

[0178] where v0 = [1,0] T , from which the symmetry axis parameters of the corresponding component d can be calculated. Figure 9 Shown is the extracted distribution of component symmetry axes;

[0179] Step 5: Sampling is established based on the input grid. The sampling points are divided into two types: reserved area sampling points and invalid area sampling points. The reserved area sampling points are established at the center points of the 4-layer grid cells inward from the boundary of the reserved area and the 3-layer grid cells near the intersection. The invalid area sampling points are established at the center points of the 4-layer grid cells outward from the boundary of the reserved area. The set of reserved area sampling points obtained is P in , the invalid area sampling point set is P dby Based on the above sampling point set, the following discrete shape approximation energy function can be established, which is the error value (in units) between the component and the retained mesh area during the approximation process:

[0180]

[0181]

[0182]

[0183]

[0184] Where Ω0(P) is the piecewise characteristic function, which is 1 when the point is inside the reserved region, 0 when it is on the boundary, and -1 when it is outside the reserved region. Φ(P) is the global implicit function representation that uses an approximate max function to aggregate the components, and H(x) is the approximate Heaviside function as a filter, with a threshold of ε = 10 -4 The above shape approximation energy function is optimized using the Subplex method. The entire optimization process is as follows:

[0185] Initialize all component half-width parameters as follows:

[0186]

[0187] (3) Using the Subplex method to iteratively optimize E Appr Minimization problem of

[0188] (4) Output the component width parameter. The other parameters of the component have been obtained in step 4.

[0189] The shape approximation result obtained according to the above optimization process is shown in Figure 10.

[0190] The shape optimization method for movable deformable components can be described as follows:

[0191] Step 1: Use the component parameters obtained from the above shape approximation method as input and expand the degrees of freedom. Convert the original straight-axis symmetric component into a straight-axis asymmetric component and set the endpoints of the symmetry axis on both sides to be variable in the parameterized direction. The definition of the straight-axis asymmetric component is as follows:

[0192]

[0193]

[0194]

[0195]

[0196]

[0197]

[0198]

[0199]

[0200] In the above definition, the parameter θ of the input component i is i , x oi ,y oi , L i Just bring in the original value. For the width parameter set on the component boundary, it can be initialized to

[0201] Step 2: Set the mechanical boundary conditions on the background grid, i.e., external force vectors, displacement constraints, etc. For the problem of minimizing the maximum stress with volume constraints, the following optimization formula can be established:

[0202] Minimizeσ max (Φ)

[0203] st

[0204] K(Φ)U=F

[0205]

[0206] Where Φ is the global implicit function representation of the asymmetric component, f v To preserve the maximum area fraction of the material, V0 is the total area of the design domain, that is, the total area of the original background grid, V(Φ) is the total area of the region satisfying Φ(x,y)>0, K(Φ) is the global stiffness matrix under the current component distribution, U is the displacement vector, and F is the force vector. The maximum von Mises stress σ max The P-norm global stress assessment strategy can be approximated as follows:

[0207]

[0208]

[0209]

[0210] σ j =DB e U e

[0211] The unit geometry matrix B in the above column e , D is the elastic coefficient matrix, U e The components of the displacement vector U at the nodes of element e obtained by finite element calculation, the initial element stiffness matrix Where t is the unit thickness, A e The initial unit stiffness matrix combined with the implicit function projection can be used to calculate the unit stiffness matrix K in the current state. e (Φ). The mapping of implicit functions to grid space is specifically in the following form:

[0212]

[0213]

[0214]

[0215] After the above calculation, the stiffness matrix K(Φ) can be obtained by the stiffness matrix K of all elements. e (Φ). The above optimization model is optimized and solved using the augmented Lagrangian method. The optimization process is as follows:

[0216] (1) Input n parameters of the direct axis symmetric component and convert them into the direct axis asymmetric component according to the conversion method given in step 1. Given the displacement and force boundary conditions, if the number of grid points is N, the displacement boundary condition is recorded as the vector The force boundary conditions are recorded as vectors Iteration number i = 1;

[0217] (2) Calculate the mapping value V from the current component implicit function Φ to the grid space e , and then calculate the area mapping value and unit stiffness matrix mapping value V(Φ) and K e (Φ);

[0218] (3) Assemble the unit stiffness matrix in (2) into a global stiffness matrix, calculate the finite element analysis, solve K(Φ)U=F, and further calculate the energy function value σ max ;

[0219] (4) Using the augmented Lagrangian method to perform the i-th optimization;

[0220] (5) If the stress function value of the i-th optimization output is and the stress function value outputted at the i-1th time satisfy Then stop the iteration and execute step (6), otherwise i increases by 1 and returns to step (2);

[0221] (6) Output the parameters of n direct-axis asymmetric components.

[0222] When the volume fraction is set to 34%, the boundary conditions are as follows Figure 11 When shown, the optimization result is Figure 12 The results shown in Figure 13 、 14 The initial stress field and the optimized stress field are shown in Figure 3, respectively.

Claims

1. A method for optimizing the transformation of movable deformable components based on pixel-wise topology optimization results, characterized by: The method comprises the following steps: Step 1: Input the 2D regular quadrilateral mesh generated by the structural optimization method and define the mesh unit set contained in the reserved area and invalid area. The reserved area is X in , the invalid area is X bdy , retain the cell set on the region boundary as X, and calculate X by the 8SSEDT method in and X bdy The signed distance field on the region, in addition to the maximum allowed approximation error D max ; Step 2: A flux-based skeletonization algorithm is used to accurately extract the skeleton structure of the retained area. The skeleton line is composed of grid cells with a single unit width. A continuous parameterized relationship between the boundary line and the skeleton line is established based on the parameter relationship between the boundary points and the skeleton points. Then, in the discrete case, the involved quantities are calculated by difference method. Step 3: Based on the parameterized representation in step 2, the distribution of bilateral symmetry lines can be further obtained. Then, based on the local symmetry analysis, a sparse optimization formula and a maximum approximation error constraint are established. The solution is optimized using the augmented Lagrangian optimizer, and the segmentation points on the bilateral symmetry lines are extracted using the threshold to achieve decomposition of the bilateral symmetry lines. Step 4: Extract the component symmetry axis corresponding to each segment on the bilateral symmetry through principal component analysis, and calculate its position parameters; Step 5: Under the premise of fixing the component symmetry axis position parameters obtained in step 4, the component symmetry axis position and distribution are calculated by using the variational shape approximation method; Step 6: Using component parameters and mechanical boundary conditions as input, establish a model that minimizes the maximum von Mises stress. Use the augmented Lagrangian optimizer combined with finite element analysis to iteratively optimize the shape. When the convergence condition is reached, obtain a component layout that eliminates stress concentration.

2. The method for transforming and optimizing a movable deformable component based on pixel-wise topology optimization results according to claim 1, characterized in that: In step 1, the 8SSEDT distance field calculation method is used to calculate the signed distance in the two-dimensional area. In step 2, a flux-based skeletonization algorithm is used to generate a skeleton, and this structure is used as a reference to extract the component symmetry axis distribution. In step 3, the theory of bilateral symmetry lines and bilateral symmetry is established, and based on this, a sparse optimization model for minimizing bilateral symmetry line segmentation is designed. The segmentation points are then obtained by solving the augmented Lagrangian optimizer. In step 4, principal component analysis is used to estimate the component symmetry axis position. In step 5, a variational shape approximation optimization model is given.

3. The method for transforming and optimizing movable deformable components based on pixel-wise topology optimization results according to claim 1, characterized in that: In step 6, a shape optimization framework based on component description is established, and a specific optimization model is established for the problem of minimizing the maximum von Mises stress. Finally, the optimal parameters of the component are obtained through an iterative optimization process.

Images