Additive manufacturing path pretreatment method, system and equipment facing topological optimization design and medium
The method synchronizes unit density and orientation optimization with K-means clustering and NURBS surface fitting to address anisotropic issues in carbon fiber composites, enhancing additive manufacturing precision and structural reliability.
Patent Information
- Application Number
- CN202510342155.1
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-03-21
- Publication Date
- 2025-07-08
AI Technical Summary
Traditional topology optimization methods for carbon fiber reinforced composite materials fail to consider their anisotropic properties, leading to inconsistencies between optimization results and actual structural performance, especially in additive manufacturing, and result in jagged boundaries that complicate processing.
A method involving synchronous optimization of unit density and orientation, using K-means clustering and NURBS surface fitting to smooth boundaries and ensure consistency between optimization results and additive manufacturing processes.
Enhances the precision and reliability of additive manufacturing by ensuring structural performance consistency and smoothing jagged boundaries, thereby improving print quality and structural integrity.
Smart Images

Figure CN120278003A_ABST
Abstract
Description
Technical Field
[0001] The present invention belongs to the technical field of topology optimization and additive manufacturing, and particularly relates to a method, system, device and medium for preprocessing additive manufacturing paths for topology optimization design. Background Art
[0002] Carbon fiber reinforced composites have significant anisotropy, and their structures and properties vary significantly in different directions. Traditional topology optimization methods usually do not fully consider the anisotropic characteristics of materials, resulting in inconsistent optimization results with the structural properties required in practical applications, especially in the additive manufacturing process. In addition, due to the difference between the printing angle and the optimization result in the additive manufacturing process, the optimized design often cannot be fully realized. The serrated boundaries of the optimization results also pose difficulties for subsequent processing and manufacturing. Existing technologies usually focus on the matching problem between optimization design and manufacturing, but there is still a lack of an efficient solution for how to effectively ensure their consistency, especially when the optimization results involve different printing angles.
[0003] Therefore, it is urgent to solve the above problems. Summary of the Invention
[0004] Object of the Invention: The object of the present invention is to provide a method for preprocessing additive manufacturing paths for topology optimization design, which can effectively detect the boundary points of point clouds and achieve the consistency between the optimization result and the additive manufacturing process.
[0005] The second object of the present invention is to provide a system for preprocessing additive manufacturing paths for topology optimization design.
[0006] The third object of the present invention is to provide an electronic device.
[0007] The fourth object of the present invention is to provide a computer storage medium.
[0008] Technical Solution: To achieve the above objects, the present invention discloses a method for preprocessing additive manufacturing paths for topology optimization design, including the following steps:
[0009] S1: During the topology optimization process, considering the anisotropy of carbon fiber reinforced composites, synchronously optimize the element density and element angle, that is, synchronously optimize the fiber orientation. The output optimization result includes the best topology structure diagram and the density and angle of each discrete element to ensure that the structural performance of the optimization result meets the actual requirements;
[0010] S2: For the discrete angle result θ obtained during the optimization process, use the K-means clustering algorithm to process it, divide it into multiple regions, and obtain the dividing line of each region;
[0011] S3: For the zigzag boundary caused by mesh division in the optimal topological mesh diagram, use the NURBS surface fitting technology to smooth the boundary;
[0012] S4: Combine the smoothed boundary line and the regional demarcation obtained by K-means clustering to generate the boundary line of each region, and fill each region with straight lines at corresponding angles to generate the corresponding Gcoode code.
[0013] Optionally, the step S1 specifically includes the following steps:
[0014] S11: First, define the design domain size, force model, support conditions, and material properties, discretize the design domain into four-node rectangular elements, and obtain the initial element stiffness matrix and the global stiffness matrix according to the initial angle and initial density of the elements;
[0015] The compliance constitutive matrix of the carbon fiber reinforced composite material is S0 and the stiffness constitutive matrix is C0 as follows:
[0016]
[0017] C0 = S0 -1
[0018] Among them, E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, v 12 is the Poisson's ratio, G 12 is the in-plane shear modulus;
[0019] When the fiber local coordinate system forms an angle θ e with the global coordinate system, the stiffness constitutive matrix C(θ e ) of the e-th element is obtained by rotating the stiffness constitutive matrix C0 of the material:
[0020] C(θ e ) = T(θ e )C0T(θ e ) T
[0021] Among them, the rotation matrix is:
[0022]
[0023] The element stiffness matrix k e (θ e ) is:
[0024]
[0025] Among them, B is the strain-displacement matrix, and Ω e is the integration region of the element;
[0026] The assembled global stiffness matrix \(K(x,\theta)\) is as follows:
[0027]
[0028] where \(p\) is the penalty factor;
[0029] S12: Based on the Solid Isotropic Material with Penalization (SIMP) method, the design variables for topology optimization are the element density and the element angle. The optimization objective is to minimize the compliance \(c\) of the structure by updating the element density and the element angle, thereby optimizing the force-bearing model. The mathematical model for topology optimization is as follows:
[0030] find \(x = [x_1,x 2,..., x e,... x n T ,\(\theta = [\theta_1,\theta 2,..., \theta e,..., \theta n T
[0031] min \(c(x,\theta)=F T U(x,\theta)
[0032] subject \(KU = F
[0033]
[0034] 0\leq x min \leq x e \leq1
[0035] where \(x\) represents the element density matrix, \(x e is the density value of the \(e\)-th element, \(\theta\) represents the element angle matrix, \(\theta e is the angle value of the \(e\)-th element, \(n\) is the number of elements, \(c\) is the compliance of the structure, \(F\) is the vector of nodal forces, \(U\) is the vector of nodal displacements, \(K\) is the global stiffness matrix, which is related to the fiber angle arrangement; \(f\) is the volume ratio specified in the design domain, \(V e is the element volume, \(V_0\) is the volume of the design domain, \(x min is the minimum element density to avoid singularity problems in the global stiffness matrix \(K\); \(V\) is the optimized overall volume;
[0036] S13: Calculate the sensitivity of the compliance \(c\) of the optimized target structure to the design variables and perform sensitivity filtering;
[0037] The relationship between the compliance \(c\) of the structure and the parameters of each element is as follows:
[0038]
[0039] where \(u e is the displacement vector of element \(e\);
[0040] The sensitivity of the structural flexibility c with respect to the element density x e is:
[0041]
[0042] The sensitivity of the structural flexibility c with respect to the element angle θ e is:
[0043]
[0044] Sensitivity filtering is performed on the density, and the density filter function is:
[0045]
[0046] where γ is a positive number to avoid division by zero, and N e is the neighborhood of the element x e with volume v e The neighborhood is defined as:
[0047] N e = {j: dist(e,g) ≤ r min}
[0048] where dist(e,g) is the distance between the centers of element e and element g, and r min is the neighborhood size, and the weight factor H eg is:
[0049] H eg = r min - dist(e,g)
[0050] To make the adjacent element angle transition smooth, sensitivity filtering is performed on the element angle:
[0051]
[0052] S14: Update the element density and element angle respectively based on the filtered sensitivities through the optimality criterion and the gradient descent method until the convergence condition is reached, and output the optimal topology configuration and the discrete density and angle results.
[0053] The element density update adopts the Optimality Criteria (OC) method. Using the filtered sensitivity information of the objective function and constraints for the design variables, determine the adjustment direction and step size of the design variables, and introduce Lagrange multipliers to handle the constraints. Update the design variables iteratively to find the structural design that optimizes the objective function;
[0054] The unit angle is iteratively updated by the gradient descent method, and the conjugate mapping is performed on the sensitivity to accelerate the convergence speed. The step size factor h and the angle movement limit m are determined to control the speed of angle update. The update criterion is as follows:
[0055]
[0056] Optionally, the step S2 specifically includes the following steps:
[0057] S21: For the valid region in the unit density matrix x that is greater than or equal to 0.5, the valid angle data is extracted by constructing a mask matrix; then, the similarity between the valid angle data is calculated, and a similarity matrix is constructed based on the angle difference using the Gaussian kernel function. The expression of the Gaussian kernel function is:
[0058]
[0059] where K ij is the similarity between the i-th and j-th data points, θ i and θ j represent angle data, and σ is the bandwidth parameter of the kernel function; by calculating the Laplacian matrix of this similarity matrix, eigenvalue decomposition is performed, and the first K eigenvectors are selected;
[0060] S22: Using the first K eigenvectors, the K-means algorithm is used to divide the data into K clustering regions. Each region represents a clustering result, and the clustering result is mapped back to the original unit angle matrix θ according to the spatial position of the original angle data; in the K-means algorithm, data division is performed by minimizing the squared distance from each point to the center of its belonging cluster. The objective function J of the K-means algorithm is:
[0061]
[0062] where r ik is an indicator variable, and the indicator variable is used to represent whether the data point θ i belongs to cluster k, μ k is the center of cluster k, |θ i -μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of units;
[0063] S23: For each clustering region, calculate the mean of all angle values within the clustering region as the unified angle value of the clustering region, and update the angle values of all units in the clustering region to this unified angle value;
[0064] S24: Construct a binary mask for each angle value using the bwlabel function, perform connectivity analysis, and divide it into independent connected regions to ensure that each region has a unique region number;
[0065] S25: For each independent region, if it has a different number from the adjacent cells, generate a boundary point; then sort and connect based on the distance and direction between the boundary points to obtain the boundary line segments between adjacent regions.
[0066] Optionally, step S3 specifically includes the following steps:
[0067] S31: The unit density matrix x is a matrix of rows×cols. The grid coordinates are re-divided using the meshgrid function to ensure that the coordinates fall at the center of each grid. The calculation method is as follows:
[0068] X = (0.5:cols - 0.5), Y = (0.5:rows - 0.5)
[0069] X and Y represent the horizontal and vertical coordinates of the grid center points. Then, use the element values in the unit density matrix as the Z coordinate, that is, the density value of each grid point; then generate a four-dimensional control point matrix C, which contains the X, Y, and Z coordinates of the grid points, as well as the homogeneous coordinate, defaulting to 1. The elements of the four-dimensional control point matrix are expressed as:
[0070] C(i,j,1) = X ij ,C(i,j,2) = Y ij ,C(i,j,3) = Z ij ,C(i,j,4) = 1
[0071] S32: Use the NURBS surface fitting technique to smooth the jagged boundaries in the topology optimization result. The surface degrees in the u and v directions are both set to 3, and the knot vectors of the NURBS surface are generated. The specific knot vector formula is as follows:
[0072] N(u,v) = (N1(u,v),N2(u,v),dots,N t (u,v))
[0073] Among them, t is the number of knots, u and v are the parameter spaces; use the augknt function to generate the knot vectors in the u and v directions; then pass the control point matrix and the knot vectors into the nrbmak function to generate the fitted NURBS surface model, and perform visualization through the nrbplot function to display the fitted result in a three-dimensional coordinate system;
[0074] S33: Surface evaluation and plane intersection: To obtain boundary points that satisfy the volume constraint, set the plane intersection height to Z = 0.5 and evaluate the NURBS surface using a denser grid. The surface evaluation is achieved through the following formula:
[0075] P(u, v) = (X(u, v), Y(u, v), Z(u, v))
[0076] where P(u, v) is the coordinate value of the NURBS surface at parameters u and v. Then, find the points close to Z = 0.5 through the tolerance range, satisfying the following condition:
[0077] |Z ij - 0.5| < tolerance
[0078] These points will be used as candidate points for the intersection plane;
[0079] S34: Process the extracted points. First, calculate the distance between adjacent points to ensure that the connected line segments are not too long. Calculate the Euclidean distance between two points:
[0080]
[0081] If d ij is less than the set maximum distance threshold d max , then these two points are considered to belong to the same curve segment. This step ensures the removal of small and unnecessary details while retaining the valid line segments connecting adjacent points;
[0082] S35: Use the spline interpolation method to smoothly fit each valid curve segment. For each set of valid points, use the spline interpolation function to fit the curve. The interpolation formula is:
[0083]
[0084] where S(t) is the spline interpolation curve, B i (t) is the basis function, c i is the interpolation coefficient, and l is the number of interpolation points. By refining the interpolation points, a smoother curve is obtained. Finally, draw the fitted curve to ensure the smoothness of the boundary;
[0085] S36: Perform curve closing. Calculate the distance between the start and end points of the curve:
[0086]
[0087] If d closure is less than the set closing threshold d closure-threshold , then connect the start and end points to complete the closing operation, ensuring the integrity and connectivity of the boundary.
[0088] Optionally, step S4 specifically includes the following steps:
[0089] S41: Check whether the boundary line of each region obtained by K-means clustering is close to being closed. If the distance between the starting point and the ending point of the line segment is less than the set threshold, it is considered that the region is close to being closed and is directly stored as a closed region; if it is not closed, connect the starting point and the ending point of the boundary line to the nearest points on the boundary to gradually improve the boundary until the region is completely closed;
[0090] S42: Calculate the bounding box of each region. Determine the minimum rectangular frame of the region by obtaining the maximum and minimum horizontal and vertical coordinates of the boundary; based on the given angle and line width, determine the slope and intercept to generate a series of filling lines;
[0091] S43: Clip the filling lines to the outer boundary. If there is an inner boundary in the region, further clip the filling lines to ensure that the filling lines do not cross the inner boundary and avoid unreasonable geometric structures;
[0092] S44: After generating and clipping the filling lines, sort them according to the distance between the line segments to ensure that the filling line segments are arranged in the appropriate order, obtain the accurate starting and ending point coordinates, and reduce the empty running path of the printer;
[0093] S45: According to the coordinates of the points on the boundary line and the sorted filling lines, calculate the distance between each two points and the extrusion amount, and generate the corresponding Gcoode code in combination with the printer model and the used filament.
[0094] Based on the same inventive concept, the present invention discloses an additive manufacturing path preprocessing system for topology optimization design, which is characterized by including:
[0095] A topology optimization module, which is used to consider the anisotropy of carbon fiber reinforced composite materials during the topology optimization process, synchronously optimize the element density and element angle, that is, synchronously optimize the fiber orientation, and the output optimization results include the optimal topology structure diagram and the density and angle of each discrete element to ensure that the structural performance of the optimization results meets the actual requirements;
[0096] A demarcation line processing module, which is used to process the discrete angle result θ obtained during the optimization process by using the K-means clustering algorithm, divide it into multiple regions, and obtain the demarcation line of each region;
[0097] A boundary processing module, which is used to smooth the jagged boundary caused by mesh division in the optimal topology grid diagram by using the NURBS surface fitting technology;
[0098] The filling and coding module is used to generate the boundary line of each region by combining the smoothed boundary line and the region boundaries obtained by K-means clustering, and fill straight lines at corresponding angles within each region to generate corresponding Gcode.
[0099] Optionally, in the topology optimization module, the size of the design domain, the force model, the support conditions, and the material properties are first defined. The design domain is discretized into four-node rectangular elements. Based on the initial angles and initial densities of the elements, the initial element stiffness matrix and the global stiffness matrix are obtained.
[0100] The compliance constitutive matrix of the carbon fiber reinforced composite material is S0 and the stiffness constitutive matrix is C0 as follows:
[0101]
[0102] C0 = S0 -1
[0103] Among them, E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, ν 12 is the Poisson's ratio, G 12 is the in-plane shear modulus;
[0104] When the local fiber coordinate system forms an angle θ e with the global coordinate system, the stiffness constitutive matrix C(θ e ) of the e-th element is obtained by rotating the stiffness constitutive matrix C0 of the material:
[0105] C(θ e ) = T(θ e )C0T(θ e ) T
[0106] Among them, the rotation matrix is:
[0107]
[0108] The element stiffness matrix k e (θ e ) is:
[0109]
[0110] Among them, B is the strain-displacement matrix, and Ω e is the integration region of the element;
[0111] The assembled global stiffness matrix K(x, θ) is:
[0112]
[0113] Among them, p is the penalty factor;
[0114] Based on the solid isotropic material with penalization (SIMP) method, the design variables for topology optimization are the element density and the element angle. The optimization objective is to minimize the compliance c of the structure by updating the element density and the element angle, so as to optimize the force-bearing model. The mathematical model of topology optimization is as follows:
[0115] find x=[x1,x 2,..., x e,... x n T ,θ=[θ1,θ 2,..., θ e,..., θ n T
[0116] min c(x,θ)=F T U(x,θ)
[0117] subject KU=F
[0118]
[0119] 0≤x min ≤x e ≤1
[0120] where x represents the element density matrix, x e is the density value of the e-th element, θ represents the element angle matrix, θ e is the angle value of the e-th element, n is the number of elements, c is the structural compliance, F is the vector of nodal forces, U is the vector of nodal displacements, K is the global stiffness matrix, which is related to the fiber angle arrangement; f is the volume ratio specified in the design domain, V e is the element volume, V0 is the volume of the design domain, x min is the minimum element density, to avoid the singularity problem in the global stiffness matrix K; V is the optimized overall volume;
[0121] Calculate the sensitivity of the optimized objective structural compliance c to the design variables and perform sensitivity filtering;
[0122] The relationship between the structural compliance c and the parameters of each element is as follows:
[0123]
[0124] where u e is the displacement vector of element e;
[0125] The sensitivity of the structural compliance c to the element density x e is:
[0126]
[0127] The sensitivity of the structural flexibility c with respect to the element angle θ e is as follows:
[0128]
[0129] Perform sensitivity filtering on the density, and the density filter function is:
[0130]
[0131] where γ is a positive number to avoid division by zero, and N e is the neighborhood of the element x e with volume v e The neighborhood is defined as:
[0132] N e = {j: dist(e,g) ≤ r min}
[0133] where dist(e,g) is the distance between the center of element e and element g, and r min is the neighborhood size, and the weight factor H eg is:
[0134] H eg = r min - dist(e,g)
[0135] To make the transition of adjacent element angles smooth, perform sensitivity filtering on the element angles:
[0136]
[0137] Based on the filtered sensitivities, update the element density and element angles respectively by the optimality criterion and the gradient descent method until the convergence condition is reached, and output the optimal topological configuration and the discrete density and angle results.
[0138] The element density update uses the Optimality Criteria (OC) method. Utilize the filtered sensitivity information of the objective function and constraints for the design variables to determine the adjustment direction and step size of the design variables, and introduce Lagrange multipliers to handle the constraints. Update the design variables iteratively to find the structural design that optimizes the objective function;
[0139] The element angles are iteratively updated by the gradient descent method. Perform conjugate mapping on the sensitivities to accelerate the convergence speed, determine the step size factor h and the angle movement limit m to control the speed of angle update, and the update criterion is:
[0140]
[0141] Optionally, for the effective region in the unit density matrix x that is greater than or equal to 0.5, the boundary line processing module extracts the effective angle data by constructing a mask matrix; then calculates the similarity between the effective angle data, and constructs a similarity matrix based on the angle difference using the Gaussian kernel function. The expression of the Gaussian kernel function is:
[0142]
[0143] where K ij is the similarity between the i-th and j-th data points, θ i and θ j represent angle data, and σ is the bandwidth parameter of the kernel function; by calculating the Laplacian matrix of the similarity matrix, performing eigenvalue decomposition, and selecting the first K eigenvectors;
[0144] Using the first K eigenvectors, the K-means algorithm is used to divide the data into K clustering regions. Each region represents a clustering result, and the clustering result is mapped back to the original unit angle matrix θ according to the spatial position of the original angle data; in the K-means algorithm, data division is performed by minimizing the squared distance from each point to the center of its belonging cluster. The objective function J of the K-means algorithm is:
[0145]
[0146] where r ik is an indicator variable, and the indicator variable is used to represent whether the data point θ i belongs to cluster k, μ k is the center of cluster k, |θ i -μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of units;
[0147] For each clustering region, calculate the mean of all angle values within the clustering region as the unified angle value of the clustering region, and update the angle values of all units in the clustering region to this unified angle value;
[0148] Construct a binary mask for each angle value through the bwlabel function, and perform connectivity analysis to divide it into independent connected regions, ensuring that each region has a unique region number;
[0149] For each independent region, if it has a different number from the adjacent cells, generate a boundary point; then sort and connect based on the distance and direction between the boundary points to obtain the boundary line segment between adjacent regions.
[0150] An electronic device of the present invention includes a processor and a memory,
[0151] A memory for storing a computer program which is run by a processor to execute the above-mentioned preprocessing method for additive manufacturing paths for topology optimization design.
[0152] A computer storage medium of the present invention stores a computer program which is run by a processor to execute the above-mentioned preprocessing method for additive manufacturing paths for topology optimization design.
[0153] Beneficial effects: Compared with the prior art, the present invention has the following remarkable advantages: The present invention realizes the consistency between the optimization result and the actual manufacturing under the consideration of anisotropy and the angle constraint of additive manufacturing; The present invention solves the problem of the zigzag boundary in the optimization result through K-means clustering and NURBS surface fitting technologies, improving the printing accuracy and the reliability of the structural performance. Description of the Drawings
[0154] Figure 1 is a flowchart of the present invention;
[0155] Figure 2 is a structural schematic diagram of the design object of the present invention;
[0156] Figure 3 is a schematic diagram of the angle density synchronous optimization result of the design object of the present invention;
[0157] Figure 4 is a schematic diagram of the demarcation of the independent regions after K-means clustering in the post-processing of the design object of the present invention;
[0158] Figure 5 is a schematic diagram of the smooth structural boundary obtained by NURBS surface fitting in the post-processing of the design object of the present invention;
[0159] Figure 6 is a schematic diagram of the boundary of each region in the post-processing of the design object of the present invention;
[0160] Figure 7 is a schematic diagram of the corresponding angle straight filling lines in the post-processing of the design object of the present invention;
[0161] Figure 8 is a schematic diagram of the G-code generated in the post-processing of the design object of the present invention;
[0162] Figure 9 is a schematic diagram of the actual printing result in the post-processing of the design object of the present invention. Detailed Embodiments
[0163] The technical solution of the present invention will be further described below with reference to the drawings.
[0164] Embodiment 1
[0165] As Figure 1 shown, a pre - processing method for additive manufacturing paths for topology - optimized design includes the following steps:
[0166] S1: During the topology optimization process, considering the anisotropy of carbon fiber - reinforced composite materials, synchronously optimize the element density and element angle, that is, synchronously optimize the fiber orientation. The output optimization results include the optimal topology structure diagram and the density and angle of each discrete element to ensure that the structural performance of the optimization results meets the actual requirements.
[0167] Step S1 specifically includes the following steps:
[0168] S11: First, define the design domain size, force model, support conditions, and material properties. Discretize the design domain into four - node rectangular elements. According to the initial angle and initial density of the elements, obtain the initial element stiffness matrix and the global stiffness matrix;
[0169] The compliance constitutive matrix of the carbon fiber - reinforced composite material is S0 and the stiffness constitutive matrix is C0 as follows:
[0170]
[0171] C0 = S0 -1
[0172] where E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, v 12 is the Poisson's ratio, G 12 is the in - plane shear modulus;
[0173] When the local fiber coordinate system forms an angle θ e with the global coordinate system, the stiffness constitutive matrix C(θ e ) of the e - th element is obtained by rotating the stiffness constitutive matrix C0 of the material:
[0174] C(θ e ) = T(θ e )C0T(θ e ) T
[0175] where the rotation matrix is:
[0176]
[0177] The element stiffness matrix k e (θ e ) is:
[0178]
[0179] where B is the strain - displacement matrix, Ω eis the integration region of the element;
[0180] The overall stiffness matrix K(x,θ) is obtained by assembly as follows:
[0181]
[0182] where p is the penalty factor.
[0183] S12: Based on the Solid Isotropic Material with Penalization (SIMP) method, the design variables for topology optimization are the element density and the element angle. The optimization goal is to minimize the compliance c of the structure by updating the element density and the element angle, thereby optimizing the force-bearing model. The mathematical model of topology optimization is as follows:
[0184] find x=[x1,x 2,..., x e,... x n T ,θ=[θ1,θ 2,..., θ e,..., θ n T
[0185] min c(x,θ)=F T U(x,θ)
[0186] subject KU=F
[0187]
[0188] 0≤x min ≤x e ≤1
[0189] where x represents the element density matrix, and x e is the density value of the e-th element, θ represents the element angle matrix, and θ e is the angle value of the e-th element, n is the number of elements, c is the structural compliance, F is the vector of nodal forces, U is the vector of nodal displacements, K is the overall stiffness matrix, which is related to the fiber angle arrangement; f is the volume ratio specified in the design domain, V e is the element volume, V0 is the design domain volume, x min is the minimum element density, to avoid singular problems in the global stiffness matrix K; V is the optimized overall volume.
[0190] S13: Calculate the sensitivity of the optimized objective structural compliance c to the design variables. To avoid the checkerboard phenomenon, sensitivity filtering is performed;
[0191] The relationship between the structural compliance c and the parameters of each element is as follows:
[0192]
[0193] where u e is the displacement vector of element e;
[0194] The sensitivity of the structural compliance c with respect to the element density x e is:
[0195]
[0196] The sensitivity of the structural compliance c with respect to the element angle θ e is:
[0197]
[0198] Perform sensitivity filtering on the density, and the density filter function is:
[0199]
[0200] where γ = 10 -3 is a small positive number to avoid division by zero, N e is the neighborhood of the element x e with volume v e and the neighborhood is defined as:
[0201] N e = {j: dist(e, g) ≤ r min}
[0202] where dist(e, g) is the distance between the center of element e and element g, and r min is the neighborhood size, and the weight factor H eg is:
[0203] H eg = r min - dist(e, g)
[0204] To make the adjacent element angles transition smoothly, perform sensitivity filtering on the element angles:
[0205]
[0206] The penalty value p = 3 and the minimum filtering radius r min = 1.5 can be used for sensitivity filtering.
[0207] S14: Update the element density and element angles respectively based on the filtered sensitivities through the Optimality Criterion (OC) and the gradient descent method until the convergence condition is reached, and output the optimal topological configuration and the discrete density and angle results.
[0208] The unit density is updated by using the optimization criterion method (OC), which uses the filtered sensitivity information of the objective function and the constraints on the design variables to determine the adjustment direction and step size of the design variables, and introduces the Lagrange multiplier to deal with the constraints. The design variables are updated iteratively to find the structural design that optimizes the objective function.
[0209] The unit angle is iteratively updated by the gradient descent method, and the sensitivity is conjugate mapped to speed up the convergence speed, determine the step factor h and the angle movement limit m, and control the speed of angle update. The update criterion is:
[0210]
[0211] S2: The discrete angle result θ obtained during the optimization process is processed using the K-means clustering algorithm to divide it into multiple regions and obtain the boundary line of each region.
[0212] The step S2 specifically includes the following steps:
[0213] S21: For the valid area greater than or equal to 0.5 in the unit density matrix x, extract the valid angle data by constructing a mask matrix; then calculate the similarity between the valid angle data, and use the Gaussian kernel function to construct a similarity matrix based on the angle difference. The expression of the Gaussian kernel function is:
[0214]
[0215] Among them, K ij is the similarity between the i-th and j-th data points, θ i and θ j represents the angle data, σ is the bandwidth parameter of the kernel function; by calculating the Laplace matrix of the similarity matrix, eigenvalue decomposition is performed and the first K eigenvectors are selected.
[0216] S22: Using the first K eigenvectors, the K-means algorithm is used to divide the data into K clustering regions. Each region represents a clustering result. The clustering result is mapped back to the original unit angle matrix θ according to the spatial position of the original angle data. In the K-means algorithm, data division is performed by minimizing the square distance from each point to the center of its cluster. The objective function J of the K-means algorithm is:
[0217]
[0218] Among them, r ik is an indicator variable, which is used to represent the data point θ i Whether it belongs to cluster k, μk is the center of cluster k, |θ i - μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of cells.
[0219] S23: For each clustering region, calculate the mean of all angle values within the clustering region as the unified angle value of the clustering region, and update the angle values of all cells in the clustering region to this unified angle value.
[0220] S24: Construct a binary mask for each angle value through the bwlabel function, perform connectivity analysis, and divide it into independent connected regions to ensure that each region has a unique region number.
[0221] S25: For each independent region, if it has a different number from the adjacent cells (right, bottom, left, top), generate a boundary point; then sort and connect based on the distance and direction between the boundary points to obtain the boundary line segment between adjacent regions.
[0222] S3: For the jagged boundary caused by grid division in the optimal topological grid graph, use the NURBS (Non-Uniform Rational B-Spline) surface fitting technology to smooth the boundary.
[0223] Step S3 specifically includes the following steps:
[0224] S31: Grid coordinate redivision and generation of the control point matrix: The cell density matrix x is a matrix of rows×cols. The grid coordinates are redivided through the meshgrid function to ensure that the coordinates fall at the center of each grid. The calculation method is:
[0225] X = (0.5:cols - 0.5), Y = (0.5:rows - 0.5)
[0226] X and Y represent the horizontal and vertical coordinates of the grid center point. Then, use the element values in the cell density matrix as the Z coordinate, that is, the density value of each grid point; then generate a four-dimensional control point matrix C, which contains the X, Y, and Z coordinates of the grid point, as well as the homogeneous coordinate, defaulting to 1. The elements of the four-dimensional control point matrix are expressed as:
[0227] C(i,j,1) = X ij , C(i,j,2) = Y ij , C(i,j,3) = Z ij , C(i,j,4) = 1.
[0228] S32: Use the NURBS surface fitting technique to smooth the jagged boundaries in the topology optimization results. Set the surface degrees in both the u and v directions to 3, and generate the knot vectors of the NURBS surface. The specific knot vector formula is as follows:
[0229] N(u,v) = (N1(u,v), N2(u,v), dots, N t (u,v))
[0230] where t is the number of knots, and u and v are the parameter spaces; use the augknt function to generate the knot vectors in the u and v directions; then pass the control point matrix and knot vectors into the nrbmak function to generate the fitted NURBS surface model, and perform visualization through the nrbplot function to display the fitted results in a three-dimensional coordinate system.
[0231] S33: Surface evaluation and plane cutting: To obtain the boundary points that satisfy the volume constraint, set the plane cutting height to Z = 0.5, and evaluate the NURBS surface using a denser grid. The evaluation of the surface is achieved through the following formula:
[0232] P(u,v) = (X(u,v), Y(u,v), Z(u,v))
[0233] where P(u,v) is the coordinate value of the NURBS surface at parameters u and v, and then find the points close to Z = 0.5 through the tolerance range tolerance, satisfying the following condition:
[0234] |Z ij - 0.5| < tolerance
[0235] These points will be used as candidate points for the cutting plane.
[0236] S34: Distance calculation between points and valid point screening: Process the extracted points. First, calculate the distance between adjacent points to ensure that the connected line segments are not too long. Calculate the Euclidean distance between two points:
[0237]
[0238] If d ij is less than the set maximum distance threshold d max , then consider these two points to belong to the same curve segment. This step ensures the removal of small and unnecessary details while retaining the valid line segments connecting adjacent points.
[0239] S35: Use the spline interpolation method to perform smooth fitting on each valid curve segment. For each set of valid points, use the spline interpolation function to fit the curve. The interpolation formula is:
[0240]
[0241] Among them, S(t) is the spline interpolation curve, B i (t) is the basis function, c i is the interpolation coefficient, and l is the number of interpolation points; by refining the interpolation points (for example, generating 1000 points for each curve segment), a smoother curve is obtained; finally, the fitted curve is drawn to ensure the smoothness of the boundary.
[0242] S36: Curve closing process: To ensure a complete boundary curve, a curve closing process is required. Calculate the distance between the start and end points of the curve:
[0243]
[0244] If d closure is less than the set closing threshold d closure-threshold , then connect the start and end points to complete the closing operation, ensuring the integrity and connectivity of the boundary.
[0245] S4. Region and boundary line generation: Combine the smoothed boundary line and the region boundaries obtained by K-means clustering to generate the boundary lines of each region, and fill each region with straight lines at corresponding angles.
[0246] Step S4 specifically includes the following steps:
[0247] S41: Check whether the region boundary line obtained by each K-means clustering is close to being closed. If the distance between the start and end points of the line segment is less than the set threshold, it is considered that the region is close to being closed and is directly stored as a closed region; if it is not closed, connect the start and end points of the boundary line by searching for the nearest points to the boundary until the region is completely closed.
[0248] S42: Calculate the bounding box of each region. Determine the minimum rectangular frame of the region by obtaining the maximum and minimum horizontal and vertical coordinates of the boundary; based on the given angle and line width, determine the slope and intercept to generate a series of filling lines.
[0249] S43: Clip the filling lines to the outer boundary. If there is an inner boundary in the region, further clip the filling lines to ensure that the filling lines do not cross the inner boundary and avoid unreasonable geometric structures.
[0250] S44: After generating and clipping the filling lines, sort them according to the distance between the line segments to ensure that the filling line segments are arranged in the appropriate order, obtain the accurate start and end coordinates, and reduce the idle running path of the printer.
[0251] S45: Calculate the distance between each pair of points and the extrusion amount based on the coordinates of the points on the boundary line and the sorted filling line. Combine the printer model and the used filament to generate the corresponding Gcode.
[0252] As Figure 2 shown, in this embodiment, the design object is a cantilever beam model. The left end of the cantilever beam model is fixed. In the figure, P is the application point of the load, located at the lower right corner, with a magnitude of 1 KN. The structural domain has a length L = 80 mm and a width H = 50 mm. The design domain is discretized using an 80×50 four-node rectangular mesh, that is, the mesh size r is 1 mm; the volume constraint is 30%, and the iteration is 400 times.
[0253] As Figure 3 shown, the synchronous optimization result of the cantilever beam angle density. Due to angle discretization, actual manufacturing cannot guarantee consistency with the optimization result. The boundary line of the independent region after K-means clustering is as Figure 4 shown, the smooth structural boundary obtained by NURBS surface fitting is as Figure 5 shown. The boundary of each region is as shown in 6. After filling with straight lines at the corresponding angles, it is as shown in 7. The generated G-code and the actual printing result are as Figure 8 and Figure 9 shown.
[0254] Embodiment 2: Based on the same inventive concept, the present invention discloses an additive manufacturing path pre-processing system for topology optimization design, including:
[0255] A topology optimization module, which is used to consider the anisotropy of carbon fiber reinforced composite materials during the topology optimization process, perform synchronous optimization of unit density and unit angle, that is, synchronous optimization of fiber orientation, and the output optimization result includes the best topology structure diagram and the density and angle of each discrete unit to ensure that the structural performance of the optimization result meets the actual requirements.
[0256] In the topology optimization module, the size of the design domain, the force model, the support condition, and the material properties are first defined. The design domain is discretized into four-node rectangular elements. According to the initial angle and initial density of the elements, the initial element stiffness matrix and the global stiffness matrix are obtained.
[0257] The compliance constitutive matrix of the carbon fiber reinforced composite material is S0 and the stiffness constitutive matrix is C0 as follows:
[0258]
[0259] C0 = S0 -1
[0260] where E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, v12 is the Poisson's ratio, and G 12 is the in-plane shear modulus;
[0261] When the local coordinate system of the fiber makes an angle θ with the global coordinate system e , the stiffness constitutive matrix C(θ e ) of the e-th element is obtained by rotating the stiffness constitutive matrix C0 of the material:
[0262] C(θ e ) = T(θ e )C0T(θ e ) T
[0263] where the rotation matrix is:
[0264]
[0265] The element stiffness matrix k e (θ e ) is:
[0266]
[0267] where B is the strain-displacement matrix and Ω e is the integration region of the element;
[0268] The assembled global stiffness matrix K(x, θ) is:
[0269]
[0270] where p is the penalty factor;
[0271] Based on the Solid Isotropic Material with Penalization (SIMP) method, the design variables for topology optimization are the element density and the element angle. The optimization objective is to minimize the compliance c of the structure by updating the element density and the element angle, thereby optimizing the force-bearing model. The mathematical model for topology optimization is:
[0272] find x = [x1, x 2,..., x e,... x n T , θ = [θ1, θ 2,..., θ e,..., θ n T
[0273] min c(x, θ) = F T U(x, θ)
[0274] subject KU = F
[0275]
[0276] 0 ≤ x min ≤ x e ≤ 1
[0277] where x represents the element density matrix, and x e is the density value of the e-th element, θ represents the element angle matrix, and θ e is the angle value of the e-th element, n is the number of elements, c is the structural flexibility, F is the vector of nodal forces, U is the vector of nodal displacements, K is the global stiffness matrix, which is related to the fiber angle arrangement; f is the volume ratio specified in the design domain, V e is the element volume, V0 is the volume of the design domain, and x min is the minimum element density, to avoid the singularity problem in the global stiffness matrix K; V is the optimized overall volume.
[0278] Calculate the sensitivity of the optimized objective structural flexibility c to the design variables, and to avoid the checkerboard phenomenon, perform sensitivity filtering;
[0279] The relationship between the structural flexibility c and each element parameter is:
[0280]
[0281] where u e is the displacement vector of element e;
[0282] The sensitivity of the structural flexibility c to the element density x e is:
[0283]
[0284] The sensitivity of the structural flexibility c to the element angle θ e is:
[0285]
[0286] Perform sensitivity filtering on the density, and the density filter function is:
[0287]
[0288] where γ = 10 -3 is a small positive number to avoid division by zero, N e is the neighborhood of the element x e with volume v e , and the neighborhood is defined as:
[0289] N e={j: dist(e, g) ≤ r min}
[0290] where dist(e, g) is the distance between the center of unit e and unit g, and r min is the neighborhood size, and the weight factor H eg is as follows:
[0291] H eg = r min - dist(e, g)
[0292] To make the angle transition between adjacent units smooth, the sensitivity filtering of unit angles is performed:
[0293]
[0294] Based on the filtered sensitivity, the unit density and unit angle are updated respectively through the Optimality Criterion (OC) and the gradient descent method until the convergence condition is reached, and the optimal topological configuration and the discrete density and angle results are output.
[0295] The unit density update adopts the Optimality Criteria (OC) method. Using the filtered sensitivity information of the objective function and constraint conditions for the design variables, the adjustment direction and step size of the design variables are determined, and the Lagrange multiplier is introduced to handle the constraint conditions. The design variables are iteratively updated to find the structural design that optimizes the objective function.
[0296] The unit angle is iteratively updated through the gradient descent method. The sensitivity is conjugate mapped to accelerate the convergence speed. The step size factor h and the angle movement limit m are determined to control the speed of angle update. The update criterion is:
[0297]
[0298] The boundary line processing module is used to process the discrete angle result θ obtained in the optimization process by using the K-means clustering algorithm, divide it into multiple regions, and obtain the boundary lines of each region.
[0299] In the boundary line processing module, for the valid region in the unit density matrix x that is greater than or equal to 0.5, the valid angle data is extracted by constructing a mask matrix; then the similarity between the valid angle data is calculated, and a similarity matrix is constructed based on the angle difference using the Gaussian kernel function. The expression of the Gaussian kernel function is:
[0300]
[0301] where K ij is the similarity between the i-th and j-th data points, and θi and θ j represent angular data, and σ is the bandwidth parameter of the kernel function; by calculating the Laplacian matrix of the similarity matrix and performing eigenvalue decomposition, the first K eigenvectors are selected.
[0302] Using the first K eigenvectors, the K-means algorithm is used to divide the data into K clustering regions, and each region represents a clustering result. The clustering result is mapped back to the original unit angle matrix θ according to the spatial position of the original angular data; in the K-means algorithm, data division is performed by minimizing the squared distance from each point to the center of its belonging cluster, and the objective function J of the K-means algorithm is:
[0303]
[0304] where r ik is an indicator variable, and the indicator variable is used to represent whether the data point θ i belongs to cluster k, μ k is the center of cluster k, |θ i - μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of cells.
[0305] For each clustering region, calculate the mean value of all angular values within the clustering region as the unified angular value of the clustering region, and update the angular values of all cells in the clustering region to this unified angular value.
[0306] Construct a binary mask for each angular value through the bwlabel function, and perform connectivity analysis to divide it into independent connected regions to ensure that each region has a unique region number.
[0307] For each independent region, if it has a different number from the adjacent cells (right, bottom, left, top), generate a boundary point; then sort and connect based on the distance and direction between the boundary points to obtain the boundary line segment between adjacent regions.
[0308] Boundary processing module, used to smooth the jagged boundary caused by grid division in the optimal topological grid map by using NURBS surface fitting technology.
[0309] Grid coordinate redivision and generation of control point matrix in the boundary processing module: The unit density matrix x is a matrix of rows × cols. The grid coordinates are redivided through the meshgrid function to ensure that the coordinates fall at the center position of each grid, and the calculation method is:
[0310] X = (0.5:cols - 0.5), Y = (0.5:rows - 0.5)
[0311] X and Y represent the horizontal and vertical coordinates of the center point of the grid. Then, the element values in the cell density matrix are used as the Z coordinate, that is, the density value of each grid point. After that, a four-dimensional control point matrix C is generated, which contains the X, Y, and Z coordinates of the grid points, as well as the homogeneous coordinate, defaulting to 1. The elements of the four-dimensional control point matrix are expressed as:
[0312] C(i,j,1) = X ij ,C(i,j,2) = Y ij ,C(i,j,3) = Z ij ,C(i,j,4) = 1
[0313] The NURBS surface fitting technology is used to smooth the zigzag boundary in the topology optimization result. The surface degrees in the u and v directions are both set to 3, and the knot vectors of the NURBS surface are generated. The specific knot vector formula is as follows:
[0314] N(u,v) = (N1(u,v),N2(u,v),dots,N t (u,v))
[0315] where t is the number of knots, and u and v are the parameter spaces. The augknt function is used to generate the knot vectors in the u and v directions. Then, the control point matrix and the knot vectors are passed into the nrbmak function to generate the fitted NURBS surface model, and the nrbplot function is used for visualization to display the fitted result in a three-dimensional coordinate system.
[0316] Surface evaluation and plane cutting: To obtain the boundary points that meet the volume constraint, the plane cutting height is set to Z = 0.5, and a denser grid is used to evaluate the NURBS surface. The evaluation of the surface is achieved through the following formula:
[0317] P(u,v) = (X(u,v),Y(u,v),Z(u,v))
[0318] where P(u,v) is the coordinate value of the NURBS surface at the parameters u and v. Then, the points close to Z = 0.5 are found through the tolerance range tolerance, satisfying the following condition:
[0319] |Z ij - 0.5| < tolerance
[0320] These points will be used as the candidate points for the cutting plane.
[0321] Calculation of the distance between points and screening of valid points: The extracted points are processed. First, the distance between adjacent points is calculated to ensure that the connected line segments are not too long. The Euclidean distance between two points is calculated:
[0322]
[0323] If d ij is less than the set maximum distance threshold d max , then it is considered that these two points belong to the same curve segment. This step ensures the removal of small and unnecessary details while retaining the effective line segments connecting adjacent points.
[0324] Use the spline interpolation method to smoothly fit each valid curve segment. For each set of valid points, use the spline interpolation function to fit the curve. The interpolation formula is:
[0325]
[0326] where S(t) is the spline interpolation curve, B i (t) is the basis function, c i is the interpolation coefficient, and l is the number of interpolation points; by refining the interpolation points (such as generating 1000 points for each curve segment), a smoother curve is obtained; finally, the fitted curve is drawn to ensure the smoothness of the boundary.
[0327] Curve closing process: To ensure a complete boundary curve, a curve closing process is required. Calculate the distance between the start and end points of the curve:
[0328]
[0329] If d closure is less than the set closing threshold d closure-threshold , then connect the start and end points to complete the closing operation, ensuring the integrity and connectivity of the boundary.
[0330] The filling and coding module is used to combine the smoothed boundary line and the region boundaries obtained by K-means clustering to generate the boundary lines of each region, and fill each region with straight lines at corresponding angles to generate the corresponding Gcoode code.
[0331] In the filling and coding module, check whether the region boundary line obtained by each K-means clustering is close to being closed. If the distance between the start and end points of the line segment is less than the set threshold, it is considered that the region is close to being closed and is directly stored as a closed region; if it is not closed, then connect the start and end points of the boundary line by searching for the nearest points to the boundary to gradually improve the boundary until the region is completely closed.
[0332] Calculate the bounding box of each region. By obtaining the maximum and minimum horizontal and vertical coordinates of the boundary, determine the minimum rectangular frame of the region; based on the given angle and line width, determine the slope and intercept to generate a series of filling lines.
[0333] Trim the filling lines to the outer boundary. If there is an inner boundary in the area, further trim the filling lines to ensure that the filling lines do not cross the inner boundary and avoid unreasonable geometric structures.
[0334] After generating and trimming the filling lines, sort them according to the distance between line segments to ensure that the filling line segments are arranged in the appropriate order, obtain accurate starting and ending coordinates, and reduce the idle running path of the printer.
[0335] According to the coordinates of the points on the boundary line and the sorted filling lines, calculate the distance between each two points and the extrusion amount, and generate the corresponding Gcoode code in combination with the printer model and the used filament.
[0336] Embodiment 3
[0337] Corresponding to the method of Embodiment 1 of the present invention, Embodiment 3 of the present invention further provides an electronic device.
[0338] In this Embodiment 3, the electronic device includes: at least one communication bus, at least one processor, at least one memory, at least one network interface, and at least one peripheral interface. The memory contains programs and data.
[0339] The communication bus can be a communication device that transmits data between components inside the electronic device, such as an internal bus (CPU and memory bus), an external bus (universal serial bus port, peripheral component interconnect express port, etc.).
[0340] The memory may include high-speed RAM memory and may also include non-volatile memory, such as at least one disk memory.
[0341] The processor calls the programs and data stored in the memory to execute a preprocessing method for additive manufacturing paths for topology optimization design provided in Embodiment 1 of the present invention.
[0342] The peripheral interface is used to connect to peripherals. Peripherals are external devices, which may include but are not limited to keyboards, displays, cursor control devices (such as mice, touchpads, or touchscreens), video input devices, etc.
[0343] The network interface provides wired or wireless communication related to an external network (such as the Internet, intranet, local area network, mobile communication network, etc.).
[0344] Embodiment 4
[0345] Corresponding to the method of Embodiment 1 of the present invention, Embodiment 4 of the present invention further provides a computer storage medium for data collection and reception. The computer storage medium stores a computer program, which is run by a processor to execute an additive manufacturing path preprocessing method for topology optimization design provided in Embodiment 1 of the present invention.
[0346] In each embodiment of the present invention, each functional unit can be integrated in a processing unit, or each unit can exist physically alone, or two or more units can be integrated in one unit. The above-mentioned integrated unit can be implemented in the form of hardware, or in the form of a software functional unit, or in the form of a combination of software and hardware.
[0347] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in each embodiment of the present invention. The foregoing storage medium includes: various media that can store program codes, such as a mobile hard disk, a USB flash drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disc.
Claims
1. An additive manufacturing path preprocessing method for topology optimization design, characterized in that It includes the following steps: S1: During the topology optimization process, considering the anisotropy of carbon fiber reinforced composite materials, synchronously optimize the element density and element angle. The output optimization results include the optimal topology structure diagram, as well as the density and angle of each discrete element; S2: For the discrete angle result θ obtained during the optimization process, use the K-means clustering algorithm to process it, divide it into multiple regions, and obtain the boundary line of each region; S3: For the zigzag boundary caused by mesh division in the optimal topology grid diagram, use the NURBS surface fitting technology to smooth the boundary; S4: Combine the smoothed boundary line and the region boundaries obtained by K-means clustering to generate the boundary line of each region, and fill straight lines with corresponding angles in each region to generate the corresponding Gcoode code.
2. The pre - processing method for additive manufacturing paths for topology - optimization design according to claim 1, characterized in that, The specific steps of step S1 include the following steps: S11: First, define the design domain size, force model, support conditions, and material properties. Discretize the design domain into four-node rectangular elements, and obtain the initial element stiffness matrix and the global stiffness matrix according to the initial angle and initial density of the elements; The compliance constitutive matrix of carbon fiber reinforced composite materials is S0 and the stiffness constitutive matrix is C0 as follows: Among them, E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, v 12 is the Poisson's ratio, and G 12 is the in-plane shear modulus; When the local coordinate system of the fiber makes an angle θ with the global coordinate system e the stiffness constitutive matrix C(θ e ) of the e-th element is obtained by rotating the stiffness constitutive matrix C0 of the material as follows: C(θ e ) = T(θ e )C0T(θ e ) T Among them, the rotation matrix is: Element stiffness matrix k e (θ e ) is as follows: where B is the strain-displacement matrix and Ω e is the integration domain of the element; The assembled global stiffness matrix K(x,θ) is: Among them, p is the penalty factor; S12: Based on the solid isotropic material penalization method, determine that the design variables of topology optimization are element density and element angle. The optimization goal is to minimize the compliance c of the structure by updating the element density and element angle, so as to optimize the force model. The mathematical model of topology optimization is: find x=[x1,x 2,..., x e,... x n T ,θ=[θ1,θ 2,..., θ e,..., θ n T min c(x,θ)=F T U(x,θ) subject to KU = F 0 ≤ x min ≤ x e ≤ 1 where x represents the element density matrix, and x e is the density value of the e-th element, θ represents the element angle matrix, and θ e is the angle value of the e-th element, n is the number of elements, c is the structural flexibility, F is the vector of nodal forces, U is the vector of nodal displacements, K is the global stiffness matrix, which is related to the fiber angle arrangement; f is the volume ratio specified in the design domain, V e is the element volume, V0 is the volume of the design domain, x min is the minimum element density to avoid the singularity problem in the global stiffness matrix K; V is the optimized overall volume; S13: Calculate the sensitivity of the compliance c of the optimized target structure to the design variables, and perform sensitivity filtering; The relationship between the compliance c of the structure and the parameters of each element is: where u e is the displacement vector of element e; The sensitivity of the structural flexibility c with respect to the element density x e is as follows: The sensitivity of the structural flexibility c with respect to the element angle θ e is as follows: Perform sensitivity filtering on the density, and the density filter function is: where γ is a positive number to avoid division by zero, N e is the neighborhood of an element x e with volume v e and the neighborhood is defined as: N e = {j: dist(e, g) ≤ r min} where dist(e,g) is the distance between the center of cell e and cell g, and r min is the neighborhood size, and the weight factor H eg is as follows: H eg = r min - dist(e, g) In order to make the transition of adjacent element angles smooth, perform sensitivity filtering on the element angles: S14: Based on the filtered sensitivities, update the element density and element angle respectively through the optimality criterion and the gradient descent method until the convergence condition is reached, and output the optimal topology configuration and the discrete density and angle results. The element density is updated using the optimality criteria method (OC). Using the filtered sensitivity information of the objective function and constraints for the design variables, determine the adjustment direction and step size of the design variables, and introduce Lagrange multipliers to handle the constraints. Iteratively update the design variables to find the structural design that optimizes the objective function; The element angle is iteratively updated by the gradient descent method. Perform conjugate mapping on the sensitivity to accelerate the convergence speed. Determine the step size factor h and the angle movement limit m to control the speed of angle update. The update criterion is:
3. The pre - processing method for additive manufacturing paths for topology - optimization design according to claim 1, wherein The specific steps of step S2 include the following steps: S21: For the valid region in the unit density matrix x that is greater than or equal to 0.5, extract the valid angular data by constructing a mask matrix; then calculate the similarity between the valid angular data, and construct a similarity matrix based on the angular difference using the Gaussian kernel function. The expression of the Gaussian kernel function is: Among them, K ij is the similarity between the i-th and j-th data points, θ i and θ j represent angular data, and σ is the bandwidth parameter of the kernel function; by calculating the Laplacian matrix of the similarity matrix, performing eigenvalue decomposition, and selecting the top K eigenvectors; S22: Use the first K eigenvectors and the K-means algorithm to divide the data into K clustering regions. Each region represents a clustering result, and the clustering result is mapped back to the original unit angular matrix θ according to the spatial position of the original angular data; in the K-means algorithm, the data is divided by minimizing the squared distance from each point to the center of its belonging cluster. The objective function J of the K-means algorithm is: where r ik is an indicator variable that is used to indicate whether the data point θ i belongs to cluster k, μ k is the center of cluster k, |θ i - μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of cells; S23: For each clustering region, calculate the mean of all angular values within the clustering region as the unified angular value of the clustering region, and update the angular values of all units in the clustering region to this unified angular value; S24: Construct a binary mask for each angular value through the bwlabel function, and perform connectivity analysis to divide it into independent connected regions, ensuring that each region has a unique region number; S25: For each independent region, if it has a different number from the adjacent cells, generate a boundary point; then sort and connect based on the distance and direction between the boundary points to obtain the boundary line segment between adjacent regions.
4. The preprocessing method for additive manufacturing path for topology optimization design according to claim 1, characterized in that The specific steps of step S3 are as follows: S31: The unit density matrix x is a matrix of rows×cols. The grid coordinates are re-divided through the meshgrid function to ensure that the coordinates fall at the center of each grid. The calculation method is: X = (0.5:cols - 0.5), Y = (0.5:rows - 0.5) X and Y represent the horizontal and vertical coordinates of the grid center point, and then use the element values in the unit density matrix as the Z coordinate, that is, the density value of each grid point; then generate a four-dimensional control point matrix C, which contains the X, Y, and Z coordinates of the grid points, as well as the homogeneous coordinate, which is defaulted to 1. The elements of the four-dimensional control point matrix are expressed as: C(i,j,1) = X ij , C(i,j,2) = Y ij , C(i,j,3) = Z ij , C(i,j,4) = 1 S32: Use the NURBS surface fitting technology to smooth the zigzag boundary in the topology optimization result. The surface degrees in the u and v directions are both set to 3, and the knot vectors of the NURBS surface are generated. The specific knot vector formula is as follows: N(u, v) = (N1(u, v), N2(u, v), dots, N t (u, v)) where t is the number of knots, and u and v are the parameter spaces; use the augknt function to generate the knot vectors in the u and v directions; then pass the control point matrix and the knot vectors into the nrbmak function to generate the fitted NURBS surface model, and perform visualization through the nrbplot function to display the fitted result in a three-dimensional coordinate system; S33: Surface evaluation and plane cutting: To obtain the boundary points that meet the volume constraint, set the plane cutting height to Z = 0.5, and evaluate the NURBS surface using a denser grid. The evaluation of the surface is achieved through the following formula: P(u, v) = (X(u, v), Y(u, v), Z(u, v)) Among them, P(u, v) is the coordinate value of the NURBS surface at parameters u and v. Then, points close to Z = 0.5 are found through the tolerance range tolerance, satisfying the following conditions: |Z ij -0.5 | < tolerance These points will be used as candidate points for the cutting plane; S34: Process the extracted points. First, calculate the distance between adjacent points to ensure that the connected line segments are not too long. Calculate the Euclidean distance between two points: If d ij is less than the set maximum distance threshold d max , it is considered that these two points belong to the same curve segment. This step ensures the removal of small and unnecessary details while retaining the effective line segments connecting adjacent points; S35: Use the spline interpolation method to perform smooth fitting on each valid curve segment. For each set of valid points, use the spline interpolation function to fit the curve. The interpolation formula is: Among them, S(t) is the spline interpolation curve, B i (t) is the basis function, c i is the interpolation coefficient, and l is the number of interpolation points; by refining the interpolation points, a smoother curve is obtained; finally, the fitted curve is drawn to ensure the smoothness of the boundary; S36: Perform curve closing processing. Calculate the distance between the start and end points of the curve: If d closure is less than the set closing threshold d closure-threshold , the head and tail points are connected to complete the closing operation, ensuring the integrity and connectivity of the boundary.
5. The pre - processing method for additive manufacturing path for topology - optimization design according to claim 1, characterized in that, The specific steps of step S4 are as follows: S41: Check whether the region boundary line obtained by each K-means clustering is close to being closed. If the distance between the start and end points of the line segment is less than the set threshold, it is considered that the region is close to being closed and is directly stored as a closed region; if it is not closed, connect the start and end points of the boundary line to the nearest points on the boundary by searching, and gradually improve the boundary until the region is completely closed; S42: Calculate the bounding box of each region. Determine the minimum rectangular frame of the region by obtaining the maximum and minimum horizontal and vertical coordinates of the boundary; based on the given angle and line width, determine the slope and intercept, and generate a series of filling lines; S43: Clip the filling lines to the outer boundary. If there is an inner boundary in the region, further clip the filling lines to ensure that the filling lines do not cross the inner boundary and avoid unreasonable geometric structures; S44: After generating and clipping the filling lines, sort them according to the distance between the line segments to ensure that the filling line segments are arranged in the appropriate order, obtain the accurate start and end coordinates, and reduce the empty running path of the printer; S45: According to the coordinates of the points on the boundary line and the sorted filling lines, calculate the distance and extrusion amount between each two points, and generate the corresponding Gcoode code in combination with the printer model and the used filament.
6. An additive manufacturing path preprocessing system for topology optimization design, characterized in that Including: A topology optimization module, which is used to consider the anisotropy of carbon fiber reinforced composites during the topology optimization process, perform synchronous optimization of the element density and element angle, that is, synchronous optimization of the fiber orientation. The output optimization results include the best topology structure diagram and the density and angle of each discrete element to ensure that the structural performance of the optimization results meets the actual requirements; A boundary line processing module, which is used to process the discrete angle result θ obtained during the optimization process by using the K-means clustering algorithm, divide it into multiple regions, and obtain the boundary line of each region; A boundary processing module, which is used to smooth the jagged boundary caused by mesh division in the best topology mesh diagram by using the NURBS surface fitting technology; A filling code generation module, which is used to generate the boundary line of each region in combination with the smoothed boundary line and the region boundary obtained by K-means clustering, and fill straight lines with corresponding angles in each region to generate the corresponding Gcoode code.
7. The additive manufacturing path preprocessing system for topology optimization design according to claim 6, wherein, In the topological optimization module, the size of the design domain, the force model, the support conditions, and the material properties are first defined. The design domain is discretized into four-node rectangular elements. Based on the initial angles and initial densities of the elements, the initial element stiffness matrix and the global stiffness matrix are obtained. The compliance constitutive matrix of the carbon fiber reinforced composite material is S0 and the stiffness constitutive matrix is C0 as follows: Among them, E1 is the elastic modulus along the fiber direction, E2 is the elastic modulus perpendicular to the fiber direction, v 12 is the Poisson's ratio, and G 12 is the in-plane shear modulus; When the local coordinate system of the fiber makes an angle θ with the global coordinate system e the stiffness constitutive matrix C(θ e ) of the e-th element is obtained by rotating the stiffness constitutive matrix C0 of the material as follows: C(θ e ) = T(θ e )C0T(θ e ) T Among them, the rotation matrix is: Element stiffness matrix k e (θ e ) is as follows: where B is the strain-displacement matrix and Ω e is the integration domain of the element; The assembled global stiffness matrix K(x,θ) is: Among them, p is the penalty factor; Based on the solid isotropic material penalization method, the design variables for topological optimization are the element density and the element angle. The optimization goal is to minimize the compliance c of the structure by updating the element density and the element angle, so as to optimize the force model. The mathematical model of topological optimization is: find x=[x1,x2,...,x e ,...x n T ,θ=[θ1,θ2,...,θ e ,...,θ n T min c(x,θ)=F T U(x,θ) subject KU=F 0 ≤ x min ≤ x e ≤ 1 where x represents the element density matrix, and x e is the density value of the e-th element, θ represents the element angle matrix, and θ e is the angle value of the e-th element, n is the number of elements, c is the structural flexibility, F is the vector of nodal forces, U is the vector of nodal displacements, K is the global stiffness matrix, which is related to the fiber angle arrangement; f is the volume ratio specified in the design domain, V e is the element volume, V0 is the volume of the design domain, and x min is the minimum element density to avoid the singularity problem in the global stiffness matrix K; V is the optimized overall volume; Calculate the sensitivity of the compliance c of the optimized target structure to the design variables and perform sensitivity filtering. The relationship between the compliance c of the structure and the parameters of each element is: where u e is the displacement vector of element e; The sensitivity of the structural flexibility c with respect to the element density x e is as follows: The sensitivity of the structural flexibility c with respect to the element angle θ e is as follows: Perform sensitivity filtering on the density. The density filter function is: where γ is a positive number to avoid division by zero, N e is the neighborhood of an element x e with volume v e and the neighborhood is defined as: N e = {j: dist(e, g) ≤ r min} where dist(e,g) is the distance between the center of unit e and unit g, and r min is the neighborhood size, and the weight factor H eg is as follows: H eg = r min - dist(e, g) In order to make the transition of adjacent element angles smooth, perform sensitivity filtering on the element angles: Based on the filtered sensitivities, update the element density and the element angle respectively by the optimality criterion and the gradient descent method until the convergence condition is reached, and output the optimal topological configuration and the discrete density and angle results. The element density update adopts the Optimality Criteria (OC) method. Using the filtered sensitivity information of the objective function and the constraint conditions for the design variables, determine the adjustment direction and step size of the design variables, and introduce the Lagrange multiplier to handle the constraint conditions. Iteratively update the design variables to find the structural design that optimizes the objective function. The element angle is iteratively updated by the gradient descent method. Perform conjugate mapping on the sensitivity to accelerate the convergence speed. Determine the step size factor h and the angle movement limit m to control the speed of angle update. The update criterion is:
8. An additive manufacturing path preprocessing system for topology optimization design according to claim 6, characterized in that, In the boundary line processing module, for the effective region in the element density matrix x that is greater than or equal to 0.5, extract the effective angle data by constructing a mask matrix; then calculate the similarity between the effective angle data, and use the Gaussian kernel function to construct a similarity matrix based on the angle difference. The expression of the Gaussian kernel function is: where K ij is the similarity between the i-th and j-th data points, and θ i and θ j represent angular data, and σ is the bandwidth parameter of the kernel function; by calculating the Laplacian matrix of the similarity matrix, performing eigenvalue decomposition, and selecting the top K eigenvectors; Using the first K eigenvectors, use the K-means algorithm to divide the data into K clustering regions. Each region represents a clustering result. The clustering result is mapped back to the original element angle matrix θ according to the spatial position of the original angle data; in the K-means algorithm, the data is divided by minimizing the squared distance from each point to the center of its cluster. The objective function J of the K-means algorithm is: where r ik is an indicator variable, and the indicator variable is used to represent whether the data point θ i belongs to the cluster k, μ k is the center of the cluster k, |θ i - μ k | is the Euclidean distance from the data point to the cluster center, and n is the number of cells; For each clustering region, calculate the mean value of all angle values in the clustering region as the unified angle value of the clustering region, and update the angle values of all elements in the clustering region to this unified angle value. Construct a binary mask for each angle value through the bwlabel function and perform connectivity analysis to divide it into independent connected regions to ensure that each region has a unique region number. For each independent region, if it has a different number from the adjacent cells, generate a demarcation point; then sort and connect based on the distance and direction between the demarcation points to obtain the demarcation line segment between adjacent regions.
9. An electronic device, characterized in that, including a processor and a memory, The memory is used to store a computer program, and the computer program is run by the processor to execute the method according to any one of claims 1-5.
10. A computer storage medium, characterized in that, The computer storage medium stores a computer program, and the computer program is run by the processor to execute the method according to any one of claims 1-5.
Citation Information
Cited By
Topological optimization method and device based on CutFEM and SIMP, medium and equipment
CN121234745A
A topology optimization method, device, medium and equipment based on CutFEM and SIMP
CN121234745B