A bidirectional asymptotic isogeometric topology optimization method based on volume subdivision

Abstract

This invention discloses a bidirectional asymptotic isogeometric topology optimization method based on volume subdivision. The method first initializes the density of all hexahedral elements in the input hexahedral mesh solid model M, and then subdivides M by k1 times to obtain the optimized model M1. Next, a three-variable Bézier spline is constructed, and external loads on M are applied to their corresponding splines. The displacement matrix U is solved using isogeometric analysis. The initial stiffness matrix of the Bézier spline corresponding to each subdivided element in M1 is then constructed. Finally, a judgment is made based on the target volume fraction and the convergence formula. If the target volume fraction and the convergence formula are satisfied, the iteration ends, and the topology optimization model is obtained; otherwise, the volume fraction for the next iteration is updated, the element sensitivity of optimized M1 is calculated, and the density of each element in M1 is updated, and the iteration is repeated. This invention efficiently performs isogeometric topology optimization and obtains topology optimization results with clearly defined spline representations.

Description

Technical Field

[0001] This invention belongs to the field of isogeometric analysis and structural optimization applications, and relates to a bidirectional asymptotic topology optimization method based on variable deletion rate, specifically a bidirectional asymptotic isogeometric topology optimization method based on volume subdivision. Background Technology

[0002] Isogeometric analysis is a novel numerical computation method that unifies design and analysis. It can more accurately represent the geometric features of a model, achieving higher precision. This enables seamless integration of CAD / CAE, while simultaneously improving computational accuracy and efficiency.

[0003] Structural optimization refers to the use of various methods and techniques in the design and analysis of engineering structures to achieve optimal shape, material distribution, and load-bearing capacity while meeting performance requirements. Its aim is to improve structural efficiency, reduce material costs, and simultaneously ensure structural safety and stability. Traditional structural optimization results require post-processing before being imported into CAD systems; therefore, seamlessly integrating design and analysis, as well as handling complex constraints, presents significant challenges for structural optimization.

[0004] However, existing approximation-type volume subdivision algorithms cause the model shape to deviate further and further from the initial mesh after multiple iterations. Furthermore, the SIMP topology optimization method suffers from unclear boundaries and intermediate density cells, failing to meet the requirement of clear structural boundaries for volume subdivision.

[0005] Therefore, this invention presents a bidirectional progressive topology optimization method based on interpolation-type volume subdivision, which efficiently obtains topology optimization results with clear boundaries that can be directly imported into CAD systems using the input hexahedral model. Summary of the Invention

[0006] The purpose of this invention is to achieve the integration of geometric representation and simulation analysis, and to improve the efficiency of the computational method and the accuracy of the resulting model. It proposes a variable deletion rate bidirectional asymptotic topology optimization method based on interpolation volume subdivision. This addresses the problem of unclear structure and the presence of intermediate density units in the SIMP topology optimization method. Furthermore, by utilizing the interpolation volume subdivision method, a high-quality topology optimization model that can be directly imported into a CAD system can be generated.

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

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

[0009] Step 2: Subdivide M into optimized model M1 by k1 times, k1≥1. After each subdivision, the density of the subdivided unit inherits the density of the hexahedral unit of the previous level, and initialize the volume fraction of optimized model M1 before entering the iteration.

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

[0011] Step 4: Based on the global stiffness matrix and right-hand side terms of the Bézier spline after the applied load, perform isogeometric analysis to solve for the displacement matrix U. Construct the initial stiffness matrix K of the Bézier spline corresponding to each subdivided element 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 proceed to step 6. If the target volume fraction V is satisfied... * If the convergence formula is not satisfied, proceed to step 6. If the convergence formula is satisfied, end the iteration and proceed to step 9.

[0013] Step 6: If the target volume V of the current iteration i With V * When they are equal, fix the target volume for subsequent iterations as 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 The sensitivity of each element in the optimization model M1 is calculated using the displacement matrix U.

[0015] Step 8: Utilize the current volume fraction V i Target volume fraction V * The sensitivity of each unit in the optimized model M1 is updated, and the density of each unit in the optimized model M1 is updated. The iteration count is incremented by 1, and the process returns to step 4.

[0016] Step 9: Delete all the elements with the lowest density in the optimization model M1, and then construct a Bézier spline for each of the remaining elements 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 rules:

[0018] Step 2.1: First, assume an internal point C in a model. Its neighborhood points can be represented as three types of points: D, K, and B. Point D is directly connected to C; point K is not directly connected to C but belongs to the same plane; and point B is not directly connected to C and does not belong to the same plane but belongs to the same volume. Now, the formulas for calculating virtual points B, K, D, and C are given below, and their virtual points are denoted 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 Satisfying V K When replacing point K, use C and V K Construct a point set from corresponding points in two adjacent bodies whose faces are common faces. When this point set is subdivided into cubic uniform B-splines, V K The corresponding cubic uniform B-spline subdivision limit point coincides with point K.

[0021] (3)V D Satisfying V D When replacing point D, in the case of D and V D In the region corresponding to the four hexahedrons whose edges share a common edge, V D When performing Catmull-Clark subdivision of the surface around the center, V D The corresponding limit point coincides with point D.

[0022] (4)V C Satisfying V C Replace point C with point V C When performing Catmull-Clark volume subdivision, V C The corresponding limit point coincides with point C.

[0023] Step 2.2: After completing the virtual point calculations for points B, K, D, and C, we now have the virtual points corresponding to each point in the adjacent point set of point C. We can then use the Catmull-Clark volume subdivision rules to calculate the contribution values ​​of the added volume points, face points, and edge points to the region corresponding to each virtual point:

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

[0025] (2) For each adjacent surface of point C, when the cubic uniform B-spline subdivision algorithm is applied to each point in the surface, 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 follows: when the Catmull-Clark surface subdivision acts on the two endpoints of the edge, the contribution value is the average of the corresponding limit points of the two points.

[0027] Step 2.3: After completing the above process using points inside each body as centers, the contributions of each body, face, and edge point are summed and averaged to obtain the final newly inserted body points, face points, and edge points. Then, by performing topological connections according to the Catmull-Clark volume subdivision topology rules, the optimized model after one subdivision is obtained.

[0028] As a preferred option, the formula for the volume subdivision limit point in step 3 is as follows:

[0029]

[0030] Among them, v 1 This represents the coordinates of vertex v after one subdivision. Similarly, f j 1 , Let represent the coordinates of the j-th edge point, face point, and block point of point v after one subdivision, respectively, and n represent the number of edges connected to vertex v. Using the limit point formula above, the generated Bézier volume's set interpolated hexahedral mesh model M can be directly subdivided to generate the Catmull-Clark subdivision limit volume.

[0031] As a preferred option, the formula for determining convergence in step 5 is as follows:

[0032]

[0033] Among them, C k-i+1 τ is the compliance value at the (k-i+1)th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of iterations before and after participating in the convergence judgment.

[0034] As a preferred option, the volume fraction update strategy in step 6 is as follows:

[0035] Step 6.1: Calculate the evolution ratio ER for the current iteration:

[0036]

[0037] Among them, ER max ER minThese represent 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 option, the specific calculation process for sensitivity in step 7 is as follows:

[0041] Step 7.1: Calculate the initial unit sensitivity for the current iteration using the following formula:

[0042]

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

[0044] Step 7.2: Filter each unit using a sensitivity filtering function:

[0045]

[0046] Where M is the filtration radius r min The number of all units within the range, α j0 For the unit sensitivity within the filtering range, r ij Let w(r) be the distance from the center of element i to the center of element j. ij ) is the weighting factor, with a value of r. min -r ij .

[0047] Step 7.3: Combining the cell sensitivity from the previous iteration, average the filtered cell sensitivity to obtain the final cell sensitivity: Where k represents the current iteration number.

[0048] As a preferred option, the cell density update strategy in step 8 is as follows:

[0049] Step 8.1: Sort all unit sensitivities in descending order, and denote α. max α represents the maximum unit sensitivity in the current iteration. min To determine the minimum element sensitivity for the current iteration, the minimum element density is set to ρ. min , indicates an empty cell.

[0050] Step 8.2, combining α max α min and volume fraction V iA binary search strategy is used 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 unit density from ρ min Updated to 1.

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

[0053] The above scheme combines a bidirectional progressive topology optimization method based on a "soft-kill" approach with an isogeometric topology optimization framework based on interpolation volume subdivision and unstructured three-variable spline construction. This enables efficient isogeometric topology optimization of input hexahedral meshes, yielding topology optimization results with clear boundaries that can be directly imported into CAD systems. This improves the quality and efficiency of existing optimization results within the multi-resolution isogeometric topology optimization framework. Attached Figure Description

[0054] Figure 1 This is a flowchart of the present invention;

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

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

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

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

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

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

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

[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 results of isogeometric analysis and topology optimization on M when the target volume fraction is 50%;

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

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

[0066] Figure 4(d) shows the spline model of the results of isogeometric analysis on M and topology optimization on the optimization model M1 when the target volume fraction is 50%. Detailed Implementation

[0067] The invention will now be further described with reference to the accompanying drawings.

[0068] like Figure 1 As shown, the specific steps of the bidirectional asymptotic isogeometric topology optimization method based on volume subdivision are as follows:

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

[0070] Specifically, density ρ e A value of 1 indicates that the unit contains solid material, and its density is... ρ e is ρ min This indicates that the cell does not contain any physical material and will be deleted after the iteration ends.

[0071] Step 2: To ensure that the subdivision result does not deviate too much from the initial model, this invention uses the Catmull-Clark interpolation subdivision method to perform a first-order volume subdivision on the original hexahedral mesh solid model M to obtain solid model M1 as the optimized model. To construct the design domain of the optimized model, the density of the subdivided elements inherits the density of the hexahedral elements of the previous level. The volume fraction of the optimized model M1 is then initialized before entering the iteration.

[0072] Specifically, during volume subdivision, one hexahedral mesh element is subdivided into eight hexahedral mesh sub-elements. In step 2, the subdivided solid model M1 corresponding to M is generated using the following interpolation subdivision rules:

[0073] Step 2.1: First, assume an internal point C in a model. Its neighborhood points can be represented as three types of points: D, K, and B. Point D is directly connected to C; point K is not directly connected to C but belongs to the same plane; and point B is not directly connected to C and does not belong to the same plane but belongs to the same volume. Now, the formulas for calculating virtual points B, K, D, and C are given below, and their virtual points are denoted 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 Satisfying V K When replacing point K, as shown in Figure 3(a), let C and V be the points. K The corresponding set of points in two adjacent bodies whose faces are common faces is {B1, V}. K When performing cubic uniform B-spline subdivisions on this point set, V K The corresponding cubic uniform B-spline subdivision limit point coincides with point K. Specifically, there are...

[0076] (3)V D Satisfying V D When replacing point D, as shown in Figure 3(b), with D and V... D In the region corresponding to the four hexahedrons whose edges share a common edge, V D When performing Catmull-Clark subdivision of the surface around the center, V D The corresponding limit point coincides with point D. Therefore, according to the limit point formula for Catmull-Clark surface subdivision, V can be obtained. D :

[0077]

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

[0079] (4)V C Satisfying V C Replace point C with point V C When performing Catmull-Clark volume subdivision, V C The corresponding limit point coincides with point C. Therefore, according to the limit point formula of Catmull-Clark body subdivision, V can be obtained. C :

[0080]

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

[0082] Step 2.2: After completing the virtual point calculations for points B, K, D, and C, as shown in Figure 3(c), we now have the virtual points corresponding to each point in the adjacency set of point C. We can then use the Catmull-Clark volume subdivision rules to calculate the contribution values ​​of each virtual point to the volume points, face points, and edge points added to the corresponding region.

[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 solid point of the hexahedron is calculated as follows:

[0084]

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

[0086] Let the set of virtual points of the current adjacent face be {p}. f1 ,p f2 ,p f3 ,p f4}, and record For virtual point p fi For the corresponding limit point, the contribution of this virtual point set to the adjacent face point is calculated as follows:

[0087]

[0088] (3) For each adjacent edge of point C, the contribution value is calculated as the average of the corresponding limit points of the two points when the Catmull-Clark surface subdivision acts on the two endpoints of the edge. As shown in Figure 3(f), let p e1 ,p e2 Let be a virtual point on the current adjacent edge, and record . When the Catmull-Clark surface subdivision is applied at point p ei If the limit point is reached, 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 using points inside each body as centers, the contributions of each body, face, and edge point are summed and averaged to obtain the final newly inserted body points, face points, and edge points. Then, topological connections are performed according to the Catmull-Clark volume subdivision topology rules to obtain the subdivided model.

[0091] Step 3: To improve computational efficiency and save memory space, this invention constructs a three-variable Bézier spline body as an approximation of the limit subdivision body of M based on the volume subdivision limit point formula. The load of model M is then applied to its corresponding Bézier spline body.

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

[0093]

[0094] Among them, v 1 This represents the coordinates of vertex v after one subdivision. Similarly, f j 1 , Let represent the coordinates of the j-th edge point, face point, and block point of point v after one subdivision, respectively, and n represent the number of edges connected to vertex v. Using the limit point formula above, the generated Bézier volume's set interpolated hexahedral mesh model M can be directly subdivided to generate the Catmull-Clark subdivision limit volume.

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

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

[0097]

[0098] Where B is the strain matrix and D is the elasticity matrix, 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 element, |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 proceed to step 6. If the target volume fraction V is satisfied... * If the convergence formula is not satisfied, proceed to step 6. If the convergence formula is satisfied, end the iteration and proceed to step 9.

[0100] Specifically, the traditional BESO method terminates after the volume constraint is satisfied, but the structure may still have room for optimization. This invention introduces a convergence criterion formula to ensure that the iteration terminates only when the structural compliance value tends to stabilize after 10 consecutive iterations 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 compliance value at the (k-i+1)th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of iterations before and after participating in the convergence judgment, which are set to 0.1% and 5 respectively.

[0104] Step 6: If the volume fraction V of the current iteration i With target volume fraction V * When the values ​​are equal, in order to stabilize the optimization results, this invention chooses to fix the subsequent iteration volume fraction as 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 evolution ratio ER for the current iteration. To avoid numerical instability caused by deleting too many units when the number of iterations is large enough, this invention uses an adaptive function to update the evolution ratio for each iteration. The target volume fraction V is taken. * =0.5; take the maximum evolutionary ratio ER max =0.02; take the minimum evolutionary 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 element in the optimized model M1 is calculated using the displacement matrix U, and then filtered. The filtered element sensitivities are then averaged using the historical element sensitivities to obtain the final element sensitivities.

[0111] Specifically, the "chessboard" phenomenon is common in topology optimization, referring to the periodic distribution of the presence or absence of elements in the structural topology region. This can be addressed using element sensitivity filtering methods. Furthermore, the traditional BESO method struggles to converge. Multiple calculations have verified that averaging the current iteration sensitivity using historical iteration sensitivity effectively solves the convergence problem.

[0112] The specific calculation process for sensitivity in step 7 is as follows:

[0113] Step 7.1: Calculate the initial unit sensitivity for the current iteration using the following formula:

[0114]

[0115] Among them, u i To optimize the displacement of the i-th element in model M1 corresponding to the displacement matrix U, ρ i For the unit density, K i,0 Let p be the initial stiffness matrix of the element, and p be the penalty factor, which is usually taken as p = 2.

[0116] Step 7.2: Filter each unit using a sensitivity filtering function:

[0117]

[0118] Where M is the filtration radius r min The number of all units within the range, α j0 For the unit sensitivity within the filtering range, r ij Let w(r) be the distance from the center of element i to the center of element j. ij ) is the weighting factor, which is usually taken as r. min -r ij .

[0119] Step 7.3: Averaging the filtered unit sensitivities based on the unit sensitivities from the previous iteration: Where k represents the current iteration number.

[0120] Step 8: Utilize the current volume fraction V i Final target volume fraction V * The sensitivity of each unit in the optimized model M1 is updated, and the density of each unit in the optimized model M1 is updated. The iteration count is incremented 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, and denote α. max α represents the maximum unit sensitivity in the current iteration. min To determine the minimum element sensitivity for the current iteration, the minimum element density is set to ρ. min , indicates an empty cell.

[0123] Step 8.2, combining α max α min and volume fraction V i A binary search strategy is used 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 unit density from ρ min Updated to 1. Where ρ min =0.0001, ρ max =1.0.

[0125] Step 9: Delete all the elements with the lowest density in the optimization model M1, and then construct a Bézier volume for each remaining element to obtain the final topology optimization model of the solid model.

[0126] The following is an example:

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

[0128] Figures 4(c) and 4(d) show the mesh model and spline model of the topology optimization results performed on the optimization model M1 when the target volume fraction is 50%, respectively, after performing isogeometric analysis on M.

[0129] This invention can generate interpolated spline bodies from an input hexahedral model and perform isogeometric topology optimization on the spline bodies. It generates optimized results that satisfy clear structural boundaries and lack intermediate density elements, and under the same input conditions, this method can perform optimization faster.

Claims

1. A bidirectional asymptotic isogeometric topology optimization method based on volume subdivision, characterized in that, Includes the following steps: Step 1: Initialize the density of all hexahedral elements in the input hexahedral mesh solid model M, perform k1 volume subdivision 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 limit point formula of volume subdivision, and apply the external load on model M to its corresponding Bézier spline body; Step 3: Based on the global stiffness matrix and right-hand side terms of the Bézier spline after the applied load, solve for the displacement matrix U using isogeometric analysis; construct the initial stiffness matrix K of the Bézier spline corresponding to each subdivided element in M1. i,0 ; Step 4: Determine the volume fraction V of the current iterative optimization model M1. i Does it meet the target volume fraction V? * If the condition is not met, proceed to step 5; otherwise, determine whether the convergence formula is met. If the convergence formula is not satisfied, proceed to step 5; otherwise, end the iteration and obtain the final topology optimization model of the entity model. Step 5: If the target volume V of the current iteration i With V * When they are equal, fix the target volume fraction for subsequent iterations as V. * Otherwise, update the volume fraction V for the next iteration. i+1 ; Step 6: Based on the initial stiffness matrix K i,0 Given the displacement matrix U, calculate the element sensitivity of the optimized model M1, update the element density of the optimized model M1, increment the iteration count by 1, and return to step 4.

2. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 1, characterized in that, In step 1, the density of each subdivided unit after each volume subdivision inherits the density of its parent hexahedral unit.

3. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 2, characterized in that, The specific implementation process of the optimized model M1 obtained in step 2 is as follows: Step 2.1: Assume an interior point C in a model, and its neighborhood points are represented by three types: D, K, and B. Point D is directly connected to C, point K is not directly connected to C but belongs to the same plane, and point B is not directly connected to C and not on the same plane but belongs to the same volume. Give the formulas for calculating the virtual points of B, K, D, and C, and denote their virtual points as V. B V K V D V C ; Step 2.2: After completing the virtual point calculations for points B, K, D, and C, we obtain the virtual points corresponding to each point in the adjacent domain point set of point C. Then, using the Catmull-Clark volume subdivision rule, we give the contribution values ​​of the volume points, face points, and edge points added to the region corresponding to each virtual point. Step 2.3: After completing the above process with each point inside the body as the center, the contribution of each body, face, and edge point is summed and averaged to obtain the final newly inserted body point, face point, and edge point; then, topological connections are performed according to the topological rules of Catmull-Clark body subdivision to obtain the optimized model after one subdivision.

4. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 3, characterized in that, In step 2.1, the calculation formulas for the virtual points of B, K, D, and C are given respectively, and their virtual points are denoted 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 Satisfying V K When replacing point K, use C and V K When the corresponding construction point set in two adjacent volumes whose faces are common faces is subjected to cubic uniform B-spline subdivision, V K The corresponding cubic uniform B-spline subdivision limit point coincides with point K. V D Satisfying V D When replacing point D, in the case of D and V D In the region corresponding to the four hexahedrons whose edges share a common edge, V D When performing Catmull-Clark surface subdivision around the center of the surface, V D The corresponding limit point coincides with point D; V C Satisfying V C Replace point C with point V C When performing Catmull-Clark volume subdivision, V C The corresponding limit point coincides with point C.

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

6. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 5, characterized in that, The specific convergence formula is as follows: Among them, C k-i+1 τ is the compliance value at the (k-i+1)th iteration, k is the current iteration number, τ is the convergence tolerance, and N is the number of iterations before and after participating in the convergence judgment.

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

8. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 7, characterized in that, In step 5, the volume fraction V for the next iteration is updated. i+1 Specifically as follows: Step 6.1: Calculate the evolution ratio ER for the current iteration: Among them, ER max ER min These are the maximum and minimum evolution 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 asymptotic isogeometric topology optimization method based on volume subdivision according to claim 8, characterized in that, The specific calculation process for the unit sensitivity in step 6 is as follows: Calculate the initial unit sensitivity for the current iteration using the following formula: Among them, u i To optimize the displacement of the i-th element in model M1 corresponding to the displacement matrix U, ρ i For the unit density, K i,0 Here, p is the initial stiffness matrix of the element, and p is the penalty factor. Each unit is filtered using a sensitivity filtering function: Where M is the filtration radius r min The number of all units within the range, α j0 For the unit sensitivity within the filtering range, r ij Let w(r) be the distance from the center of element i to the center of element j. ij ) is the weighting factor, with a value of r. min -r ij ; Based on the cell sensitivity from the previous iteration, the filtered cell sensitivity is averaged to obtain the final cell sensitivity: Where k represents the current iteration number.

10. The bidirectional asymptotic isogeometric topology optimization method based on volume subdivision according to claim 9, characterized in that, The specific process for updating and optimizing the density of each element in model M1 in step 6 is as follows: Sort all unit sensitivities in descending order, and denote α. max α min These are the maximum and minimum unit sensitivities in the current iteration, respectively. The minimum unit density is set to ρ. min , indicates an empty cell; Combining α max α min and volume fraction V i A binary search strategy is used to find the sensitivity threshold α. th ; Sensitivity α i ≤α th The cell density is updated from 1 to ρ min Sensitivity α i >α th The unit density from ρ min Updated to 1.

Images