Body subdivision-based bidirectional progressive isogeometric topological optimization method

Abstract

The invention discloses a two-way progressive isogeometric topology optimization method based on volume subdivision, which comprises the following steps of: firstly, initializing the density of all hexahedron units input into a hexahedron mesh entity model M, and carrying out k1-time volume subdivision on the M to obtain an optimization model M1; secondly, constructing a three-variable Bezier spline body, applying an external load on the M to the corresponding spline body, and solving a displacement matrix U through isogeometric analysis; and constructing an initial stiffness matrix of the Bezier spline body corresponding to each subdivision unit in the M1. Finally, judging according to the target volume fraction and the convergence formula, and if the target volume fraction and the convergence formula are met, ending iteration to obtain a topological optimization model; otherwise, updating the volume fraction of the next iteration, calculating and optimizing the unit sensitivity of the M1, updating the density of each unit of the M1, and repeating the iteration. According to the method, isogeometric topological optimization is efficiently carried out, and a topological optimization result expressed by a spline with a clear boundary is obtained.

Description

Technical Field

[0001] The invention belongs to the application field of isogeometric analysis and structural optimization, and relates to a bidirectional progressive topology optimization method based on a variable deletion rate, and specifically to a bidirectional progressive isogeometric topology optimization method based on volume subdivision. Background Art

[0002] Isogeometric analysis is a new type of numerical calculation method that unifies design and analysis. Isogeometric analysis can express the geometric features of the model more accurately and with higher precision. It makes the seamless integration of CAD / CAE possible, while improving the calculation accuracy and efficiency.

[0003] Structural optimization refers to the use of various methods and techniques to achieve the optimal shape, material distribution and load capacity of the structure while meeting performance requirements when designing and analyzing engineering structures. Its purpose is to improve the efficiency of the structure, reduce materials and costs, and ensure the safety and stability of the structure. The results of traditional structural optimization require some post-processing operations before they can be imported into the CAD system, so seamless integration of design and analysis and handling of complex constraints are also very challenging issues for structural optimization.

[0004] However, the existing approximate volume subdivision algorithm will cause the model shape to deviate more and more from the initial mesh after multiple iterations. In addition, the SIMP topology optimization method has problems such as unclear boundaries and intermediate density units, which cannot meet the requirements of volume subdivision for clear structural boundaries.

[0005] Therefore, the present invention provides a bidirectional progressive topology optimization method based on interpolation body subdivision, which utilizes the input hexahedral model to efficiently obtain a topology optimization result with clear boundaries and can be directly imported into the CAD system. Summary of the invention

[0006] The purpose of the present invention is to realize the integration of geometric representation and simulation analysis, and to improve the efficiency of the calculation method and the accuracy of the result model, and to propose a variable deletion rate bidirectional progressive topology optimization method based on interpolation volume subdivision. The problem of unclear structure and the existence of intermediate density units in the SIMP topology optimization method is solved. At the same time, the interpolation volume subdivision method can generate a high-quality topology optimization model that can be directly imported into the CAD system.

[0007] The specific steps of the present invention are as follows:

[0008] Step 1: Input the original hexahedral mesh solid model M and initialize the density of all hexahedral units in the solid model M.

[0009] Step 2: Perform k1 volume subdivisions on M to obtain the optimized model M1, k1 ≥ 1. After each volume subdivision, the density of the subdivided unit inherits the density of the hexahedral unit at the previous level, and the volume fraction of the optimized model M1 before entering the iteration is initialized.

[0010] Step 3: According to the volume subdivision limit point formula, a three-variable Bézier spline body is constructed as an approximation of the limit subdivision body of M, and the external load on the model M is applied to its corresponding Bézier spline body.

[0011] Step 4: According to the global stiffness matrix and right-hand side term of the Bézier spline body after the load is applied, use the global stiffness matrix and right-hand side term to perform isogeometric analysis to solve the displacement matrix U. Construct the initial stiffness matrix K of the Bézier spline body corresponding to each subdivision unit in the optimization model M1 i,0 .

[0012] Step 5: Determine the volume fraction V of the current iterative optimization model M1 i Has the target volume fraction V been met? * , if the target volume fraction V is not met * Then go to step 6. If the target volume fraction V is met * , then further determine whether the convergence formula is satisfied. If not, execute step 6. If the convergence formula is satisfied, end the iteration and execute step 9.

[0013] Step 6: If the target volume V of the current iteration i With V * When they are equal, the target volume of subsequent iterations is fixed to V * Otherwise, update the volume fraction V for the next iteration i+1 .

[0014] Step 7: Use the initial stiffness matrix K obtained in step 4 i,0 And the displacement matrix U is used to calculate the sensitivity of each unit of the optimization model M1.

[0015] Step 8: Use the current volume fraction V i , target volume fraction V * The density of each unit in the optimized model M1 is updated based on the sensitivity of each unit in the optimized model M1, and the number of iterations is increased by 1 and the process returns to step 4.

[0016] Step 9: Delete all the units with the smallest density in the optimized model M1, and then construct a Bézier spline body for each remaining unit to obtain the final topology optimization model of the solid model.

[0017] Preferably, in step 2, the subdivided solid model corresponding to the hexahedral mesh solid model is generated using the following interpolation subdivision rule:

[0018] Step 2.1: First, assume that there is an internal point C in a model, and its neighborhood points can be represented by three types of points: D, K, and B. Point D is a point directly connected to C, point K is a point not directly connected to C but on the same surface, and point B is a point not directly connected to C and not on the same surface but on the same body. Now follow the steps below to give the virtual point calculation formulas for B, K, D, and C, and denote their virtual points as V B 、V K 、V D 、V C .

[0019] (1)V B Keeping point B unchanged, we have V B =B.

[0020] (2) V K When V K When replacing point K, C, V K The corresponding points in the two adjacent bodies whose faces are common faces construct a point set. When the point set is subdivided by three uniform B-spline subdivisions, V K The corresponding cubic uniform B-spline subdivision limit point coincides with the K point.

[0021] (3)V D When V D When replacing point D, D In the area corresponding to the four hexahedrons whose edges are common edges, V D When performing Catmull-Clark tessellation for the center of the surface, V D The corresponding limit point coincides with point D.

[0022] (4)V C When V C Instead of point C, at point V C When performing Catmull-Clark body subdivision, V C The corresponding limit point coincides with point C.

[0023] Step 2.2, after completing the calculation of the virtual points of points B, K, D, and C, the virtual points corresponding to each point in the adjacent domain point set of point C are obtained. Then, the Catmull-Clark volume subdivision rule can be used to give the rule for calculating the contribution values ​​of the volume points, surface points, and edge points added to the area corresponding to each virtual point:

[0024] (1) For each neighbor of point C, let the set of virtual points contained in the current neighbor be {p c1 ,p c2 ,…,p c8}, then the contribution of the virtual point to the point of the hexahedron is the average value of the virtual points in the set.

[0025] (2) For each adjacent face of point C, when the cubic uniform B-spline subdivision algorithm is applied to each point in the face, the contribution value is the average value of the limit points corresponding to these points.

[0026] (3) For each adjacent edge of point C, the contribution value is calculated as the average value of the corresponding limit points of the two points when Catmull-Clark surface subdivision is applied to the two endpoints of the edge.

[0027] Step 2.3, after completing the above process with the point inside each body as the center, the contribution of each body, surface, and edge point is summed and averaged to obtain the newly inserted body point, surface point, and edge point. Then, according to the topological rules of Catmull-Clark body subdivision, the topological connection is performed to obtain the optimized model after subdivision once.

[0028] Preferably, the formula for the volume subdivision limit point in step 3 is as follows:

[0029]

[0030] Among them, v 1 Represents the coordinates of vertex v after one subdivision. Similarly, f j 1 , They represent the coordinates of the jth edge point, face point, and block point of point v after one subdivision, and n represents the number of edges connected to vertex v. Using the above limit point formula, the generated Bézier body collective interpolation hexahedral mesh model M can be directly subdivided to generate the Catmull-Clark subdivision limit body.

[0031] Preferably, the formula for judging convergence in step 5 is as follows:

[0032]

[0033] Among them, C k-i+1 is the flexibility value at the k-i+1th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of previous and next iterations involved in judging convergence.

[0034] Preferably, the volume fraction updating strategy in step 6 is as follows:

[0035] Step 6.1, calculate the current iteration evolution ratio ER:

[0036]

[0037] Among them, ER max , ER minare the maximum and minimum evolution ratios, respectively.

[0038] Step 6.2: Calculate the volume fraction for the next iteration based on the current evolution ratio:

[0039] V i+1 =V i (1-ER)

[0040] As a preferred embodiment, the specific calculation process of sensitivity in step 7 is:

[0041] Step 7.1, use the following formula to calculate the initial unit sensitivity of the current iteration:

[0042]

[0043] Among them, u i To optimize the displacement of the i-th unit in the model M1 in the displacement matrix U, ρ i is the unit density, K i,0 is the initial stiffness matrix of the unit, and p is the penalty factor.

[0044] Step 7.2: Use the sensitivity filter function to filter each unit:

[0045]

[0046] Where M is the filter radius r min The number of all units in the range, α j0 is the unit sensitivity within the filtering range, r ij is the distance from the center of unit i to the center of unit j, w(r ij ) is the weight factor, and its value is r min -r ij .

[0047] Step 7.3: Combine the unit sensitivity of the previous iteration and average the filtered unit sensitivity to obtain the final unit sensitivity: Where k represents the current iteration number.

[0048] Preferably, the cell density update strategy in step 8 is as follows:

[0049] Step 8.1: Sort all unit sensitivities in descending order, denoted by α max is the maximum unit sensitivity of the current iteration, α min The minimum cell sensitivity of the current iteration is set to ρ, and the minimum cell density is set to min , indicating an empty cell.

[0050] Step 8.2: Combine α max , α min and volume fraction V i, using a binary strategy to find the sensitivity threshold α th .

[0051] Step 8.3, Sensitivity α i ≤a th The cell density is updated from 1 to ρ min , sensitivity α i >α th The cell density is changed from ρ min Updated to 1.

[0052] The present invention has the following characteristics and beneficial effects:

[0053] By adopting the above scheme, the bidirectional progressive topology optimization method based on the "soft kill method" is combined with the isogeometric topology optimization framework based on interpolation volume subdivision and unstructured three-variable spline volume construction method, which realizes efficient isogeometric topology optimization of the input hexahedral mesh and obtains the topology optimization results with clear boundaries and spline expression that can be directly imported into the CAD system. The quality and efficiency of the existing optimization results in the multi-resolution isogeometric topology optimization framework are improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 It is a flow chart of the present invention;

[0055] Figure 2(a) is a schematic diagram of the input hexahedral mesh model;

[0056] Figure 2(b) is a schematic diagram of the force input of the hexahedral mesh model;

[0057] Figure 3(a) shows the virtual point V K Schematic diagram of calculation;

[0058] Figure 3(b) shows the virtual point V D Schematic diagram of calculation;

[0059] Figure 3(c) is a schematic diagram of a virtual point set of a hexahedral grid;

[0060] Figure 3(d) is a schematic diagram of the contribution value of the virtual point set to the hexahedron points;

[0061] Figure 3(e) is a schematic diagram of the contribution value of the virtual point set to the hexahedron surface points;

[0062] Figure 3(f) is a schematic diagram of the contribution value of the virtual point set to the edge points of the hexahedron;

[0063] Figure 4(a) shows the mesh model of the isogeometric analysis and topology optimization results on M when the target volume fraction is 50%;

[0064] Figure 4(b) shows the spline model of isogeometric analysis and topology optimization results on M when the target volume fraction is 50%;

[0065] Figure 4(c) shows the mesh model of the topology optimization results performed on the optimized model M1 after performing isogeometric analysis on M when the target volume fraction is 50%;

[0066] Figure 4(d) shows the spline model of the isogeometric analysis performed on M and the topological optimization results performed on the optimized model M1 when the target volume fraction is 50%. DETAILED DESCRIPTION

[0067] The present invention will be further described below in conjunction with the accompanying drawings.

[0068] like Figure 1 As shown in the figure, the bidirectional progressive isogeometric topology optimization method based on volume subdivision has the following specific steps:

[0069] Step 1: For the input original hexahedral mesh solid model M, the present invention uses the solid model M as the analysis model, as shown in FIG2(a). In order to initialize the design domain, the present invention uses the density ρ of all hexahedral units in the original hexahedral mesh solid model M e are all set to 1. An external load F = 1N is applied to the center of the right side of the model M, as shown in Figure 2(b).

[0070] Specifically, the density ρ e A value of 1 indicates that the unit has solid material, density ρ e is ρ min Indicates that there is no solid material in the element and it will be deleted after the iteration ends.

[0071] Step 2: In order to ensure that the subdivision result does not deviate too much from the initial model, the present invention uses the Catmull-Clark interpolation subdivision method to perform a volume subdivision on the original hexahedral mesh solid model M to obtain a solid model M1 as the optimization model. In order to construct the optimization model design domain, the density of the subdivided unit after volume subdivision inherits the density of the hexahedral unit of the previous level. And initialize the volume fraction of the optimization model M1 before entering the iteration.

[0072] Specifically, when performing volume subdivision, a hexahedral mesh unit is subdivided into 8 hexahedral mesh subunits. In step 2, the subdivided solid model M1 corresponding to M is generated using the following interpolation subdivision rule:

[0073] Step 2.1: First, assume that there is an internal point C in a model, and its neighborhood points can be represented by three types of points: D, K, and B. Point D is a point directly connected to C, point K is a point not directly connected to C but on the same surface, and point B is a point not directly connected to C and not on the same surface but on the same body. Now follow the steps below to give the virtual point calculation formulas for B, K, D, and C, and record their virtual points as V B 、VK 、V D 、V C .

[0074] (1)V B Keeping point B unchanged, we have V B =B.

[0075] (2) V K When V K When replacing point K, as shown in Figure 3(a), let C, V K The corresponding point set in the two adjacent bodies whose faces are common faces is {B1,V K ,B2}, when the point set is subdivided by three uniform B-spline, V K The corresponding cubic uniform B-spline subdivision limit points coincide with the K point.

[0076] (3)V D When V D When point D is replaced, as shown in Figure 3(b), with D and V D In the area corresponding to the four hexahedrons whose edges are common edges, V D When performing Catmull-Clark tessellation for the center of the surface, V D The corresponding limit point coincides with point D. According to the limit point formula of Catmull-Clark surface subdivision, V can be obtained D :

[0077]

[0078] Among them, v 1 Represents the coordinates of vertex v after a surface subdivision. Similarly, f j 1 They represent the coordinates of the j-th edge point and face point of point v after one subdivision, and n represents the number of edges connected to vertex v.

[0079] (4)V C When V C Instead of point C, at point V C When performing Catmull-Clark body subdivision, V C The corresponding limit point coincides with point C. According to the limit point formula of Catmull-Clark curved body subdivision, V can be obtained C :

[0080]

[0081] Among them, v 1 Represents the coordinates of vertex v after one subdivision. Similarly, f j 1 , They represent the coordinates of the j-th edge point, surface point, and block point of point v after one subdivision, and n represents the number of edges connected to vertex v.

[0082] Step 2.2: After completing the calculation of the virtual points of points B, K, D, and C, as shown in Figure 3(c), the virtual points corresponding to each point in the adjacent domain point set of point C are obtained. The Catmull-Clark volume subdivision rule can be used to give a rule for calculating the contribution value of each virtual point to the volume point, surface point, and edge point added to the corresponding area.

[0083] (1) For each neighboring body of point C, as shown in Figure 3(d), let the set of virtual points contained in the current neighboring body be {p c1 ,p c2 ,…,p c8}, the contribution of the virtual point to the hexahedron point is calculated as follows:

[0084]

[0085] (2) For each adjacent face of point C, when the cubic uniform B-spline subdivision algorithm is applied to each point in the face, the contribution value is the average value of the limit points corresponding to these points. As shown in Figure 3(e), let P be a virtual point on the current adjacent face. When the cubic uniform B-spline algorithm is applied to point P, the corresponding limit point is in It is a virtual point corresponding to the point that is not directly connected to C and is not on the same surface but belongs to the same entity.

[0086] At this time, let the virtual point set of the current adjacent face be {p f1 ,p f2 ,p f3 ,p f4}, and remember is a virtual point p fi The corresponding limit point, the contribution of the virtual point set to the adjacent surface point is calculated as follows:

[0087]

[0088] (3) For each adjacent edge of point C, the contribution value is calculated as the average value of the corresponding limit points of the two points when the Catmull-Clark surface subdivision is applied to the two endpoints of the edge. As shown in Figure 3(f), let p e1 ,p e2 is a virtual point on the current adjacent edge, and When Catmull-Clark tessellation is applied to point, p ei When the virtual point is the limit point, the contribution of the virtual point to the adjacent edge point is calculated as follows:

[0089]

[0090] Step 2.3, after completing the above process with the point inside each body as the center, the contribution of each body, surface, and edge point is summed and averaged to obtain the newly inserted body point, surface point, and edge point. Then, according to the topological rules of Catmull-Clark body subdivision, the topological connection is performed to obtain the result model after subdivision once.

[0091] Step 3: In order to improve computational efficiency and save memory space, the present invention constructs a three-variable Bézier spline body as an approximation of the limit subdivision body of M according to the body subdivision limit point formula, and applies the load of the model M to its corresponding Bézier spline body.

[0092] Specifically, the volume subdivision limit point formula used in step 3 is:

[0093]

[0094] Among them, v 1 Represents the coordinates of vertex v after one subdivision. Similarly, f j 1 , They represent the coordinates of the jth edge point, face point, and block point of point v after one subdivision, and n represents the number of edges connected to vertex v. Using the above limit point formula, the generated Bézier body collective interpolation hexahedral mesh model M can be directly subdivided to generate the Catmull-Clark subdivision limit body.

[0095] Step 4: Construct the global stiffness matrix K and right-hand term F of the Bézier spline body corresponding to the model M according to the applied load; use the global stiffness matrix K and right-hand term F to perform isogeometric analysis to solve the global displacement matrix U; construct the initial stiffness matrix K of the Bézier spline body corresponding to each subdivided unit in the optimization model M1 i,0 , which is used for subsequent unit sensitivity calculations.

[0096] Specifically, in step 4, K i,0 The calculation method is:

[0097]

[0098] Where B is the strain matrix, D is the elastic matrix, and the matrix B = (B1, B1, ..., B m), m is the number of basis functions of the Bézier spline, x, y and z represent the three spatial directions of the hexahedral unit, |J i | represents the Jacobian determinant of the hexahedral element corresponding to the Bézier spline.

[0099] Step 5: Determine the volume fraction V of the current iterative optimization model M1 i Has the target volume fraction V been met? * , if the target volume fraction V is not met * Then go to step 6. If the target volume fraction V is met * , then further determine whether the convergence formula is satisfied. If not, execute step 6. If the convergence formula is satisfied, end the iteration and execute step 9.

[0100] Specifically, the traditional BESO method ends after satisfying the volume constraint, while the structure may still have room for optimization. The present invention chooses to introduce a convergence judgment formula to ensure that the iteration is terminated when the structural flexibility value of 10 consecutive iterations tends to be stable while satisfying the volume constraint.

[0101] The convergence formula in step 5 is as follows:

[0102]

[0103] Among them, C k-i+1 is the flexibility value at the k-i+1th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of previous and next iterations involved in judging convergence, which are set to 0.1% and 5 respectively.

[0104] Step 6: If the volume fraction V of the current iteration i and the target volume fraction V * When the volume fraction is equal to V, in order to stabilize the optimization result, the present invention selects a fixed subsequent iteration volume fraction V * Otherwise, update the volume fraction V for the next iteration i+1 .

[0105] Specifically, the volume fraction update strategy in step 6 is as follows:

[0106] Step 6.1, calculate the current iterative evolution ratio ER. In order to avoid the situation where too many units are deleted when the number of iterations is large enough to cause numerical instability, the present invention uses an adaptive function to update the iterative evolution ratio each time. Take the target volume fraction V * =0.5; take the maximum evolution ratio ER max =0.02; take the minimum evolution ratio ER min =0.01:

[0107]

[0108] Step 6.2: Calculate the volume fraction for the next iteration based on the current evolution ratio:

[0109] V i+1 =V i (1-ER)

[0110] Step 7: Use the initial stiffness matrix K obtained in step 4 i,0 The sensitivity of each unit of the optimization model M1 is calculated with the displacement matrix U, and each unit sensitivity is filtered. The filtered unit sensitivity is then averaged in combination with the unit historical sensitivity to obtain the final unit sensitivity.

[0111] Specifically, the common "checkerboard" phenomenon in topology optimization refers to the periodic distribution of the presence or absence of units in the topological area of ​​the structure, which can be solved by the unit sensitivity filtering method. In addition, the traditional BESO method is difficult to converge. After multiple calculations and verifications, the convergence problem can be solved well by averaging the current iteration sensitivity using the historical iteration sensitivity.

[0112] The specific calculation process of sensitivity in step 7 is:

[0113] Step 7.1, use the following formula to calculate the initial unit sensitivity of the current iteration:

[0114]

[0115] Among them, u i To optimize the displacement of the i-th unit in the model M1 in the displacement matrix U, ρ i is the unit density, K i,0 is the initial stiffness matrix of the unit, p is the penalty factor, usually p=2.

[0116] Step 7.2: Use the sensitivity filter function to filter each unit:

[0117]

[0118] Where M is the filter radius r min The number of all units in the range, α j0 is the unit sensitivity within the filtering range, r ij is the distance from the center of unit i to the center of unit j, w(r ij ) is the weight factor, usually r min -r ij .

[0119] Step 7.3: Combine the unit sensitivity of the previous iteration and average the filtered unit sensitivity: Where k represents the current iteration number.

[0120] Step 8: Use the current volume fraction V i , the final target volume fraction V * The density of each unit in the optimized model M1 is updated based on the sensitivity of each unit in the optimized model M1, and the number of iterations is increased by 1 and the process returns to step 4.

[0121] Specifically, the cell density update strategy in step 8 is as follows:

[0122] Step 8.1: Sort all unit sensitivities in descending order, denoted by α max is the maximum unit sensitivity of the current iteration, α min The minimum cell sensitivity of the current iteration is set to ρ, and the minimum cell density is set to min , indicating an empty cell.

[0123] Step 8.2: Combine α max , α min and volume fraction V i , using a binary strategy to find the sensitivity threshold α th .

[0124] Step 8.3, Sensitivity α i ≤a th The cell density is updated from 1 to ρ min , sensitivity α i >α th The cell density is changed from ρ min Updated to 1. min =0.0001,ρ max =1.0.

[0125] Step 9: Delete all the units with the smallest density in the optimized model M1, and then construct a Bézier body for each remaining unit to obtain the final topology optimization model of the solid model.

[0126] The examples are as follows:

[0127] FIG4(a) and FIG4(b) respectively correspond to the mesh model and spline model of the isogeometric analysis and topology optimization results on M when the target volume fraction is 50%.

[0128] Figure 4(c) and Figure 4(d) correspond to the mesh model and spline model of the isogeometric analysis performed on M and the topological optimization results performed on the optimized model M1 when the target volume fraction is 50%, respectively.

[0129] The present invention can generate an interpolation spline body for an input hexahedral model and perform isogeometric topology optimization on the spline body, thereby generating an optimization result that satisfies clear structural boundaries and does not have intermediate density units, and the method can perform optimization faster under the same input.

Claims

1. A bidirectional progressive isogeometric topology optimization method based on volume subdivision, characterized in that: The following steps are involved: Step 1, initialize the density of all hexahedral units in the input hexahedral mesh solid model M, perform k1 volume subdivisions on M to obtain the optimized model M1, and initialize the volume fraction of M1 before entering the iteration; Step 2: construct a three-variable Bézier spline body according to the body subdivision limit point formula, and apply the external load on the model M to its corresponding Bézier spline body; Step 3: According to the global stiffness matrix and right-hand term of the Bézier spline body after the load is applied, use isogeometric analysis to solve the displacement matrix U; construct the initial stiffness matrix K of the Bézier spline body corresponding to each subdivided unit in M1 i,0 ; Step 4: Determine the volume fraction V of the current iterative optimization model M1 i Whether the target volume fraction V is met * , if it does not meet the requirements, execute step 5, otherwise determine whether the convergence formula is met; If the convergence formula is not satisfied, execute step 5, otherwise terminate the iteration and obtain the final topology optimization model of the solid model; Step 5: If the target volume V of the current iteration i With V * When they are equal, the target volume fraction of subsequent iterations is fixed to V * ; Otherwise update the volume fraction V for the next iteration i+1 ; Step 6: Based on the initial stiffness matrix K i,0 and displacement matrix U, calculate the unit sensitivity of the optimization model M1, and update the density of each unit of the optimization model M1. The number of iterations is increased by 1 and returns to step 4.

2. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 1, characterized in that: In the step 1, after each volume subdivision, the density of the subdivided unit inherits the density of the hexahedral unit at the previous level.

3. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 2, characterized in that: The specific implementation process of obtaining the optimization model M1 in step 2 is as follows: Step 2.1, suppose there is an internal point C in a model, and its neighboring points are represented by three types of points: D, K and B; among them, point D is a point directly connected to C, point K is a point not directly connected to C but on the same surface, and point B is a point not directly connected to C and not on the same surface but on the same body; give the virtual point calculation formulas of B, K, D, and C respectively, and record their virtual points as V B 、V K 、V D 、V C ; Step 2.2, after completing the calculation of the virtual points of the above points B, K, D, and C, the virtual points corresponding to each point in the point set of the adjacent domain of point C are obtained, and then the contribution values ​​of the body points, surface points, and edge points added to the area corresponding to each virtual point are calculated using the Catmull-Clark volume subdivision rule; Step 2.3: After completing the above process with the internal point of each body as the center, the contribution of each body, surface, and edge point is summed and averaged to obtain the newly inserted body point, surface point, and edge point; then, topological connection is performed according to the topological rules of Catmull-Clark body subdivision to obtain the optimized model after one subdivision.

4. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 3, characterized in that: In step 2.1, the virtual point calculation formulas of B, K, D, and C are given respectively, and their virtual points are recorded as V B 、V K 、V D 、V C , the specific process is as follows: V B Keeping point B unchanged, we have V B =B; V K When V K When replacing point K, C, V K The corresponding construction point set in the two adjacent bodies whose faces are common faces is subdivided by three uniform B-spline subdivisions. K The corresponding cubic uniform B-spline subdivision limit point coincides with the K point; V D When V D When replacing point D, D In the area corresponding to the four hexahedrons whose edges are common edges, V D When performing Catmull-Clark tessellation for the center of the surface, V D The corresponding limit point coincides with point D; V C When V C Instead of point C, at point V C When performing Catmull-Clark body subdivision, V C The corresponding limit point coincides with point C.

5. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 4, characterized in that: The contribution value in step 2.2 is calculated as follows: For each neighbor of point C, let the set of virtual points contained in the current neighbor be {p c1 ,p c2 ,…,p c8 }, then the contribution of the virtual point to the point of the hexahedron is the average value of the virtual points in the set; For each adjacent face of point C, when the cubic uniform B-spline subdivision algorithm is applied to each point in the face, the contribution value is the average value of the limit points corresponding to these points; For each adjacent edge of point C, the contribution value is calculated as the average value of the corresponding limit points of the two points when Catmull-Clark surface subdivision is applied to the two endpoints of the edge.

6. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 5, characterized in that: The convergence formula is as follows: Among them, C k-i+1 is the flexibility value at the k-i+1th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of previous and next iterations involved in judging convergence.

7. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 6, characterized in that: In step 4, the iteration is ended to obtain the final topology optimization model of the solid model. The specific process is: delete all units with the smallest density in the optimization model M1, construct a Bézier spline body for each remaining unit, and obtain the final topology optimization model of the solid model.

8. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 7, characterized in that: In step 5, the volume fraction V of the next iteration is updated i+1 The details are as follows: Step 6.1, calculate the current iteration evolution ratio ER: Among them, ER max , ER min are the maximum and minimum evolutionary ratios, respectively; Step 6.2: Calculate the volume fraction for the next iteration based on the current evolution ratio: V i+1 =V i (1-ER).

9. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 8, characterized in that: The specific calculation process of unit sensitivity in step 6 is: The initial unit sensitivity of the current iteration is calculated using the following formula: Among them, u i To optimize the displacement of the i-th unit in the model M1 in the displacement matrix U, ρ i is the unit density, K i,0 is the initial stiffness matrix of the unit, p is the penalty factor; Use the sensitivity filter function to filter each unit: Where M is the filter radius r min The number of all units in the range, α j0 is the unit sensitivity within the filtering range, r ij is the distance from the center of unit i to the center of unit j, w(r ij ) is the weight factor, and its value is r min -r ij ; Combined with the unit sensitivity of the previous iteration, the filtered unit sensitivity is averaged to obtain the final unit sensitivity: Where k represents the current iteration number.

10. The bidirectional progressive isogeometric topology optimization method based on volume subdivision according to claim 9, characterized in that: The specific process of updating the density of each unit of the optimization model M1 in step 6 is as follows: Sort all unit sensitivities in descending order, denoted by α max , α min They are the maximum and minimum unit sensitivities of the current iteration respectively; Set the minimum cell density to ρ min , indicating an empty cell; Combined with α max , α min and volume fraction V i , using a binary strategy to find the sensitivity threshold α th ; Sensitivity α i ≤α th The cell density is updated from 1 to ρ min , sensitivity α i >α th The cell density is changed from ρ min Updated to 1.

Images