A topological optimization and size control method for beam structures based on geometric skeleton

Abstract

A method for size control of topological optimization of beam structures based on a geometric skeleton belongs to the technical field of topological optimization of beam structures. The present invention solves the problem that the size control results of existing methods for size control of beam structures based on a geometric skeleton cannot reach the optimal level. The present invention first simulates the heat conduction behavior through partial differential equations, and then solves the temperature distribution state, and then obtains the divergence distribution based on the temperature distribution, and finally identifies the geometric skeleton by the numerical value of the divergence. By expanding the geometric skeleton with a specified size radius, establishing a control set of required characteristic sizes, comparing the formed geometry with the original structure, and constructing characteristic size constraints including the minimum entity characteristic size, the maximum entity characteristic size and the minimum blank characteristic size, the effective control of the minimum entity characteristic size, the maximum entity characteristic size and the minimum blank characteristic size can be achieved. The method of the present invention can be applied to the size control of topological optimization of beam structures.

Description

Technical Field

[0001] The present invention belongs to the technical field of beam structure topology optimization, and in particular relates to a beam structure topology optimization size control method based on a geometric skeleton. Background Art

[0002] Topology optimization is a powerful tool for designing beam structures. However, the resulting structures are often overly complex and delicate, posing significant challenges to manufacturing and hindering the transition from design to actual product. To address this issue, two key approaches are currently under consideration: developing methods for controlling feature dimensions to avoid excessive branching in structural solutions; and developing manufacturing-compatible design methods to ensure that designs meet manufacturing requirements. Feature dimension control technology offers significant advantages in terms of versatility.

[0003] Among the many developed topology optimization methods, variable density methods have strong application flexibility and stability due to their continuous design variable schemes under pixelation or voxelization. In particular, the solid isotropic material penalty model (SIMP) in the variable density method has become one of the mainstream topology optimization methods. Dimension control for this method is usually carried out in the variable filtering and step mapping stages of the optimization. However, this approach often suffers from weak control capabilities, unstable methods, and implicit expressions. In addition to the above methods, dimension control methods based on geometric skeletons have also received widespread attention. However, geometric skeleton extraction is often based on graphics methods and does not consider the continuous design variable characteristics of topology optimization. This leads to a serious loss of sensitivity information, which significantly reduces the optimality of the dimension control results. Summary of the Invention

[0004] The purpose of the present invention is to solve the problem that the existing beam structure size control method based on geometric skeleton does not take into account the continuous design variable characteristics of topological optimization, resulting in the size control result cannot reach the optimal one, and proposes a beam structure topological optimization size control method based on geometric skeleton.

[0005] The technical solution adopted by the present invention to solve the above technical problems is: a method for topological optimization and size control of a beam structure based on a geometric skeleton, the method specifically comprising the following steps:

[0006] Step 1: Mesh the design domain of the beam structure to obtain a performance analysis mesh layer; expand the boundary of the design domain of the beam structure outward, mesh the expanded area, and obtain a thermal analysis mesh layer;

[0007] The density of each unit in the performance analysis grid layer is used as the design variable of each unit, and the number of iterations is initialized to p = 1, and the design variable value of each unit is initialized;

[0008] Step 2: After filtering the design variables of each unit, substitute the filtering results into the step function to obtain the physical density field;

[0009] Step 3: Calculate the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtain the sensitivity of the target value and volume fraction constraint value respectively;

[0010] Step 4: Obtain the correlation matrix L based on the performance analysis grid layer and the thermal analysis grid layer, and then obtain the temperature field distribution based on the correlation matrix L and the physical density field;

[0011] Step 5: Obtain the divergence distribution according to the temperature field distribution;

[0012] Step 6: extracting the solid skeleton and blank skeleton of the beam structure according to the divergence distribution obtained in step 5, and expanding the solid skeleton and the blank skeleton;

[0013] Constructing minimum entity feature size constraints and maximum entity feature size constraints according to the expanded entity skeleton, constructing minimum blank feature size constraints according to the expanded blank skeleton, and obtaining minimum entity feature size constraint value, maximum entity feature size constraint value, minimum blank feature size constraint value, and sensitivity of minimum entity feature size constraint value, maximum entity feature size constraint value, and minimum blank feature size constraint value;

[0014] Step 7: The target value, volume fraction constraint value, minimum physical feature size constraint value, maximum physical feature size constraint value, minimum blank feature size constraint value, and the sensitivity of the target value, the sensitivity of the volume fraction constraint value, the sensitivity of the minimum physical feature size constraint value, the sensitivity of the maximum physical feature size constraint value, and the sensitivity of the minimum blank feature size constraint value are brought into the optimization solver to obtain the updated design variable value of each unit;

[0015] Step 8: Calculate the target change rate based on the target value; and determine whether the iteration stop condition is met based on the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value, and the minimum blank feature size constraint value;

[0016] If the iteration stop condition is met, the iteration is stopped and the physical density field of the size-controlled rear beam structure is output;

[0017] If the iteration stop condition is not met, set p=p+1 and return to step 2 using the updated design variable values ​​of each unit.

[0018] Furthermore, the specific process of step 2 is as follows:

[0019] Filter the design variables of unit j to

[0020]

[0021] Where, is the filter value corresponding to the design variable of unit j, M j The circle is centered at unit j and r f is the set of all units within the radius, ρ i is the design variable of unit i, w jo is the weight corresponding to the design variable of unit i;

[0022] The vector is composed of the filter values ​​corresponding to the design variables of each unit Then Substituting the step function H in equation (2), we can obtain the physical density field

[0023]

[0024] Where β represents the mapping strength of the step function, and η represents the threshold.

[0025] Furthermore, the weight w corresponding to the design variable of unit i ji The calculation method is:

[0026] w ji =max(0,r f -dist(j,i)) (3)

[0027] where dist(j,i) is the distance between the center of cell j and the center of cell i.

[0028] Furthermore, in step 4, the temperature field distribution is obtained according to the correlation matrix L and the physical density field. The specific process is:

[0029] Set up the partial differential equation:

[0030]

[0031] Where, represents the gradient operator, t is the time of heat conduction, and q represents the temperature;

[0032] Simplifying formula (4) according to time, we get:

[0033]

[0034] Where q (0) represents the temperature at t = 0;

[0035] The temperature in equation (5) is solved by the finite element method:

[0036]

[0037] Where Q represents the heat conduction matrix; q′ is the vector of the temperature of each unit; T is the matrix composed of the node information after converting the unit information into node information; and L is the association matrix between the performance analysis grid layer and the thermal analysis grid layer.

[0038] Furthermore, the specific process of step 5 is as follows:

[0039] Step 51: Calculate the temperature q of each unit j j Gradient in the x direction and the gradient in the y direction

[0040]

[0041] Where, Represents the gradient operators in the x and y directions, N e’ is the e'th shape function of element j; q j,e’ is the temperature of the e'th node in element j;

[0042] Step 52: Normalize the gradient calculated in step 41 to:

[0043]

[0044] Where, The normalized value of Normalized value of the intermediate variable

[0045] The divergence Δq on unit j is calculated using formula (10) j :

[0046]

[0047] Step 53: Substitute the divergence Δq on unit j j Normalized to:

[0048]

[0049] Where, is the divergence Δq j The normalized value of , where e is the base of the natural logarithm.

[0050] Furthermore, in step 6, the solid skeleton and the blank skeleton of the beam structure are extracted according to the divergence distribution obtained in step 5, and the solid skeleton and the blank skeleton are expanded; the specific process is:

[0051] The vector composed of the normalized divergence values ​​on each unit is recorded as Will Substitute the step function H into equation (12) to generate the entity skeleton s solid :

[0052]

[0053] Where, β solid is the mapping intensity, η solid is the parameter of the step function;

[0054] The minimum expansion size of the entity skeleton set The kth element in for:

[0055]

[0056] Where γ is the control parameter, e is the base of the natural logarithm, and N k It is centered on unit k and r minsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in ;

[0057] The maximum inflated size of the entity skeleton set The kth element in for:

[0058]

[0059] Where γ is the control parameter, N k ' is centered on unit k and r maxsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in ;

[0060] Will Substitute the step function H into equation (15) to generate the blank skeleton s void :

[0061]

[0062] Where, β void is the mapping intensity, η void is the parameter of the step function;

[0063] The expanded skeleton set of the blank skeleton The kth element in for:

[0064]

[0065] Where γ is the control parameter, N k " is centered on unit k and r minvoid is the set of all cells within the circle of the expansion radius, s void(i) It is the skeleton void The i-th element in .

[0066] Furthermore, the minimum physical feature size constraint is:

[0067]

[0068] Where n represents the number of discrete units, C minsolid It is the minimum entity feature size constraint value, and the superscript T represents transpose.

[0069] Furthermore, the maximum physical feature size constraint is:

[0070]

[0071] Among them, C maxsolid is the maximum solid feature size constraint value.

[0072] Furthermore, the minimum blank feature size constraint is:

[0073]

[0074] Among them, C minvoid is the minimum blank feature size constraint value.

[0075] Furthermore, the target change rate is calculated based on the target value; and whether the iteration stop condition is satisfied is determined based on the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value, and the minimum blank feature size constraint value; the specific process is:

[0076] Step 81: According to the target value obj of the cantilever beam after the current iteration (p) Calculating target rate of change

[0077]

[0078] Where, obj (p-l) is the target value of the cantilever beam after the plth iteration, obj (p-l-1) is the target value of the cantilever beam after the pl-1th iteration, |obj (p) | represents obj (p) The absolute value of l = 0, 1, 2, 3, 4;

[0079] Step 82: Set the iteration stop condition to satisfy both condition (1) and condition (2);

[0080] Condition (1): The target change rate is less than 0.001 for at least five consecutive iterations;

[0081] Condition (2): The minimum entity feature size constraint value of the current iteration is less than the threshold ε minsolid , the maximum entity feature size constraint value is less than the threshold ε maxsolid , the minimum blank feature size constraint value is less than the threshold ε minvoid .

[0082] The beneficial effects of the present invention are:

[0083] The present invention extracts the geometric skeleton. First, the heat conduction behavior is simulated by partial differential equations, and then the temperature distribution state is solved. Then, the divergence distribution is obtained based on the temperature distribution. Finally, the geometric skeleton is identified by the numerical value of the divergence. The use of partial differential equations to simulate the heat conduction process can ensure that the skeleton is continuously differentiable throughout the entire generation process, thereby retaining all sensitivity information and achieving perfect compatibility with the classic gradient optimization solver. In order to achieve the regulation of characteristic size, the present invention expands the extracted skeleton with a specified size radius, establishes a control set of required characteristic sizes, compares the formed geometry with the original structure, and constructs characteristic size constraints including the minimum entity characteristic size, the maximum entity characteristic size, and the minimum blank characteristic size. It can achieve effective regulation of the minimum entity characteristic size, the maximum entity characteristic size, and the minimum blank characteristic size, thereby achieving optimal control of the characteristic sizes of the entity domain and the blank domain. The method of the present invention can be adapted to the currently available variable density topology optimization method and has good versatility. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 It is a skeleton-generated framework diagram based on PDE;

[0085] Figure 2 It is a schematic diagram of two-layer grid;

[0086] Among them, (a) is the temperature distribution of the thermal analysis grid layer, (b) is the divergence distribution of the thermal analysis grid layer, (c) is the temperature distribution of the performance analysis grid layer, and (d) is the divergence distribution of the performance analysis grid layer;

[0087] Figure 3 It is a schematic diagram of the minimum entity feature size constraint;

[0088] Among them, (a) is the original structural entity, (b) is the expanded skeleton of the minimum entity characteristic size, and (c) is the complement of the original structural entity;

[0089] Figure 4 It is a schematic diagram of the maximum entity feature size constraint;

[0090] Where (a) is the original structural entity, (b) is the complement of the expanded skeleton with the largest entity feature size, and (c) is the original structural entity area;

[0091] Figure 5 is a schematic diagram of the minimum blank feature size constraint;

[0092] Among them, (a) is the original structural entity, (b) is the expanded skeleton with the minimum blank feature size, and (c) is the original structural entity area;

[0093] Figure 6 are the design domain and boundary conditions of the cantilever beam;

[0094] Figure 7 It is the topological optimization result of the cantilever beam with minimum entity characteristic size constraint;

[0095] Among them, (a) is the optimized structure; (b) is the solid skeleton; (c) is the expanded skeleton with the minimum solid feature size; (d) is the comparison between the complement of the optimized structure and the expanded skeleton;

[0096] Figure 8(a) shows the convergence history of the target value of the cantilever beam under the constraint of the minimum entity feature size;

[0097] Figure 8(b) shows the volume fraction convergence history of the cantilever beam under the constraint of the minimum physical feature size;

[0098] Figure 9 is the topology optimization result constrained by the maximum entity characteristic size of the cantilever beam;

[0099] Among them, (a) is the optimized structure; (b) is the solid skeleton; (c) is the expanded skeleton with the maximum solid feature size; (d) is the comparison between the complement of the optimized structure and the maximum expanded skeleton;

[0100] Figure 10(a) shows the convergence history of the target value of the cantilever beam under the constraint of the maximum physical feature size;

[0101] Figure 10(b) Volume fraction convergence history of the cantilever beam under the maximum solid feature size constraint;

[0102] Figure 11 are the design domain and boundary conditions of the MBB beam;

[0103] Figure 12 is the MBB topology optimization result with minimum and maximum entity feature size constraints;

[0104] Figure 13 is the MBB topology optimization result with minimum, maximum solid, and minimum blank feature size constraints;

[0105] Among them, (a) is the optimized structure; (b) is the solid skeleton; (c) is the blank skeleton; (d) is the expanded skeleton of the minimum solid feature size; (e) is the expanded skeleton of the maximum solid feature size; (f) is the expanded skeleton of the minimum blank feature size; (g) is the comparison between the complement of the optimized structure and the expanded skeleton of the minimum solid feature size; (h) is the comparison between the optimized structure and the expanded skeleton of the maximum solid feature size; (i) is the comparison between the optimized structure and the expanded skeleton of the minimum blank feature size;

[0106] Figure 14(a) shows the target value convergence history of the MBB beam;

[0107] Figure 14(b) shows the convergence history of the volume fraction of the MBB beam. DETAILED DESCRIPTION

[0108] Specific implementation method 1: Combination Figure 1 This embodiment describes a method for topological optimization and size control of a beam structure based on a geometric skeleton, and the method specifically includes the following steps:

[0109] Step 1: Mesh the design domain of the beam structure to obtain a performance analysis grid layer; expand the boundaries of the design domain of the beam structure outward (i.e., shift the four boundaries of the performance analysis grid layer outward by margin, and the new rectangular area covering the design domain formed by the four shifted boundaries is the expanded area), mesh the expanded area, and obtain a thermal analysis grid layer;

[0110] The thermal analysis mesh level is larger than the performance analysis mesh level, which provides a margin to avoid divergent numerical anomalies at the boundaries of the design domain.

[0111] The density of each unit in the performance analysis grid layer is used as the design variable of each unit, and the number of iterations is initialized to p = 1, and the design variable value of each unit is initialized;

[0112] Step 2: After filtering the design variables of each unit, substitute the filtering results into the step function to obtain the physical density field;

[0113] Step 3: Calculate the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtain the sensitivity of the target value and volume fraction constraint value respectively;

[0114] Step 4: Obtain the correlation matrix L based on the performance analysis grid layer and the thermal analysis grid layer, and then obtain the temperature field distribution based on the correlation matrix L and the physical density field;

[0115] Step 5: Obtain the divergence distribution according to the temperature field distribution;

[0116] Step 6: extracting the solid skeleton and blank skeleton of the beam structure according to the divergence distribution obtained in step 5 (i.e., the geometric skeleton in the present invention includes the solid skeleton and the blank skeleton), and expanding the solid skeleton and the blank skeleton;

[0117] It should be noted that in variable density topology optimization, especially SIMP, the area with density 1 in the design domain constitutes the solid domain, and the area with density 0 constitutes the blank domain. Therefore, the solid skeleton is the skeleton of the solid domain, and the blank skeleton is the skeleton of the blank domain.

[0118] Constructing minimum entity feature size constraints and maximum entity feature size constraints according to the expanded entity skeleton, constructing minimum blank feature size constraints according to the expanded blank skeleton, and obtaining minimum entity feature size constraint value, maximum entity feature size constraint value, minimum blank feature size constraint value, and sensitivity of minimum entity feature size constraint value, maximum entity feature size constraint value, and minimum blank feature size constraint value;

[0119] Step 7, the target value, volume fraction constraint value, minimum physical feature size constraint value, maximum physical feature size constraint value, minimum blank feature size constraint value, and the sensitivity of the target value, the sensitivity of the volume fraction constraint value, the sensitivity of the minimum physical feature size constraint value, the sensitivity of the maximum physical feature size constraint value, and the sensitivity of the minimum blank feature size constraint value are brought into the optimization solver (moving asymptote solver, etc.) to obtain the updated design variable values ​​of each unit;

[0120] Step 8: Calculate the target change rate based on the target value; and determine whether the iteration stop condition is met based on the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value, and the minimum blank feature size constraint value;

[0121] If the iteration stop condition is met, the iteration is stopped and the physical density field of the size-controlled beam structure is output (i.e., the output of step 2 in the current iteration);

[0122] If the iteration stop condition is not met, set p=p+1 and return to step 2 using the updated design variable values ​​of each unit.

[0123] Specific embodiment 2: This embodiment differs from specific embodiment 1 in that the specific process of step 2 is as follows:

[0124] In the SIMP optimization process, in order to avoid the checkerboard phenomenon, the design variables of unit j are filtered as

[0125]

[0126] Where, is the filter value corresponding to the design variable of unit j, M j The circle is centered at unit j and r f is the set of all units within the radius (it should be noted that to determine whether a unit is in the set M j What needs to be judged is whether the center of this unit is within the circle with unit j as the center and r as the center. f is the radius of the circle), ρ i is the design variable of unit i, w ji is the weight corresponding to the design variable of unit i;

[0127] The vector is composed of the filter values ​​corresponding to the design variables of each unit Then Substituting the step function H in equation (2), we can obtain the physical density field

[0128]

[0129] Where β represents the mapping strength of the step function, and η represents the threshold.

[0130] Other steps and parameters are the same as those in the first embodiment.

[0131] After sorting the filter values ​​corresponding to each column of cells in the design domain, a column vector is obtained In the column vector , the filtered values ​​for the cells in column 1 appear first, followed by the filtered values ​​for the cells in column 2, and so on. Elements with values ​​below the threshold η are in will be mapped to 0, and the elements with values ​​higher than the threshold η will be mapped to will be mapped to 1.

[0132] Specific embodiment 3: This embodiment differs from specific embodiment 1 or 2 in that the weight w corresponding to the design variable of the unit i ji The calculation method is:

[0133] w ji =max (0,r f -dist(j,i)) (3)

[0134] where dist(j,i) is the distance between the center of cell j and the center of cell i.

[0135] Other steps and parameters are the same as those in the first or second embodiment.

[0136] Specific embodiment 4: This embodiment differs from any one of specific embodiments 1 to 3 in that, in step 4, the temperature field distribution is obtained according to the correlation matrix L and the physical density field. The specific process is as follows:

[0137] Set up a homogeneous partial differential equation with Newman boundary conditions:

[0138]

[0139] Where, represents the gradient operator, t is the time of heat conduction, and q represents the temperature;

[0140] Simplifying formula (4) according to time, we get:

[0141]

[0142] Where q (0) represents the temperature at t = 0;

[0143] The temperature in equation (5) is solved by the finite element method:

[0144]

[0145] Where Q represents the heat conduction matrix assembled from the unit matrix; q′ is the vector composed of the temperature of each unit; T is the matrix composed of the node information after converting each unit information into node information; and L is the association matrix between the performance analysis grid layer and the thermal analysis grid layer.

[0146] The other steps and parameters are the same as those in the first to third embodiments.

[0147] It can be assumed that the cells in the performance analysis grid layer (numbered 1 to n, a total of n cells) involve variables The units in the corresponding thermal analysis grid layer (a total of m units) are numbered consecutively from i' to j', and the field after mapping to the thermal analysis grid layer is Based on this, L can be expressed as

[0148]

[0149] Where the subscripts of the matrix elements represent the row and column indices of the elements. It should be noted that the specific form of L can be adjusted according to the actual correspondence between the cell numbers of the two grid layers.

[0150] Specific embodiment 5: This embodiment differs from specific embodiments 1 to 4 in that the specific process of step 5 is as follows:

[0151] Step 51: Calculate the temperature q of each unit j j Gradient in the x direction and the gradient in the y direction

[0152]

[0153] Where, Represents the gradient operators in the x and y directions, N e’ is the e'th shape function of unit j; the present invention takes the quadrilateral unit as an example, a unit has 4 nodes and 4 shape functions, so the value of m is 4. j,e’ is the temperature of the e'th node in element j;

[0154] Step 52: To eliminate the influence of the gradient magnitude, normalize the gradient calculated in step 41 to:

[0155]

[0156] Where, The normalized value of To avoid the divisor being 0, the normalized value of Add a minimum value ∈, the intermediate variable The divergence Δq on unit j is calculated using formula (10) j :

[0157]

[0158] Step 53: Substitute the divergence Δq on unit j j Normalized to:

[0159]

[0160] Where, is the divergence Δq j The normalized value of , where e is the base of the natural logarithm.

[0161] The other steps and parameters are the same as those in the first to fourth embodiments.

[0162] Specific embodiment 6: This embodiment differs from any one of specific embodiments 1 to 5 in that, in step 6, the solid skeleton and the blank skeleton of the beam structure are extracted according to the divergence distribution obtained in step 5, and the solid skeleton and the blank skeleton are expanded; the specific process is as follows:

[0163] The vector composed of the normalized divergence values ​​on each unit is recorded as Will Substitute the step function H into equation (12) to generate the entity skeleton s solid :

[0164]

[0165] Where, β solid is the mapping intensity, η solid is the parameter of the step function;

[0166] The minimum expansion size of the entity skeleton set The kth element in for:

[0167]

[0168] Where γ is the control parameter, e is the base of the natural logarithm, and N k It is centered on unit k and r minsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in the set N k How many units are there in It represents the sum of the number of 1s;

[0169] The maximum inflated size of the entity skeleton set The kth element in for:

[0170]

[0171] Where γ is the control parameter, N k ' is centered on unit k and r maxsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in the set N k 'How many units are there in It represents the sum of the number of 1s;

[0172] Will Substitute the step function H into equation (15) to generate the blank skeleton s void :

[0173]

[0174] Where, β void is the mapping intensity, η void is the parameter of the step function;

[0175] The expanded skeleton set of the blank skeleton The kth element in for:

[0176]

[0177] Where γ is the control parameter, Nk " is centered on unit k and r minvoid is the set of all cells within the circle of the expansion radius, s void(i) It is the skeleton void The i-th element in the set N k How many units are there in "? It represents the sum of the number of 1s.

[0178] The other steps and parameters are the same as those in the first to fifth embodiments.

[0179] Specific embodiment 7: This embodiment differs from any one of specific embodiments 1 to 6 in that the minimum physical feature size constraint is:

[0180]

[0181] Where n represents the number of discrete units, C minsolid It is the minimum entity feature size constraint value, and the superscript T represents transpose.

[0182] The other steps and parameters are the same as those in the first to sixth embodiments.

[0183] The idea of ​​establishing the minimum entity feature size constraint is as follows Figure 3 As shown, the geometric skeleton generated by the method of the present invention is as follows Figure 3 As shown by the red line in (b), based on the generated geometric skeleton, an expanded skeleton is obtained (such as Figure 3 The shaded area in (b). The expanded skeleton and Figure 3 The complement of the original structure of (a) (i.e. Figure 3 The intersection of (c)) in is an empty set, so the mathematical expression of the minimum entity feature size constraint can be established.

[0184] Specific embodiment 8: This embodiment differs from any one of specific embodiments 1 to 7 in that the maximum physical feature size constraint is:

[0185]

[0186] Among them, C maxsolid is the maximum solid feature size constraint value.

[0187] The other steps and parameters are the same as those in the first to seventh embodiments.

[0188] The idea of ​​establishing the maximum entity feature size constraint is as follows Figure 4 As shown, Figure 4 The blue outline area in (b) is obtained by dilating the skeleton (red line). Figure 4 (a) in the equation can be obtained Figure 4(c) in the figure, the complement of the blue outline area is added to Figure 4 Comparing the grayscale areas in (c) with those in (c), we can see that there is no intersection, so we can establish a mathematical expression for the maximum entity feature size constraint.

[0189] Specific embodiment 9: This embodiment differs from any one of specific embodiments 1 to 8 in that the minimum blank feature size constraint is:

[0190] according to Figure 5 (a) in the equation can be obtained Figure 5 (c) in , because Figure 5 The shaded area in (b) is obtained after the blank skeleton is expanded. Figure 5 The shaded areas in (c) have no intersection, so

[0191]

[0192] Among them, C minvoid is the minimum blank feature size constraint value.

[0193] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0194] Specific embodiment ten: This embodiment differs from any one of specific embodiments one to nine in that the target change rate is calculated based on the target value; and whether the iteration stop condition is satisfied is determined based on the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value, and the minimum blank feature size constraint value; the specific process is as follows:

[0195] Step 81: According to the target value obj of the cantilever beam after the current iteration (p) Calculating target rate of change

[0196]

[0197] Where, obj (p-l) is the target value of the cantilever beam after the plth iteration, obj (p-l-1) is the target value of the cantilever beam after the pl-1th iteration, |obj (p) | represents obj (p) The absolute value of l = 0, 1, 2, 3, 4;

[0198] Step 82: Set the iteration stop condition to satisfy both condition (1) and condition (2);

[0199] Condition (1): The target change rate is less than 0.001 for at least five consecutive iterations;

[0200] Condition (2): The minimum entity feature size constraint value of the current iteration is less than the threshold εminsolid , the maximum entity feature size constraint value is less than the threshold ε maxsolid , the minimum blank feature size constraint value is less than the threshold ε minvoid .

[0201] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0202] The following examples will further illustrate the method provided by the present invention. The common parameter configurations are exactly the same. Since the topology optimization is based on the SIMP method, the penalty factor is set to 3. The elastic modulus and minimum elastic modulus of the material are 1 and 10 respectively. -9 , with a Poisson's ratio of 0.3. The initial value of β in the step function expression is set to 1 and then doubled (with an upper limit of 512) after every 50 iterations or when the target change rate falls below a specified threshold (0.001). The time t in equation (5) is calculated as follows:

[0203]

[0204] Figure 2 The margin is (max(r minsolid ,r maxsolid ,r minvoid )+1) three times, that is, the four boundaries of the performance analysis grid layer are all shifted outward by margin, Figure 2 The expansion result of (c) is as follows Figure 2 As shown in (a), Figure 2 The expansion result of (d) is as follows Figure 2 As shown in (b). The parameter γ in equations (13), (14) and (16) is equal to 256. In the process of extracting skeletons of solid areas and blank areas, the parameter γ of the step function is 0.65 (γ solid ) and 0.40(η void ). Mapping intensity β solid and β void The initial value of is 2. When the target change rate is less than 0.001 or every 50 iterations, the mapping strength is doubled (the upper limit is 64). ∈ is assigned to 0.0001. The moving asymptote method (MMA) is used for iterative updates with a step size of 0.01. When the target change rate is less than 0.001 for 5 consecutive iterations and the minimum entity, maximum entity, or minimum blank constraint value is less than ε minsolid , ε maxsolid and ε minvoid When , the convergence condition is met.

[0205] Case 1: Taking the optimization of a cantilever beam as an example, the cantilever beam is optimized with the goal of minimum flexibility (i.e., maximum stiffness). The upper limit of the volume fraction is 0.3, and the desired minimum entity size or maximum entity size needs to be achieved. The design domain and boundary conditions to be optimized are as follows: Figure 6 As shown, the left edge is completely fixed and the right edge is loaded in the middle with f ext = 1. The entire region is discretized using quadrilateral elements, with a total of 400 × 200 elements. The initial design variable ρ is 0.3, and the filter radius is 7.

[0206] In order to verify the effectiveness of the minimum entity feature size, the minimum entity size radius is 3. minsolid Set to 0.01. After 131 iterations, the optimization results of the design domain are as follows Figure 7 As shown in (a) in FIG. Using the method proposed by the present invention, it can be clearly seen that Figure 7 The skeleton shown in (b) in the figure is shown in Figure 1. In addition, the green point in the figure indicates the minimum entity feature size. Figure 7 As shown in (c) in . Figure 7 Figure 8(d) shows the difference between the optimized structure and the expanded skeleton. The figure shows only a few differences, confirming the effectiveness of the minimum entity feature size constraint. The final optimization target value is 114.73, and the volume fraction reaches its upper limit of 0.30. The target history curve in Figure 8(a) shows a fast and stable optimization process, demonstrating the powerful capability of the algorithm. The volume fraction history curve in Figure 8(b) shows only minor fluctuations, further confirming the effectiveness of the proposed method.

[0207] In order to verify the effectiveness of the maximum solid feature size constraint, the same cantilever beam was analyzed with the expansion radius set to 7 and ε maxsolid Set to 0.007, and keep other parameters unchanged. The topology optimization results are as follows Figure 9 As shown in (a) in Figure 2, comparing the green points with the black structure, it can be seen that the maximum entity feature size has been basically met. Figure 9 (b) shows the geometric skeleton of the structure. Figure 9 The expanded skeleton in (c) and Figure 9 Comparing the optimized structure of (a) in the figure, we can get Figure 9 (d) in the figure strongly demonstrates that almost the entire structure satisfies the constraints. After 135 iterations, the target value is 108.11. The volume fraction reaches its upper limit of 0.30. The historical curves of target value and volume fraction, as shown in Figures 10(a) and 10(b), demonstrate that the proposed method is effective and stable for the maximum solid feature size problem.

[0208] Case 2: The following example uses an MBB beam to demonstrate the effect of the minimum blank feature size. Minimum flexibility is used as the objective function, and the upper limit of the volume fraction is set to 0.20. The design domain and boundaries are as follows: Figure 11 The entire area is discretized into 600×200 quadrilateral elements, and an external load f is set at the center of the top edge. ext is 1. The filter radius is set to 7. The volume fraction is 0.32. The initial design variable is 0.32. For the convenience of comparison, the MBB structure with the minimum and maximum entity feature sizes is first solved, and its expansion radius is 4 and 7 respectively. The optimization results are shown in Figure 12 As shown, the objective function value is 37.54 and the volume fraction is 0.32.

[0209] In order to demonstrate the ability to control the minimum blank feature size, we add corresponding constraints based on the above case and set the minimum blank feature expansion radius to 6. Other parameters are exactly the same as the previous MBB. The optimization results are as follows: Figure 13 As shown, compared Figure 12 and Figure 13 In (a), it can be seen that the gap near the load area becomes significantly larger. The intersection area of ​​the structure will have obvious rounded corners. In addition, Figure 13 It can be clearly seen in (c) that the skeleton of the blank area is extracted, but some skeletons are still missing. This is because the size of the blank area in some local areas will change significantly, resulting in a relatively small divergence. Fortunately, these areas have met the blank feature size requirements. When the local area is small, the skeleton is clear and the feature size is effectively controlled. The final target value obtained is 38.81, which is greater than the example without the minimum blank feature size constraint, and the volume fraction is 0.32. The historical curves of Figures 14(a) and 14(b) also show the stability of the entire optimization process.

[0210] The above examples are merely illustrative of the calculation model and process of the present invention and are not intended to limit the embodiments of the present invention. Persons skilled in the art will readily appreciate that other variations or modifications based on the above description are possible. This list of embodiments is not exhaustive; however, any obvious variations or modifications derived from the technical solution of the present invention remain within the scope of protection of the present invention.

Claims

1. A method for topological optimization and size control of beam structures based on a geometric skeleton, characterized in that: The method specifically comprises the following steps: Step 1: Mesh the design domain of the beam structure to obtain a performance analysis mesh layer; expand the boundary of the design domain of the beam structure outward, mesh the expanded area, and obtain a thermal analysis mesh layer; The density of each unit in the performance analysis grid layer is used as the design variable of each unit, and the number of iterations is initialized to p = 1, and the design variable value of each unit is initialized; Step 2: After filtering the design variables of each unit, substitute the filtering results into the step function to obtain the physical density field; Step 3: Calculate the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtain the sensitivity of the target value and volume fraction constraint value respectively; Step 4: Obtain the correlation matrix L based on the performance analysis grid layer and the thermal analysis grid layer, and then obtain the temperature field distribution based on the correlation matrix L and the physical density field; The temperature field distribution is obtained according to the correlation matrix L and the physical density field. The specific process is: Set up the partial differential equation: Where, represents the gradient operator, t is the time of heat conduction, and q represents the temperature; Simplifying formula (4) according to time, we get: Where q (0) represents the temperature at t = 0; The temperature in equation (5) is solved by the finite element method: Where Q represents the heat conduction matrix; q′ is the vector composed of the temperature of each unit; T is the matrix composed of node information after converting each unit information into node information; L is the correlation matrix between the performance analysis grid layer and the thermal analysis grid layer, represents the physical density field; Step 5: Obtain the divergence distribution according to the temperature field distribution; the specific process is: Step 51: Calculate the temperature q of each unit j j Gradient in the x direction and the gradient in the y direction Where, and Represents the gradient operators in the x and y directions, N e, is the e'th shape function of element j; q j,e’ is the temperature of the e'th node in element j; Step 52: Normalize the gradient calculated in step 51 and calculate the divergence Δq on unit j using formula (10): j : in, yes The normalized value of yes The normalized value of Step 53: Substitute the divergence Δq on unit j j Normalization; Step 6: extracting the solid skeleton and blank skeleton of the beam structure according to the divergence distribution obtained in step 5, and expanding the solid skeleton and the blank skeleton; The solid skeleton and blank skeleton of the beam structure are extracted according to the divergence distribution obtained in step 5, specifically: The vector composed of the normalized divergence values ​​on each unit is recorded as Will Substitute the step function H into equation (12) to generate the entity skeleton s solid : Where, β solid is the mapping intensity, η solid is the parameter of the step function; Will Substitute the step function H into equation (15) to generate the blank skeleton s void : Where, β void is the mapping intensity, η void is the parameter of the step function; Constructing minimum entity feature size constraints and maximum entity feature size constraints according to the expanded entity skeleton, constructing minimum blank feature size constraints according to the expanded blank skeleton, and obtaining minimum entity feature size constraint value, maximum entity feature size constraint value, minimum blank feature size constraint value, and sensitivity of minimum entity feature size constraint value, maximum entity feature size constraint value, and minimum blank feature size constraint value; Step 7: The target value, volume fraction constraint value, minimum physical feature size constraint value, maximum physical feature size constraint value, minimum blank feature size constraint value, and the sensitivity of the target value, the sensitivity of the volume fraction constraint value, the sensitivity of the minimum physical feature size constraint value, the sensitivity of the maximum physical feature size constraint value, and the sensitivity of the minimum blank feature size constraint value are brought into the optimization solver to obtain the updated design variable value of each unit; Step 8: Calculate the target change rate based on the target value; and determine whether the iteration stop condition is met based on the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value, and the minimum blank feature size constraint value; If the iteration stop condition is met, the iteration is stopped and the physical density field of the size-controlled rear beam structure is output; If the iteration stop condition is not met, set p=p+1 and return to step 2 using the updated design variable values ​​of each unit.

2. The method for topological optimization and size control of a beam structure based on a geometric skeleton according to claim 1, characterized in that: The specific process of step 2 is: Filter the design variables of unit j to Where, is the filter value corresponding to the design variable of unit j, M j The circle is centered at unit j and r f is the set of all units within the radius, ρ i is the design variable of unit i, w ji is the weight corresponding to the design variable of unit i; The vector is composed of the filter values ​​corresponding to the design variables of each unit Then Substituting the step function H in equation (2), we can obtain the physical density field Where β represents the mapping strength of the step function, and η represents the threshold.

3. The method for topological optimization and size control of a beam structure based on a geometric skeleton according to claim 2, characterized in that: The weight w corresponding to the design variable of unit i ji The calculation method is: w ji =max(0,r f -dist(j,i)) (3) where dist(j,i) is the distance between the center of cell j and the center of cell i.

4. The method for topological optimization and size control of a beam structure based on a geometric skeleton according to claim 3, characterized in that: The gradient calculated in step 41 is normalized as follows: Where, yes The normalized value of yes Normalized value of the intermediate variable The divergence Δq on unit j is j Normalization, specifically: Where, is the divergence Δq j The normalized value of , where e is the base of the natural logarithm.

5. The method for topological optimization and size control of beam structures based on a geometric skeleton according to claim 4, characterized in that: In step 6, the solid skeleton and the blank skeleton are expanded; the specific process is: The minimum inflated size of the entity skeleton set The kth element in for: Where γ is the control parameter, e is the base of the natural logarithm, and N k It is centered on unit k and r minsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in ; The maximum inflated size of the entity skeleton set The kth element in for: Where γ is the control parameter, N k ' is centered on unit k and r maxsolid is the set of all cells within the circle of the expansion radius, s solid(i) Is the entity skeleton solid The i-th element in ; Expanded skeleton collection of blank skeletons The kth element in for: Where γ is the control parameter, N k " is centered on unit k and r minvoid is the set of all cells within the circle of the expansion radius, s void(i) It is the skeleton void The i-th element in .

6. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 5, characterized in that: The minimum physical feature size constraint is: Where n represents the number of discrete units, C minsolid It is the minimum entity feature size constraint value, and the superscript T represents transpose.

7. The method for topological optimization and size control of beam structures based on a geometric skeleton according to claim 6, characterized in that: The maximum physical feature size constraint is: Among them, C maxsolid is the maximum solid feature size constraint value.

8. The method for topological optimization and size control of beam structures based on a geometric skeleton according to claim 7, characterized in that: The minimum blank feature size constraint is: Among them, C minvoid is the minimum blank feature size constraint value.

9. The method for topological optimization and size control of beam structures based on a geometric skeleton according to claim 8, characterized in that: The target change rate is calculated according to the target value; and whether the iteration stop condition is satisfied is determined according to the target change rate, the minimum entity feature size constraint value, the maximum entity feature size constraint value and the minimum blank feature size constraint value; the specific process is: Step 81: According to the target value obj of the cantilever beam after the current iteration (p) Calculating target rate of change Where, obj (p-l) is the target value of the cantilever beam after the plth iteration, obj (p-l-1) is the target value of the cantilever beam after the pl-1th iteration, |obj (p) | represents obj (p) The absolute value of l = 0, 1, 2, 3, 4; Step 82: Set the iteration stop condition to satisfy both condition (1) and condition (2); Condition (1): The target change rate is less than 0.001 for at least five consecutive iterations; Condition (2): The minimum entity feature size constraint value of the current iteration is less than the threshold ε minsolid , the maximum entity feature size constraint value is less than the threshold ε maxsolid , the minimum blank feature size constraint value is less than the threshold ε minvoid .

Images