Multi-material topological optimization method based on univariate discrete interpolation scheme

Abstract

The invention relates to a multi-material topological optimization method based on a univariate discrete interpolation scheme, belongs to the technical field of multi-material topological optimization, and solves the problems of multiple design variables, low optimization efficiency, unreasonable design result and difficulty in processing and manufacturing of a designed structure in multi-material topological optimization in the prior art. Comprising the following steps: performing finite element discretization on a design domain, and applying a boundary condition and an external load; based on a univariate discrete combination interpolation method, a multi-material topological optimization model is constructed, and then the sensitivity of the multi-material topological optimization model is obtained; and on the basis of the multi-material topological optimization model and the sensitivity, updating a design variable to obtain an optimal topological density field. And efficient and reasonable multi-material topological optimization is realized.

Description

Technical Field

[0001] The invention relates to the technical field of structural optimization, and in particular to a multi-material topology optimization method based on a single variable discrete interpolation scheme. Background Art

[0002] Topological optimization refers to calculating the optimal distribution of materials through algorithms within a given design domain based on certain constraints (such as material volume, material mass, etc.), thereby obtaining an efficient structural shape, which can effectively reduce structural weight and improve structural performance.

[0003] In recent years, with the rapid development of additive manufacturing technology (such as 3D printing), structural topology optimization design has become one of the hot topics in engineering and scientific research. In the past two or three decades, topology optimization methods have become more and more diverse, covering a variety of methods such as solid isotropic material penalty (SIMP) method, level set method, bidirectional evolutionary structural optimization method, mobile deformable cavity method and phase field method. It is worth noting that among density-based methods, SIMP method is favored due to its high efficiency and easy implementation. With the continuous advancement of materials science and manufacturing technology, traditional single-material topology optimization methods can no longer meet the needs of structural versatility and composite performance, and topology optimization has gradually expanded to the field of multi-materials. Compared with single-material structures, multi-material structures can better meet complex functional requirements and significantly improve structural performance through the reasonable combination and distribution of materials. However, the topology optimization of multi-material structures still faces many challenges, especially in terms of the number of design variables, the clarity of material interfaces, and the rationality of material distribution.

[0004] In summary, the traditional multi-material topology optimization method has the problem of "combinatorial explosion" caused by the increase of design variables as the number of alternative materials increases. The number of design variables for each design unit is proportional to the number of material phases, resulting in many design variables and low optimization efficiency. The existing multi-material topology optimization step interpolation model will have unreasonable designs such as material envelopes, which greatly reduces the optimization design space. The feasible domain of the design is significantly smaller than the feasible domain of the discrete material solution, resulting in unreasonable design results and difficult to process and manufacture the designed structure. Summary of the invention

[0005] In view of the above analysis, an embodiment of the present invention aims to provide a multi-material topology optimization method based on single-variable discrete interpolation, so as to solve the problems in the existing multi-material topology optimization, such as many design variables, low optimization efficiency, unreasonable design results, and difficult processing and manufacturing of the designed structure.

[0006] An embodiment of the present invention provides a multi-material topology optimization method based on a univariate discrete interpolation scheme, comprising the following steps:

[0007] The design domain is discretized by finite element method and boundary conditions and external loads are applied;

[0008] Based on the single variable discrete combined interpolation method, a multi-material topology optimization model is constructed, and then the sensitivity of the multi-material topology optimization model is obtained;

[0009] Based on the multi-material topology optimization model and sensitivity, the design variables are updated to obtain the optimal topological density field.

[0010] Furthermore, based on the single variable discrete combined interpolation method, a multi-material topology optimization model is constructed including:

[0011] Based on the univariate characteristic function, the design variables are mapped into the topological density field of the finite element unit;

[0012] Based on the topological density field of the finite element unit, a filtered topological density field of the finite element unit is obtained, and then a normalized topological density field of the finite element unit is obtained;

[0013] Based on the normalized topological density field of the finite element unit, the elastic modulus of the candidate material is discretely interpolated to obtain the elastic modulus of the finite element unit, and then the elastic tensor of the finite element unit is obtained;

[0014] Based on the elastic tensor of the finite element unit, a multi-material topology optimization model is constructed.

[0015] Furthermore, the univariate characteristic function is expressed as:

[0016]

[0017] In the formula, It means that when the design variable of the ξth finite element is χ(ξ), The topological density of the phase candidate material, It means that when the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of the phase candidate material, η represents the penalty function, Θ Indicates the number of phases of the candidate material.

[0018] Furthermore, when the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of phase candidate materials It is expressed as:

[0019]

[0020] In the formula, Indicates is the threshold parameter for the phase candidate material, β represents the projection parameter, and tanh() represents the hyperbolic tangent function.

[0021] Furthermore, the elastic modulus and elastic tensor of the finite element unit are respectively expressed as:

[0022]

[0023] Where, E[χ(ξ)] and D[χ(ξ)] represent the elastic modulus and elastic tensor of the ξ-th finite element when the design variable is χ(ξ); denote the elastic modulus and elastic tensor of the empty material, respectively. Respectively represent Elastic modulus and elastic tensor of the candidate materials, Indicates that in the ξth finite element The weight coefficient of the candidate material.

[0024] Furthermore, the ξth finite element unit Weight coefficient of phase candidate materials It is expressed as:

[0025]

[0026] In the formula, They represent the Normalized topological density of k-phase candidate materials.

[0027] Furthermore, based on the multi-material topology optimization model and sensitivity, the gradient optimization algorithm MMA is used to update the design variables until the multi-material topology optimization model converges to obtain the optimal design variables and then the optimal topological density field.

[0028] Furthermore, when the design variables are updated by the gradient optimization algorithm MMA, the optimization parameters in the iteration are updated based on the adaptive parameter adjustment strategy; wherein the optimization parameters include penalty parameters and projection parameters.

[0029] Furthermore, the adaptive parameter adjustment strategy includes:

[0030] Before the gradient optimization algorithm MMA iteration, the optimization parameters are initialized and the counter value is initialized to 0;

[0031] During the MMA iteration of the gradient optimization algorithm, it is determined whether the projection parameter in the current iteration is less than the set projection parameter threshold. If so, the design variable in the current iteration is used as the optimal design variable; otherwise:

[0032] Determine whether the value of the counter in the current iteration is greater than 0,

[0033] If yes, determine whether the change between the objective function in the current iteration and the objective function in the previous iteration is less than the predetermined change value. If yes, set the value of the counter to 1, and do not update the optimization parameters in the current iteration; if no, do not update the optimization parameters in the current iteration;

[0034] If not, the counter value is increased by 1, and it is determined whether the counter value at this time is the set number threshold. If so, the optimization parameters in the current iteration are updated and the counter value is set to 1; if not, the optimization parameters in the current iteration are not updated.

[0035] Furthermore, the projection parameters in the optimization parameters are updated by the following formula:

[0036] β=θ×β′

[0037] Where β′ is the projection parameter in the current iteration, β is the updated projection parameter, and θ represents the first adjustment coefficient.

[0038] The penalty parameter in the optimization parameter is updated by the following formula:

[0039]

[0040] Where η′ is the penalty parameter in the current iteration, η is the updated penalty parameter, and σ represents the second adjustment coefficient.

[0041] Compared with the prior art, the present invention can achieve at least one of the following beneficial effects:

[0042] The present invention provides a multi-material topology optimization method based on a single variable discrete interpolation scheme. The method discretizes the design domain by finite elements, applies boundary conditions and external loads, and then constructs a multi-material topology optimization model based on a single variable discrete combined interpolation method, thereby obtaining the sensitivity of the multi-material topology optimization model. Based on the multi-material topology optimization model and the sensitivity, the design variables are updated to obtain the optimal topological density field, thereby solving the problems of many design variables, low optimization efficiency, unreasonable design results, and difficult processing and manufacturing of the designed structure. The method constructs a single variable characteristic function, uses a single variable to select multiple discrete attribute physical quantities, and combines the composite function of the single variable reduced-order multi-characteristic function with the corresponding material attribute for material interpolation, thereby eliminating the material envelope phenomenon existing in the existing step-type interpolation model and overcoming the disadvantage of limited design domain. The method makes the number of design variables independent of the number of candidate material phases, significantly reduces the number of design variables, and reduces the number of iterations while achieving similar or even lower objective function values, thereby improving the optimization efficiency. The method can realize the reasonable distribution of materials in the structure under the premise of ensuring the constraints of material quality or volume fraction.

[0043] In the present invention, the above-mentioned technical solutions can also be combined with each other to achieve more preferred combination solutions. Other features and advantages of the present invention will be described in the subsequent description, and some advantages can become obvious from the description, or can be understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained through the contents particularly pointed out in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] The drawings are only for the purpose of illustrating particular embodiments and are not to be considered limiting of the present invention. Like reference symbols denote like components throughout the drawings.

[0045] Figure 1 A schematic flow chart of a multi-material topology optimization method based on a univariate discrete interpolation scheme provided in Example 1 of the present invention;

[0046] Figure 2 A schematic diagram of a specific process of a multi-material topology optimization method based on a univariate discrete interpolation scheme provided in Example 1 of the present invention;

[0047] Figure 3 A schematic diagram of the conditions of a two-dimensional MBB beam in scenario 1 provided in embodiment 2 of the present invention;

[0048] Figure 4 A schematic diagram of the iterative convergence process of the objective function value and volume fraction of the MBB beam structure in scenario 1 provided by embodiment 2 of the present invention;

[0049] FIG5( a ) is a comparison diagram of the two-dimensional MBB beam topology optimization results obtained by the proposed method and the DMO method under scenario 1 provided in Example 2 of the present invention when θ=1.24, σ=1.065;

[0050] FIG5( b ) is a comparison diagram of the two-dimensional MBB beam topology optimization results obtained by the proposed method and the DMO method under scenario 1 provided in Example 2 of the present invention when θ=1.1445 and σ=1.0247;

[0051] Figure 6 A schematic diagram of the conditions of the cantilever beam structure under scenario 2 provided in embodiment 2 of the present invention;

[0052] Figure 7 A graph showing the change trend of the objective function and volume fraction with the number of iterations under 5 material quantities in scenario 2 provided in embodiment 2 of the present invention;

[0053] FIG8( a ) is a comparison diagram of the two-dimensional MBB beam topology optimization results obtained by the proposed method for five materials under scenario 2 provided in Example 2 of the present invention and the DMO method;

[0054] FIG8( b ) is a comparison diagram of the two-dimensional MBB beam topology optimization results obtained by the proposed method with three materials under scenario 1 provided in Example 2 of the present invention and the DMO method;

[0055] Fig. 9 A schematic diagram of the conditions of a bridge under scenario 3 provided in embodiment 2 of the present invention;

[0056] Fig.10 A schematic diagram of iterative convergence of the objective function value and the volume fractions of three materials in the multi-material topology optimization process of the bridge structure under scenario 3 provided in embodiment 2 of the present invention;

[0057] Fig.11 The material distribution result of the three-dimensional bridge structure after multi-material topology optimization under scenario 3 provided in Example 2 of the present invention;

[0058] Fig.12 A comparison chart of the 3D bridge structure topology optimization results under scenario 3 provided in Example 2 of the present invention and the DMO method. DETAILED DESCRIPTION

[0059] The preferred embodiments of the present invention are described in detail below in conjunction with the accompanying drawings, wherein the accompanying drawings constitute a part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, but are not used to limit the scope of the present invention.

[0060] Example 1

[0061] A specific embodiment of the present invention discloses a multi-material topology optimization method based on a univariate discrete interpolation scheme, such as Figure 1 and Figure 2 As shown, the following steps are included:

[0062] S1. Discretize the design domain by finite element method and apply boundary conditions and external loads.

[0063] Specifically, finite element discretization is performed by meshing the design domain to obtain various finite element units.

[0064] Specifically, boundary conditions define the physical behavior at the boundaries of the model; external loads refer to external forces or influences applied to the design domain.

[0065] S2. Based on the univariate discrete combined interpolation method, a multi-material topology optimization model is constructed, and then the sensitivity of the multi-material topology optimization model is obtained.

[0066] During implementation, based on the univariate discrete combined interpolation method, a multi-material topology optimization model is constructed, including:

[0067] S21. Based on the univariate characteristic function, the design variable is mapped to the topological density field of the finite element unit. That is, based on the univariate characteristic function, a single design variable is mapped to multiple topological densities, and the topological density field is represented by different topological density functions.

[0068] In specific implementation, the univariate characteristic function is expressed as:

[0069]

[0070] In the formula, It means that when the design variable of the ξth finite element is χ(ξ), The topological density of the phase candidate material, It means that when the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of the phase candidate material, η represents the penalty function, Θ Indicates the number of phases of the candidate material.

[0071] Specifically, when the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of phase candidate materials It is expressed as:

[0072]

[0073] In the formula, Indicates The threshold parameter of the phase candidate material, β represents the projection parameter, tanh() represents the hyperbolic tangent function, and H() represents the Heaviside function. Among them, the projection parameter β controls the steepness of the Heaviside function.

[0074] Specifically, the design variable χ(ξ) in the ξ-th finite element unit represents a value between 0 and 1.

[0075] The proposed method constructs multiple characteristic functions through a single optimization variable, which are associated with multiple discrete values ​​to distinguish different material phases; the design variables are converted into topological density fields of multiple candidate materials using single variable characteristic functions.

[0076] S22. Based on the topological density field of the finite element unit, obtain the topological density field of the finite element unit after filtering, and then obtain the normalized topological density field of the finite element unit.

[0077] Specifically, the topological density field of the filtered finite element unit is obtained by:

[0078] First, the Helmholtz partial differential equation, expressed as:

[0079]

[0080] In the formula, represents the filter radius, ▽ represents the Laplace differential operator, It means that after filtering, when the design variable of the ξth finite element is χ(ξ), Topological density of phase candidate materials.

[0081] Then, based on formula (3), the following relationship is obtained after finite element discretization:

[0082]

[0083] In the formula, represents the discretized unfiltered density vector, Represented as the discretized filtered density vector, K ft represents the filter matrix, T ft Represents a transformation matrix.

[0084] Specifically, the filter matrix K ft and the transformation matrix T ft It is expressed as:

[0085]

[0086] in,

[0087]

[0088]

[0089] In the formula, Ω represents the integration domain of the finite element unit, Sym. indicates that the matrix is ​​symmetric along the diagonal, They represent the x, y, and z axis coordinates of the first unit node of the finite element unit, respectively. They represent the x, y, and z axis coordinates of the second node of the finite element unit, respectively. They represent the x, y, and z axis coordinates of the third node of the finite element unit, respectively. They respectively represent the x-, y-, and z-axis coordinates of the 4th element node of the finite element element.

[0090] It should be noted that a coordinate system, such as a Cartesian coordinate system, is established based on the characteristics of the design domain and the requirements of the optimization problem; the origin of the coordinate system can be selected as a specific point in the design domain, such as a symmetrical point or a mass point, and the direction of the coordinate axis can be determined according to the symmetry or other geometric characteristics of the design domain.

[0091] Specifically, the topological density of the filtered finite element unit is mapped to the interval (0,1) through the Heaviside projection function to obtain the normalized topological density field of the finite element unit, which is expressed as:

[0092]

[0093] In the formula, It means that after filtering, when the design variable of the ξth finite element is χ(ξ), Normalized topological density of phase candidate materials.

[0094] S23. Based on the normalized topological density field of the finite element unit, the elastic modulus of the candidate material is discretely interpolated to obtain the elastic modulus of the finite element unit, and then the elastic tensor of the finite element unit is obtained.

[0095] Specifically, the elastic modulus of the candidate materials is discretely interpolated based on the DMO method, and the elastic modulus and elastic tensor of each finite element unit are the weighted sum of the candidate materials, expressed as:

[0096]

[0097] Where, E[χ(ξ)] and D[χ(ξ)] represent the elastic modulus and elastic tensor of the ξ-th finite element when the design variable is χ(ξ); denote the elastic modulus and elastic tensor of the empty material, respectively. Respectively represent Elastic modulus and elastic tensor of the candidate materials, Indicates that in the ξth finite element The weight coefficient of the candidate material.

[0098] More specifically, Elastic modulus of candidate materials Given according to the specific application scenario.

[0099] More specifically, Elastic tensor of the phase candidate material Expressed as: According to Elastic modulus of candidate materials And Poisson's ratio calculation, specifically:

[0100] If it is a two-dimensional topology optimization, then

[0101]

[0102] If it is a three-dimensional topology optimization, then

[0103]

[0104] In the formula, It is set Poisson's ratio of the candidate material.

[0105] More specifically, to prevent singularity, the elastic modulus of the empty material Usually a very small value is set, such as According to The above formula is used to calculate.

[0106] More specifically, the ξth finite element Weight coefficient of phase candidate materials It is expressed as:

[0107]

[0108] In the formula, Φ k [χ(ξ)] respectively represent the Normalized topological density of k-phase candidate materials.

[0109] More specifically, the elastic modulus of the candidate material is an inherent material property of the candidate material, which is given according to an actual application scenario.

[0110] S24. Construct a multi-material topology optimization model based on the elastic tensor of the finite element unit.

[0111] Specifically, based on the elastic tensor of the finite element unit, the finite element analysis technology is used to establish the objective function and constraint function, and construct a multi-material topology optimization model.

[0112] Optionally, the multi-material topology optimization model is an optimization model with minimization of structural flexibility as the objective function and volume fraction of candidate materials as the constraint function, which is expressed as:

[0113]

[0114] Where n represents the total number of finite element units, C[u,χ(ξ)] represents the structural flexibility of the ξth finite element under the structural displacement field u when the design variable value is χ(ξ), ε(u) represents the linear strain tensor under the structural displacement field u, T represents the transpose, and a(u,δu) represents the linear strain tensor between the structural displacement field u and the Sobolev space H. 1 The double-line performance under the virtual displacement field δu in (Ω), l(δu) represents the double-line performance under the virtual displacement field δu in Sobolev space H 1 The linear load function under the virtual displacement field δu in (Ω), δu is the Sobolev space H 1 (Ω), |Ω| represents the volume of the design domain Ω, Γ Drepresents the displacement boundary, i.e., the Dirichlet boundary condition, and g represents the displacement boundary Γ D The displacement vector on Indicates Volume constraint function for phase candidate materials, Indicates The maximum volume fraction of the phase candidate material, υ0, represents the volume of each finite element.

[0115] Specifically, in the structural displacement field u and Sobolev space H 1 The double-line performance a(u, δu) under the virtual displacement field δu in (Ω) is expressed as:

[0116]

[0117] Where ε(δu) represents the linear strain tensor under the virtual displacement field δu.

[0118] Specifically, in the Sobolev space H 1 The linear load function l(δu) under the virtual displacement field δu in (Ω) is expressed as:

[0119]

[0120] In the formula, f represents the external force on the structure, and h represents the N The traction force distributed on N Represents the force boundary, that is, the Neumann boundary condition. Among them, the external force on the structure and the force on the boundary Γ N The upper distributed traction force is set according to the actual application scenario.

[0121] Correspondingly, in the optimization model with minimization of structural flexibility as the objective function and the volume fraction of candidate materials as the constraint function, the sensitivities of the objective function and constraint function are expressed as:

[0122]

[0123] in,

[0124]

[0125] In the formula, K ft represents the filter matrix, T ft represents the transformation matrix, sech() represents the hyperbolic secant function, Represents the matrix All elements of the column, It means that when the design variable value of the ξth finite element is χ(ξ), Topological density of phase candidate materials.

[0126] Optionally, the multi-material topology optimization model is an optimization model with minimization of structural flexibility as the objective function and mass density of candidate materials as the constraint function, which is expressed as:

[0127]

[0128] In the formula, G m represents the mass density constraint function, M f represents the maximum mass density, Indicates The mass density of the candidate materials of each phase is given according to the actual application scenario.

[0129] Correspondingly, the optimization model with minimization of structural flexibility as the objective function and the mass density of candidate materials as the constraint function, the sensitivities of the objective function and constraint function are expressed as:

[0130]

[0131] It should be noted that the sensitivity of the multi-material topology optimization model is derived as follows:

[0132] First, the derivative of the objective function with respect to the design variable χ(ξ) is derived by the chain rule:

[0133]

[0134] Among them, the objective function has a weight coefficient W ξ,j The sensitivity is expressed as:

[0135]

[0136] According to the differential law, we can get:

[0137]

[0138] Substituting formula (23) into formula (22), we get:

[0139]

[0140] The first-order derivative of the elastic tensor with respect to the weight coefficient is obtained by equation (8), and it is substituted into equation (24) to obtain:

[0141]

[0142] Based on formula (11), the weight coefficient pair is obtained Sensitivity:

[0143]

[0144] Based on formula (6), we get about The partial derivative of is expressed as:

[0145]

[0146] Combining the Helmholtz partial differential filtering equation given in equations (3) and (4), we can obtain the filtered topological density field: The topological density field before filtering The sensitivity is as follows:

[0147]

[0148] According to equations (1) and (2), the partial derivative of the single variable characteristic function with respect to the design variable is expressed as:

[0149]

[0150] Substituting equations (25), (26), (27), (28) and (29) into equation (21), we can obtain the sensitivity of the objective function to the design variables:

[0151]

[0152] in,

[0153]

[0154] In the formula, K ft represents the filter matrix, T ft represents the transformation matrix, sech() represents the hyperbolic secant function, It means that when the design variable value of the ξth finite element is χ(ξ), Topological density of phase candidate materials.

[0155] Then, the sensitivity of the constraint functions to the design variables is derived.

[0156] According to the chain rule, the volume constraint function The sensitivity to the design variables can be expressed as follows:

[0157]

[0158] The volume constraint function given by formula (12) can be obtained as right The partial derivatives of are as follows:

[0159]

[0160] Substituting equations (27), (28), (29) and (33) into equation (32), we can obtain the sensitivity of the volume constraint function to the design variables:

[0161]

[0162] Similarly, the quality constraint function G can be given m Sensitivity to design variable χ(ξ):

[0163]

[0164] S3. Based on the multi-material topology optimization model and sensitivity, the design variables are updated to obtain the optimal topological density field.

[0165] During implementation, based on the multi-material topology optimization model and sensitivity, the gradient optimization algorithm MMA is used to update the design variables until the multi-material topology optimization model converges, and the optimal design variables are obtained, and then the optimal topological density field is obtained, which is substituted into the post-processing function to draw the material distribution isosurface diagram.

[0166] In specific implementation, when the design variables are updated by the gradient optimization algorithm MMA, the optimization parameters in the iteration are updated based on the adaptive parameter adjustment strategy; wherein the optimization parameters include penalty parameters and projection parameters.

[0167] Specifically, the adaptive parameter adjustment strategy is:

[0168] Before the gradient optimization algorithm MMA iterates, the optimization parameters are initialized and the value of the counter a is initialized to 0;

[0169] During the MMA iteration of the gradient optimization algorithm, it is determined whether the projection parameter β′ in the current iteration is less than the set projection parameter threshold. If yes, the design variables in the current iteration are taken as the optimal design variables; otherwise:

[0170] Determine whether the value of the counter in the current iteration is greater than 0,

[0171] If the value of the counter a in the current iteration is less than or equal to 0, then determine whether the change δ between the objective function in the current iteration and the objective function in the previous iteration is less than the predetermined change value. If yes, the value a of the counter is set to 1, and the optimization parameters in the current iteration are not updated; if no, the optimization parameters in the current iteration are not updated;

[0172] If the value of the counter in the current iteration is greater than 0, the value a of the counter is increased by 1, and it is determined whether the value a of the counter at this time is the set number threshold. If so, the optimization parameters in the current iteration are updated and the value a of the counter is set to 1; if not, the optimization parameters in the current iteration are not updated.

[0173] More specifically, the projection parameters in the optimization parameters are updated by the following formula:

[0174] β=θ×β′ (36)

[0175] Where β′ is the projection parameter in the current iteration, β is the updated projection parameter, and θ represents the first adjustment coefficient.

[0176] Specifically, the penalty parameter in the optimization parameter is updated by the following formula:

[0177]

[0178] Where η′ is the penalty parameter in the current iteration, η is the updated penalty parameter, and σ represents the second adjustment coefficient.

[0179] More specifically, the predetermined value of the change amount is set to 1e-4; and the number threshold of the counter is set to 4.

[0180] More specifically, the first adjustment coefficient and the second adjustment coefficient are set according to actual requirements. The projection parameter threshold is set according to actual requirements.

[0181] More specifically, among the optimization parameters, the penalty parameter η is initialized to 1; the projection parameter β is initialized to a small positive value; illustratively, the projection parameter β=0.3.

[0182] Specifically, based on the obtained optimal design variables, the optimal topological density field of each phase candidate material is further obtained through formula (11).

[0183] It can be understood that the optimization method in this embodiment can adaptively adjust the penalty parameters and projection parameters, thereby making the algorithm more robust and adaptable during the iteration process, avoiding local optimal solutions and numerical instability; it can achieve reasonable distribution of materials in the structure under the premise of ensuring material mass or volume fraction constraints; in the optimization design, low-stiffness materials are mainly distributed in non-critical areas, medium-stiffness materials are distributed on the load transfer path, and high-stiffness materials are distributed in the load action area and the support area, and the load transfer path is clear; the use of intermediate materials can be avoided in some structures.

[0184] Compared with the prior art, the embodiment of the present invention provides a multi-material topology optimization method based on a single variable discrete interpolation scheme, which discretizes the design domain by finite elements, applies boundary conditions and external loads, and then constructs a multi-material topology optimization model based on a single variable discrete combined interpolation method, thereby obtaining the sensitivity of the multi-material topology optimization model, and then updates the design variables based on the multi-material topology optimization model and the sensitivity to obtain the optimal topological density field, thereby solving the problems of many design variables, low optimization efficiency, unreasonable design results, and difficult processing and manufacturing of the designed structure; by constructing a single variable characteristic function, using a single variable to select multiple discrete attribute physical quantities, combining the composite function of the single variable reduced-order multi-characteristic function with the corresponding material property for material interpolation, eliminating the material envelope phenomenon existing in the existing step-type interpolation model, and overcoming the disadvantage of limited design domain; making the number of design variables independent of the number of candidate material phases, significantly reducing the number of design variables, and achieving similar or even lower objective function values, the number of iterations is reduced, thereby improving the optimization efficiency; and being able to achieve reasonable distribution of materials in the structure under the premise of ensuring material quality or volume fraction constraints.

[0185] Example 2

[0186] In order to verify the effectiveness and correctness of the method proposed in Example 1, a specific example is provided for verification.

[0187] Scenario 1: Using volume constraints, material allocation is performed by setting the volume fraction of the candidate material in the design domain to minimize the structural flexibility and verify the effectiveness of the proposed method for two-dimensional topology optimization.

[0188] like Figure 3 As shown in Figure 1, taking a two-dimensional MBB beam as an example, the parameters are set as the length unit L = 60, the external force F = 1; half of the structure is selected as the design domain, and a 120×60 quadrilateral grid is used for structural discretization; the elastic parameters of the candidate materials and the given volume fraction constraints are shown in Table 1; the first adjustment coefficient of the projection parameter β is θ = 1.24, and the second adjustment coefficient of the penalty parameter η is σ = 1.065.

[0189] Table 1 Elastic parameters of candidate materials and given volume fraction constraints

[0190]

[0191] like Figure 4 The objective function value and volume fraction of the MBB beam structure are shown in the figure. It can be seen from the figure that the objective function value first rises and then drops rapidly at the beginning of the iteration, and stabilizes after about 50 iterations; the volume fraction curve shows that the volume fractions of MAT1 and MAT2 also gradually converge during the iteration process, and finally approach the set volume constraint of 0.25.

[0192] As shown in Figures 5(a) and 5(b), the two-dimensional MBB beam topology optimization results obtained by the method proposed in Example 1 under volume constraints are compared with DMO when θ = 1.24, σ = 1.065 and θ = 1.1445, σ = 1.0247, respectively. Among them, when θ = 1.24, σ = 1.065, the number of iterations in the DMO method is 89, the objective function is 23.9965, the volume fraction of material 1 is 0.2499, and the volume fraction of material 2 is 0.2498; when θ = 1.24, σ = 1.065, the number of iterations in the method of Example 1 is 88, the objective function is 23.7622, ​​the volume fraction of material 1 is 0.2497, and the volume fraction of material 2 is 0.2497. When θ=1.1445 and σ=1.0247, the number of iterations in the DMO method is 189, the objective function is 24.0376, the volume fraction of material 1 is 0.2500, and the volume fraction of material 2 is 0.2500; when θ=1.1445 and σ=1.0247, the number of iterations in the method of Example 1 is 191, the objective function is 23.5106, the volume fraction of material 1 is 0.2500, and the volume fraction of material 2 is 0.2500.

[0193] It can be seen from Figure 5(a) and Figure 5(b) that, under the premise of satisfying the volume fraction constraint, the method proposed in Example 1 and the DMO method achieve a lower objective function flexibility value when the number of iterations is close, indicating that its optimization effect is better than the DMO method. It should be noted that the material distribution obtained by the traditional DMO method is relatively mixed, which increases the difficulty of manufacturing. On the contrary, the result obtained by the method proposed in Example 1 has a concentrated material distribution, hard materials are allocated to areas with greater load-bearing areas, and soft materials are allocated to the middle area. The optimized structure is more suitable for force transmission and bearing. In addition, for the first adjustment coefficient and the second adjustment coefficient, the smaller the value, the better the objective function value will be, and a higher number of iterations will also be required.

[0194] Scenario 2: Using mass constraints, material allocation is performed by setting the mass of candidate materials to minimize structural flexibility and verify the effectiveness of the proposed method in two-dimensional topology optimization.

[0195] like Figure 6 As shown, a cantilever beam structure is taken as an example, where the length unit L=80, the design domain is discretized by a 120×80 quadrilateral unit with four nodes, and the middle part of the right section bears a unit load F.

[0196] Considering the structural multi-material topology optimization design of the cantilever beam under a single mass constraint, there are five materials in total, and the material parameters are shown in Table 2; the mass constraint is set to 0.3.

[0197] Table 2 Material parameters

[0198]

[0199] like Figure 7 As shown in the figure, the objective function and volume fraction change with the number of iterations under 5 material quantities. It can be seen that the objective function converges after 60 iterations. MAT1 and MAT2 also do not change after 60 iterations, while MAT3, MAT4 and MAT5 converge after 50 iterations. It should be noted that the volume fraction of MAT4 gradually decreases and finally converges to 0.

[0200] like Figure 8(a) and 8(b) The comparison diagram of the method proposed in Example 1 and the DMO method in the case of five materials (MAT1 / 2 / 3 / 4 / 5) and three materials (MAT2 / 3 / 5) is shown. Among them, when there are five materials (MAT1 / 2 / 3 / 4 / 5), the number of iterations of the DMO method is 82, the objective function is 39.2550, the volume fraction of material 1 is 0.1266, the volume fraction of material 2 is 0.3500, the volume fraction of material 3 is 0.1250, the volume fraction of material 4 is 0.0032, and the volume fraction of material 5 is 0.0031; the number of iterations of the method proposed in Example 1 is 86, the objective function is 36.8620, the volume fraction of material 1 is 0.5576, the volume fraction of material 2 is 0.1869, the volume fraction of material 3 is 0.0544, the volume fraction of material 4 is 0.0000, and the volume fraction of material 5 is 0.0037. When there are three materials (MAT2 / 3 / 5), the number of iterations of the DMO method is 87, the objective function is 39.3614, the volume fraction of material 1 is 0.2487, the volume fraction of material 2 is 0.2614, and the volume fraction of material 3 is 0.0056; the number of iterations of the method proposed in Example 1 is 86, the objective function is 38.2253, the volume fraction of material 1 is 0.5324, the volume fraction of material 2 is 0.0073, and the volume fraction of material 3 is 0.0289.

[0201] It can be seen that for the topology optimization results under five materials, the optimal objective function value obtained by the DMO method is 39.2550, while the optimal objective function value obtained by the proposed method is 36.8620, which is significantly better than the DMO method. In the final topology optimization configuration, material MAT4 is completely degraded, and the material MAT5 with the largest stiffness is selected to be filled in the load application area, while the MAT3 material is filled near the support end. It should be noted that the material distribution obtained by the DMO method is relatively scattered, while the material distribution obtained by the method in this paper is concentrated. For the topology optimization results under three materials, the minimum flexibility value of the DMO method is 39.3614, while the method proposed in Example 1 is 38.2253, which also shows better performance. In the material distribution, the method in Example 1 fills MAT5 at the load position and the support end position, and the material distribution is more reasonable.

[0202] Scenario 3: Using volume constraints, material allocation is performed by setting the volume fraction of the candidate material in the design domain to minimize structural flexibility and verify the effectiveness of the proposed method in three-dimensional topology optimization.

[0203] Taking the bridge design problem as an example, the structural tetrahedral mesh is generated using the open source 3D finite element mesh generator GMSH. The candidate materials are aluminum, iron, and titanium, and the volume fractions of the three materials are 0.2, 0.15, and 0.1, respectively. The material parameters are shown in Table 3.

[0204] Table 3 Elastic parameters of candidate materials for 3D bridge structures

[0205]

[0206] like Fig. 9 As shown in the figure, the loads and boundary conditions of the bridge design problem are given. The overall size of the bridge structure is 400mm×80mm×80mm. In this example, the bridge structure is divided into the design domain Ω D and non-design domain Ω ND . Design variables are allocated in the design domain for multi-material topology optimization, while no optimization variables are allocated in the non-design domain, and the material allocation and structural configuration in this area remain unchanged. Uniform force F = 16667Pa along ξ z The negative direction of the axis acts on the non-design area of ​​the bridge Ω ND The fixed displacement area Ω at both ends under the bridge B The boundary conditions are as follows. The number of meshes in the entire domain is 63463 and the number of nodes is 14267. Based on the proposed multi-material topology optimization method, the flexibility of the bridge structure is minimized while satisfying the volume constraint by reasonably allocating candidate materials at different spatial positions of the bridge. According to the symmetry of the structure, a quarter of the structure is taken during the optimization process. In post-processing, the topological configuration and material distribution of the complete structure are obtained by mirroring.

[0207] The objective function value of the bridge structure in the multi-material topology optimization process and the iterative convergence of the volume fractions of the three materials are shown in Figure 2. Fig.10 As shown. It can be seen that the objective function value first rises rapidly in the initial iteration and then drops rapidly after the 40th iteration, and then gradually stabilizes. In addition, the volume fraction of metal aluminum drops rapidly in the early iteration and then stabilizes. The other two materials show a trend of first dropping and then rising. After about 50 iterations, the volume fractions of the three materials basically converge. This shows that during the optimization process, the distribution of materials gradually reached a balance and achieved the minimization of the objective function, verifying the effectiveness and stability of the multi-material topology optimization method.

[0208] like Fig.11 As shown in the figure, the material distribution results of the three-dimensional bridge structure after multi-material topology optimization are given. Through the display of the topology optimization results from different perspectives in the figure, it can be seen that the high-rigidity metal material titanium is concentrated in the key force-bearing area and the supporting area. The medium-rigidity metal material iron is mainly distributed on the load transfer path, while the low-rigidity metal material aluminum is mainly distributed in non-critical areas and plays a connecting role. The obtained topological configuration and material distribution of the bridge structure not only improve the stiffness of the structure, but also enhance the load transfer path, and the resulting material distribution has clear and smooth boundaries. It can be considered that the multi-material topology optimization method proposed in this paper can effectively realize the reasonable distribution of materials in the three-dimensional structure to minimize the structural flexibility, verifying the feasibility of this method in the multi-material topology optimization design of complex three-dimensional structures.

[0209] Fig.12 A comparison chart of the results of 3D bridge structure topology optimization using the traditional DMO method and the method proposed in Example 1 is given. Among them, the number of iterations of the DMO method is 89, the objective function is 1.7675e-4, the volume fraction of material 1 is 0.1998, the volume fraction of material 2 is 0.1498, and the volume fraction of material 3 is 0.1000; the number of iterations of the method proposed in Example 1 is 79, the objective function is 1.4198, the volume fraction of material 1 is 0.2032, the volume fraction of material 2 is 0.1494, and the volume fraction of material 3 is 0.0999.

[0210] It can be seen that the results obtained by the method proposed in Example 1 are relatively consistent with the topological configuration obtained by the DMO method. In terms of material distribution, the two methods distribute high-rigidity titanium materials in the load-bearing area. The difference is that the method proposed in Example 1 distributes the material iron in the bridge top area, while the DMO method distributes aluminum materials. In addition, the method proposed in Example 1 distributes high-rigidity material titanium in the boundary fixed support area, while the DMO method also distributes aluminum materials. From the perspective of the load transfer path, the material distribution obtained by the method proposed in Example 1 is more favorable. It should be noted that the method proposed in Example 1 achieves a lower objective function value while using fewer iterations. This proves the superiority of the method proposed in Example 1.

[0211] Those skilled in the art will appreciate that all or part of the processes of the above-mentioned embodiments can be implemented by instructing related hardware through a computer program, and the program can be stored in a computer-readable storage medium, wherein the computer-readable storage medium is a disk, an optical disk, a read-only storage memory, or a random access memory, etc.

[0212] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by any technician familiar with the technical field within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention.

Claims

1. A multi-material topology optimization method based on a univariate discrete interpolation scheme, characterized in that: The following steps are involved: The design domain is discretized by finite element method and boundary conditions and external loads are applied; Based on the single variable discrete combined interpolation method, a multi-material topology optimization model is constructed, and then the sensitivity of the multi-material topology optimization model is obtained; Based on the multi-material topology optimization model and sensitivity, the design variables are updated to obtain the optimal topological density field.

2. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 1 is characterized in that: The multi-material topology optimization model is constructed based on the single variable discrete combination interpolation method, including: Based on the univariate characteristic function, the design variables are mapped into the topological density field of the finite element unit; Based on the topological density field of the finite element unit, a filtered topological density field of the finite element unit is obtained, and then a normalized topological density field of the finite element unit is obtained; Based on the normalized topological density field of the finite element unit, the elastic modulus of the candidate material is discretely interpolated to obtain the elastic modulus of the finite element unit, and then the elastic tensor of the finite element unit is obtained; Based on the elastic tensor of the finite element unit, a multi-material topology optimization model is constructed.

3. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 2 is characterized in that: The univariate characteristic function is expressed as: In the formula, It means that when the design variable of the ξth finite element is χ(ξ), The topological density of the phase candidate material, It means that when the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of the phase candidate material, η represents the penalty function, and Θ represents the number of phases of the candidate material.

4. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 3 is characterized in that: When the design variable of the ξth finite element is χ(ξ), Heaviside mapping function of phase candidate materials It is expressed as: In the formula, Indicates is the threshold parameter of the phase candidate material, β represents the projection parameter, and tanh( ) represents the hyperbolic tangent function.

5. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 4 is characterized in that: The elastic modulus and elastic tensor of the finite element unit are expressed as: Where, E[χ(ξ)] and D[χ(ξ)] represent the elastic modulus and elastic tensor of the ξ-th finite element when the design variable is χ(ξ); denote the elastic modulus and elastic tensor of the empty material, respectively. Respectively represent Elastic modulus and elastic tensor of the candidate materials, Indicates that in the ξth finite element The weight coefficient of the candidate material.

6. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 5 is characterized in that: The ξth finite element Weight coefficient of phase candidate materials It is expressed as: In the formula, Φ k [χ(ξ)] respectively represent the Normalized topological density of k-phase candidate materials.

7. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 1 is characterized in that: Based on the multi-material topology optimization model and sensitivity, the gradient optimization algorithm MMA is used to update the design variables until the multi-material topology optimization model converges to obtain the optimal design variables and then the optimal topological density field.

8. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 7 is characterized in that: When the design variables are updated by the gradient optimization algorithm MMA, the optimization parameters in the iteration are updated based on the adaptive parameter adjustment strategy; wherein the optimization parameters include penalty parameters and projection parameters.

9. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 7, characterized in that: The adaptive parameter adjustment strategy includes: Before the gradient optimization algorithm MMA iteration, the optimization parameters are initialized and the counter value is initialized to 0; During the MMA iteration of the gradient optimization algorithm, it is determined whether the projection parameter in the current iteration is less than the set projection parameter threshold. If so, the design variable in the current iteration is used as the optimal design variable; otherwise: Determine whether the value of the counter in the current iteration is greater than 0, If yes, determine whether the change between the objective function in the current iteration and the objective function in the previous iteration is less than the predetermined change value. If yes, set the value of the counter to 1, and do not update the optimization parameters in the current iteration; if no, do not update the optimization parameters in the current iteration; If not, the counter value is increased by 1, and it is determined whether the counter value at this time is the set number threshold. If so, the optimization parameters in the current iteration are updated and the counter value is set to 1; if not, the optimization parameters in the current iteration are not updated.

10. The multi-material topology optimization method based on a univariate discrete interpolation scheme according to claim 9 is characterized in that: The projection parameters in the optimization parameters are updated by the following formula: β=θ×β′ Where β′ is the projection parameter in the current iteration, β is the updated projection parameter, and θ represents the first adjustment coefficient. The penalty parameter in the optimization parameter is updated by the following formula: Where η′ is the penalty parameter in the current iteration, η is the updated penalty parameter, and σ represents the second adjustment coefficient.

Images