Geometric skeleton-based topological optimization size regulation and control method for beam type structure

Abstract

The invention discloses a method for regulating and controlling the topological optimization size of a beam type structure based on a geometric framework, and belongs to the technical field of topological optimization of beam type structures. According to the method, the problem that the size regulation result of an existing beam type structure size regulation method based on a geometric framework cannot be optimal is solved. The method comprises the steps of firstly simulating a heat conduction behavior through a partial differential equation, then solving a temperature distribution state, then obtaining divergence distribution based on temperature distribution, and finally identifying a geometric skeleton through the numerical value of divergence. A geometric skeleton is expanded according to a specified size radius, a contrast set of required feature sizes is established, the formed geometry is compared with an original structure, feature size constraints including the minimum entity feature size, the maximum entity feature size and the minimum blank feature size are constructed, and the feature size constraints are obtained. And effective regulation and control of the minimum entity feature size, the maximum entity feature size and the minimum blank feature size can be realized. The method can be applied to topological optimization size regulation and control of the beam type structure.

Description

Technical Field

[0001] The 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] Topological optimization is a powerful tool for beam structure design, but the structures it designs are often too complex and slender, which poses a huge challenge to processing and manufacturing, and hinders the transformation of design solutions into actual products. In order to solve this problem, the current focus is on the following two aspects: on the one hand, developing feature size control methods to avoid too many branches in the structural solution; on the other hand, developing design methods compatible with manufacturing processes so that the design solution can meet the processing requirements. Feature size control technology has significant advantages in versatility.

[0003] Among the many developed topology optimization methods, the variable density method has strong application flexibility and stability due to the scheme of continuous design variables 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. The size control for this method is usually carried out for the variable filtering and step mapping links in the optimization. However, this method often has the characteristics of weak control ability, unstable method, and implicit expression. In addition to the above methods, size control methods based on geometric skeletons have also received widespread attention. However, geometric skeleton extraction is often based on graphics methods, without considering the continuous design variable characteristics of topology optimization, which leads to a serious loss of sensitivity information, thereby significantly reducing the optimality of the size 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, meshing the design domain of the beam structure to obtain a performance analysis mesh layer; expanding the boundary of the design domain of the beam structure outward, meshing the expanded area, and obtaining 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 p is initialized to 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, calculating the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtaining the sensitivity of the target value and the volume fraction constraint value respectively;

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

[0011] Step 5, obtaining the divergence distribution according to the temperature field distribution;

[0012] Step 6: extracting the solid skeleton and the 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 values, maximum entity feature size constraint values, minimum blank feature size constraint values, and the sensitivity of the minimum entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value;

[0014] Step 7, bringing the target value, volume fraction constraint value, minimum entity feature size constraint value, maximum entity 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 entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value into the optimization solver to obtain the updated design variable values ​​of each unit;

[0015] Step 8, calculating the target change rate according to the target value; and judging whether the iteration stop condition is met 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;

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

[0017] If the iteration stop condition is not met, p=p+1 is set, and the execution returns 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] In the formula, 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 cells 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 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, and the specific process is:

[0029] Set up the partial differential equation:

[0030]

[0031] In the formula, 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] In the formula, 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 composed of the temperature of each unit; T is the matrix composed of the node information after converting the unit information into node information; 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 unit j respectively j The gradient in the x direction and the gradient in the y direction

[0040]

[0041] In the formula, Represents the gradient operators in the x and y directions respectively, 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] In the formula, The normalized value of The normalized value of the intermediate variable

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

[0046]

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

[0048]

[0049] In the formula, 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] In the formula, β 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] In the formula, γ is the control parameter, e is the base of the natural logarithm, and N k The unit k is the center and r is the 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] In the formula, γ 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] In the formula, β 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] In the formula, γ 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 entity feature size constraint is:

[0067]

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

[0069] Furthermore, the maximum entity 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 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:

[0076] Step 81: According to the target value obj of the cantilever beam after the current iteration (p) Calculate 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) | indicates obj (p) The absolute value of l = 0, 1, 2, 3, 4;

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

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

[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 simulates the heat conduction behavior by 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. The partial differential equations are used to simulate the heat conduction process, which 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 the characteristic size, the present invention expands the extracted skeleton with a specified size radius, establishes a control set of the required characteristic size, compares the formed geometry with the original structure, and constructs a characteristic size constraint including the minimum entity characteristic size, the maximum entity characteristic size, and the minimum blank characteristic size, which 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 size of the entity domain and the blank domain. The method of the present invention can adapt to the current existing variable density topology optimization method and has good versatility. BRIEF DESCRIPTION OF THE DRAWINGS

[0084] Figure 1 It is a framework diagram generated by the skeleton 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] Among them, (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 constraints;

[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 entity feature size;

[0098] Fig. 9 It is the topological 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 entity feature size;

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

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

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

[0104] Fig.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 The present embodiment describes a method for controlling the size of a beam structure based on a geometric skeleton topology optimization, 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 boundary of the design domain of the beam structure outward (that is, the four boundaries of the performance analysis grid layer are all translated outward by margin, and the new rectangular area covering the design domain surrounded by the four translated 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 design domain boundaries.

[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 p is initialized to 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, calculating the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtaining the sensitivity of the target value and the volume fraction constraint value respectively;

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

[0115] Step 5, obtaining 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 (that is, 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 values, maximum entity feature size constraint values, minimum blank feature size constraint values, and the sensitivity of the minimum entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value;

[0119] Step 7, bringing the target value, volume fraction constraint value, minimum entity feature size constraint value, maximum entity 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 entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value into the optimization solver (moving asymptote method solver, etc.) to obtain the updated design variable values ​​of each unit;

[0120] Step 8, calculating the target change rate according to the target value; and judging whether the iteration stop condition is met 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;

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

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

[0123] Specific implementation method 2: This implementation method is different from the specific implementation method 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] In the formula, 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 obtain the physical density field

[0128]

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

[0130] The 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 η in will be mapped to 0, and elements with values ​​higher than the threshold η will be mapped to will be mapped to 1.

[0132] Specific implementation method three: This implementation method is different from specific implementation methods one or two 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] The other steps and parameters are the same as those in the first or second embodiment.

[0136] Specific implementation method 4: This implementation method is different from any one of specific implementation methods 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] In the formula, 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] In the formula, 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 by 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; 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 Specific Embodiments 1 to 3.

[0147] It can be assumed that the cells (numbered 1 to n, a total of n cells) in the performance analysis grid layer involve variables The cells in the corresponding thermal analysis grid layer (a total of m cells) 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] In the formula, the subscripts of the matrix elements represent the indexes of the rows and columns of the elements. It should be noted that the specific form of L can be adjusted according to the actual correspondence of the unit numbers of the two layers of grids.

[0150] Specific implementation method 5: This implementation method is different from specific implementation methods 1 to 4 in that the specific process of step 5 is as follows:

[0151] Step 51, calculate the temperature q of unit j respectively j The gradient in the x direction and the gradient in the y direction

[0152]

[0153] In the formula, Represents the gradient operators in the x and y directions respectively, 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. q j,e’ is the temperature of the e'th node in element j;

[0154] Step 52: In order to eliminate the influence of the gradient magnitude, the gradient calculated in step 41 is normalized to:

[0155]

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

[0157]

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

[0159]

[0160] In the formula, 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 Specific Embodiments 1 to 4.

[0162] Specific implementation method 6: This implementation method is different from any one of specific implementation methods 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:

[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] In the formula, β 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] In the formula, γ is the control parameter, e is the base of the natural logarithm, and N k The unit k is the center and r is the 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] In the formula, γ 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] In the formula, β 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] In the formula, γ 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 1s.

[0178] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0179] Specific implementation method 7: This implementation method is different from any one of specific implementation methods 1 to 6 in that the minimum entity feature size constraint is:

[0180]

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

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

[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 implementation example eight: This implementation example is different from any one of specific implementation examples one to seven in that the maximum entity 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 Specific Embodiments 1 to 7.

[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 4In (c), the complement of the blue outline area is Figure 4 By comparing the grayscale area in (c) with that in (c), it can be seen that there is no intersection, so a mathematical expression for the maximum entity feature size constraint can be established.

[0189] Specific implementation method 9: This implementation method is different from any one of specific implementation methods 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 area in (c) has 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 implementation method ten: This implementation method is different from the specific implementation methods one to nine 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:

[0195] Step 81: According to the target value obj of the cantilever beam after the current iteration (p) Calculate 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) | indicates obj (p) The absolute value of l = 0, 1, 2, 3, 4;

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

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

[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 example 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 , Poisson's ratio is 0.3. The initial value of β in the step function expression is set to 1, and then doubled (up to 512) after every 50 iterations or when the target change rate is lower than the specified threshold (0.001). The time t in equation (5) is calculated by the following formula:

[0203]

[0204] Figure 2 The margin is (max(r minsolid ,r maxsolid ,r minvoid )+1), 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 The parameter γ in equations (13), (14) and (16) is equal to 256. In the skeleton extraction process of the solid area and the blank area, the parameter γ of the step function is 0.65 (γ solid ) and 0.40(η void ). Mapping strength β 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 will be doubled (the upper limit is 64). ∈ is assigned to 0.0001. The moving asymptote method (MMA) is used for iterative update, and its step size is set to 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 , 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 minimum flexibility (i.e., maximum stiffness) as the goal. 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 a load f ext = 1. The entire area 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 is set to 0.01. After 131 iterations, the optimization results of the design domain are as follows Figure 7 As shown in (a) of FIG. Utilizing the method proposed by the present invention, it can be clearly seen that Figure 7 The skeleton shown in (b) in the figure, in addition, the green point in the figure indicates the minimum entity feature size. Figure 7 As shown in (c) in . Figure 7 (d) shows the difference between the optimized structure and the expanded skeleton. There are only a few differences in the figure, which confirms 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 historical curve of the target in Figure 8(a) shows a fast and stable optimization process, showing the powerful ability of the algorithm. The volume fraction history curve in Figure 8(b) shows only slight fluctuations, further confirming the effectiveness of the method of the present invention.

[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 Fig. 9 As shown in (a) in the figure, comparing the green dots with the black structure, it can be seen that the maximum entity feature size has been basically met. Fig. 9 (b) shows the geometric skeleton of the structure. Fig. 9 The expanded skeleton in (c) and Fig. 9 Comparing with the optimized structure of (a), we can get Fig. 9 (d) in Fig. 10 strongly proves that almost the entire structure satisfies the constraints. After 135 iterations, the target value is 108.11. The volume fraction reaches the upper limit of 0.30. As shown in Fig. 10(a) and Fig. 10(b), the historical curves of the target value and the volume fraction show that the method proposed in the present invention is effective and stable for the maximum entity feature size problem.

[0208] Case 2: The following example will use the MBB beam to demonstrate the effect of the minimum blank feature size. The 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 Fig.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 Fig.12 As shown, the objective function value is 37.54 and the volume fraction is 0.32.

[0209] In order to reflect 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: Fig.13 As shown, compared Fig.12 and Fig.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, Fig.13 It can be clearly seen in (c) that the skeleton of the blank area is extracted, but part of the skeleton is still missing. This is because the size of the blank area will change significantly in some local areas, 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 is 38.81, which is larger than the example without the minimum blank feature size constraint, and the volume fraction is 0.32. The historical curves of Figure 14 (a) and Figure 14 (b) also show the stability of the entire optimization process.

[0210] The above calculation examples of the present invention are only used to explain the calculation model and calculation process of the present invention in detail, and are not intended to limit the implementation methods of the present invention. For ordinary technicians in the relevant field, other different forms of changes or modifications can be made based on the above description. It is impossible to list all the implementation methods here. All obvious changes or modifications derived from the technical solution of the present invention are still 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, meshing the design domain of the beam structure to obtain a performance analysis mesh layer; expanding the boundary of the design domain of the beam structure outward, meshing the expanded area, and obtaining 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 p is initialized to 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, calculating the target value and volume fraction constraint value of the beam structure according to the physical density field, and obtaining the sensitivity of the target value and the volume fraction constraint value respectively; Step 4: Obtain the correlation matrix L according to the performance analysis grid layer and the thermal analysis grid layer, and then obtain the temperature field distribution according to the correlation matrix L and the physical density field; Step 5, obtaining the divergence distribution according to the temperature field distribution; Step 6: extracting the solid skeleton and the blank skeleton of the beam structure according to the divergence distribution obtained in step 5, and expanding the solid skeleton and the blank skeleton; 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 values, maximum entity feature size constraint values, minimum blank feature size constraint values, and the sensitivity of the minimum entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value; Step 7, bringing the target value, volume fraction constraint value, minimum entity feature size constraint value, maximum entity 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 entity feature size constraint value, the sensitivity of the maximum entity feature size constraint value, and the sensitivity of the minimum blank feature size constraint value into the optimization solver to obtain the updated design variable values ​​of each unit; Step 8, calculating the target change rate according to the target value; and judging whether the iteration stop condition is met 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; If the iteration stop condition is met, the iteration is stopped and the physical density field of the beam structure after size control is output; If the iteration stop condition is not met, p=p+1 is set, and the execution returns to step 2 using the updated design variable values ​​of each unit.

2. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 1 is characterized in that: The specific process of step 2 is: Filter the design variables of unit j to In the formula, 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 cells 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 obtain the physical density field Where β represents the mapping strength of the step function and η represents the threshold.

3. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 2 is 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 controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 3 is characterized 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: Set up the partial differential equation: In the formula, represents the gradient operator, t is the time of heat conduction, and q represents the temperature; Simplifying formula (4) according to time, we get: In the formula, 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 the node information after converting the unit information into node information; L is the association matrix between the performance analysis grid layer and the thermal analysis grid layer.

5. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 4 is characterized in that: The specific process of step 5 is as follows: Step 51, calculate the temperature q of unit j respectively j The gradient in the x direction and the gradient in the y direction In the formula, and Represents the gradient operators in the x and y directions respectively, 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 41 to: In the formula, yes The normalized value of yes The normalized value of the intermediate variable The divergence Δq on cell j is calculated using formula (10): j : Step 53: Substitute the divergence Δq on unit j j Normalized to: In the formula, is the divergence Δq j The normalized value of , where e is the base of the natural logarithm.

6. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 5 is characterized 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: 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 : In the formula, β solid is the mapping intensity, η solid is the parameter of the step function; The minimum expansion size of the entity skeleton set The kth element in for: In the formula, γ is the control parameter, e is the base of the natural logarithm, and N k The unit k is the center and r is the 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: In the formula, γ 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 ; Will Substitute the step function H into equation (15) to generate the blank skeleton s void : In the formula, β void is the mapping intensity, η void is the parameter of the step function; The expanded skeleton set of the blank skeleton The kth element in for: In the formula, γ 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 .

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

8. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 7 is characterized in that: The maximum entity feature size constraint is: Among them, C maxsolid is the maximum solid feature size constraint value.

9. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 8, characterized in that: The minimum blank feature size constraint is: Among them, C minvoid is the minimum blank feature size constraint value.

10. The method for controlling the size of a beam structure topology optimization based on a geometric skeleton according to claim 9, characterized in that: The target change rate is calculated according to the target value; and whether the iteration stop condition is met 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) Calculate 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) | indicates obj (p) The absolute value of l = 0, 1, 2, 3, 4; Step 82, setting the iteration stop condition to satisfy both condition (1) and condition (2); Condition (1): The target change rate for at least five consecutive iterations is less than 0.001; 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