A multi-material structure topology optimization method

Abstract

The application discloses a multi-material structure topology optimization method, which innovatively applies Tsai-Wu criterion to multi-material structure topology optimization design, introduces a double-filter interface identification strategy, realizes accurate extraction of an interface phase field, and constructs a gradient interface finite element model based on a material proportion, then establishes an interface tension and compression state adaptive constraint method, and further adopts an adaptive threshold strategy to dynamically control the constraint tightness, so that full compression of the multi-material interface is realized at the algorithm level.

Description

Technical Field

[0001] This application relates to the field of structural optimization design, specifically to a topology optimization method for multi-material structures. Background Technology

[0002] Multimaterial structures can fully leverage the performance advantages of different materials, effectively improving overall performance while achieving structural lightweighting, thus attracting significant attention in modern engineering fields such as aerospace and automotive industries. With breakthroughs in multimaterial additive manufacturing technology, multimaterial topology optimization has become a key technology for realizing integrated innovative design of "materials-structures".

[0003] However, the interface, as a transition zone of mechanical properties between heterogeneous materials, is a potential weak point in load-bearing multimaterial structures, and its failure behavior directly restricts the overall reliability of the structure. Studies have shown that interfaces generally exhibit significant tensile-compressive asymmetry, meaning their compressive strength is much higher than their tensile strength, making them highly susceptible to debonding or cracking under tensile stress. Therefore, in the design of multimaterial structures, achieving a "fully compressed" or "compressive stress-dominated" layout at the interface is a key approach to avoid interface failure risks and improve the load-bearing reliability of the structure. Currently, optimization methods considering this characteristic can be mainly divided into two categories: one is an indirect method based on interface strain energy control, which guides the optimization process by introducing a penalty term related to the interface tensile energy into the design objective; the other is a direct method based on interface stress constraints, which explicitly controls the mechanical behavior of the interface by introducing an independent strength criterion. While the former offers flexibility, its optimization results are sensitive to weight parameters and it is difficult to strictly guarantee that the interface is under full pressure. The latter has clear physical meaning and more deterministic control, but its complex constitutive relationship leads to high computational costs. Other stress constraint methods that use simplified strength criteria are often limited by the theory itself and cannot fully realize the design of the interface under full pressure while ensuring numerical efficiency.

[0004] Therefore, developing an active interface stress constraint mechanism that combines numerical efficiency and control completeness within the existing topology optimization framework has become a key technological bottleneck that urgently needs to be overcome in achieving high-reliability multi-material structure design. However, existing technologies still have shortcomings in terms of theoretical completeness and control accuracy, and there is a lack of an efficient optimization design method that can keep the interface basically under pressure. Summary of the Invention

[0005] The purpose of this application is to improve the above-mentioned defects of the prior art to a certain extent, and to provide a multi-material structure topology optimization method, which is beneficial to realize full compression of multi-material interfaces in the optimization results.

[0006] In order to overcome the deficiencies of the prior art and achieve the purpose of this invention, the following technical solution is adopted to solve the problem: Scheme 1 involves a multi-material structure topology optimization method, which includes steps such as discretizing the structure within the design domain, establishing material subdomains, initializing the design variable set, and performing iterative design until the optimal design variable set is obtained. In the step of establishing material subdomains, a method is constructed... A material subdomain, The total number of materials is multi-material, and the design unit for each design domain corresponds to... Each material subdomain element corresponds to a design variable. ,in, Number the design unit. Number the corresponding material subdomain. =1…… All design variables Should meet In the iterative design process until the optimal set of design variables is obtained, each iterative design performs the following steps: S1: Based on the initial design variable set or the design variable set formed in the previous iteration, for Filtration is performed to obtain the unit filtration density. , forming a filter density field ; S2: Using the Heaviside function to filter cell density Mapping is performed to obtain the unit physical density. Forming a physical density field :

[0007] in, The parameter controlling the steepness of the Heaviside function is initially set to 2, and incremented every 25 iterations for a total of five increments, eventually reaching a value of 32. =0.5; S3: Obtain the element stiffness matrix of the design domain ; S4: Based on the element stiffness matrix of the design domain Assemble the global stiffness matrix according to the number. Then, according to the static equilibrium equation Calculate structural displacement ,in This is the equivalent load vector; S5.1: For unit physical density Filtration is performed to obtain the secondary unit filtration density. This forms a secondary filtration density field. ; S5.2: Using the Heaviside function to filter the density of the quadratic unit Mapping is performed to obtain the cell expansion density. , forming an expansion density field :

[0008] in, The initial value is 2, and it is incremented every 25 iterations, for a total of five increments, eventually reaching 32. =0.2; S5.3: Calculate the values ​​of two adjacent materials Interface phase field and interface density :

[0009]

[0010]

[0011] in, and Characterizing the two adjacent materials at the interface, Indicates the first Such materials, Indicates the first Such materials, =1…… , =1…… ; S5.4: Calculate the values ​​of two adjacent materials Element stress in the interface domain :

[0012]

[0013] in, for The elastic coefficient matrix of the material interface, For the first The elastic modulus of the material, For the first The elastic modulus of the material, For element displacement vectors, For geometric matrices, , , , , , The stress components of the interface domain element stress. , , The global rectangular coordinate system in which the interface domain unit is located ( coordinate axes; S5.5: Obtain the interfacial tensile / compressive coefficient :

[0014]

[0015]

[0016]

[0017]

[0018]

[0019] Among them, and These represent the tensile and compressive strengths of the interface element along axis 1, respectively. and These represent the tensile and compressive strengths of the interface element in the axial direction 2, respectively. and These represent the tensile and compressive strengths of the interface element along axis 3, respectively. , , These represent the shear strengths in the yz, zx, and xy planes, respectively, with axes 1, 2, and 3 forming the interface coordinate system. The three-dimensional axes of the interface coordinate system and the global rectangular coordinate system are transformed by two rotations around the axes. First, the rotations around the global coordinate system are transformed by two rotations around the axes. Axis rotation The result after rotating around Axis rotation The sign of the rotation angle follows the right-hand coordinate system rule, and the corresponding relationship of the rotation is as follows: , , The expressions for the direction cosines between each axis are shown in the table below:

[0020] S5.6: Use the Heaviside function to obtain the interface tensile / compressive state coefficients. :

[0021] in, The parameters that control how closely the Heaviside function approximates 0 and 1 are... Left side of the middle The natural constant is represented by the subscript on the right. Number the design units; S5.7: Establish interface tension / compression state constraints:

[0022]

[0023] in, Represents the entire interface field. It is the first design domain Volume of each unit It is the first The maximum allowable structural volume fraction of a material subdomain. Refers to the first Material subdomains of a type of material Indicates the constraint threshold; S6: Construct a mathematical model that minimizes structural flexibility while considering interface and volume constraints:

[0024]

[0025]

[0026]

[0027]

[0028] in, It is the structural compliance function. It is the interface tension / compression state adaptive constraint method function. It is the first Volume constraints on the material subdomains of a given material; S7: Perform sensitivity analysis; S8: The moving asymptote algorithm is used to solve the optimization problem, update the design variables of each unit, and obtain and store the set of design variables formed in this iteration; S9: Determine whether the current iteration meets the exit condition. If the current iteration meets the exit condition, exit the iteration and record the set of design variables formed in this iteration as the optimal design variable set. If the exit condition is not met, proceed to the next iteration. In the above iterative design steps, each step can start execution after obtaining the required input.

[0029] Option 2 is based on Option 1, but includes S5.8 between S5.7 and S6: Constructing adaptive constraints for interface tension and compression states. In the first 30 iterations, The value is fixed at the initial value. Then, every 10 iterations, it is updated in the following way:

[0030] in, and This represents the constraint threshold that has been updated twice consecutively. When =0, The value of is , for In the interface domain The integral value on the surface represents the integral value of the tension element. for In the interface domain The integral value on the surface represents the integral value of the compressed element.

[0031] Option 3 is based on Option 2, wherein, in S7, Relative to unit physical density The derivative expression is as follows:

[0032] in,

[0033]

[0034]

[0035] in, For geometric matrices,

[0036]

[0037] The interaxial direction cosines with respect to the element physical density The derivative is:

[0038]

[0039]

[0040]

[0041] The interface gradient is related to the unit physical density. The derivative is

[0042]

[0043] According to the chain rule, we can obtain Relative to unit design variables Derivative.

[0044] Option 4 is based on Option 3, wherein in S5.5:

[0045]

[0046]

[0047]

[0048] in, Strength reduction factor ( <1).

[0049] Option 5 is based on Option 4, where, in S3,

[0050]

[0051]

[0052] in, The material penalty factor has a value of 4. Indicates except the first Other than this type of material Such materials, For the first The elastic modulus of the material, For the first The unit stiffness matrix of a material.

[0053] Scheme 6 is based on Scheme 4. In S9, the method to determine whether the current iteration meets the exit conditions is to compare the design variable set formed in the current iteration with the design variable set formed in the previous iteration, or to compare the objective function value obtained in the current iteration with the objective function value obtained in the previous iteration. In addition, the interface tension and compression constraints and volume constraints are both satisfied, or the total number of iterations reaches 400.

[0054] Scheme 7 is based on Scheme 6, where the exit condition for iteration is that the absolute value of the difference between all design variables in the design variable set formed in this iteration and the corresponding design variables in the design variable set formed in the previous iteration is less than 0.001, or the total number of iterations reaches 300.

[0055] Option 8 is based on Option 1, where, in S6:

[0056] in, Set to 10 -6 .

[0057] As can be seen from the above description of this application, compared with the prior art, this application has the following beneficial effects: The Tsai-Wu criterion, a classic strength criterion, has long been used for failure assessment of composite materials. However, due to the lack of equivalent stress modeling and conversion methods within the framework of multi-material, continuum topology optimization, it cannot be directly applied to the topology optimization design of multi-material structures. This application introduces a quadratic filtered density field to accurately identify the complex interfaces between multi-material regions. Then, based on the interface density gradient and combined with the stress rotation axis formula, it innovatively derives the "interface equivalent stress" applicable to discrete elements in topology optimization, thereby extracting the input with clear physical meaning required by the Tsai-Wu criterion. It is precisely through the pioneering "interface equivalent stress" modeling method of this application that the Tsai-Wu criterion can be applied to the topology optimization of multi-material structures: This application converts the discrete, heterogeneous material distribution and macroscopic strain field in topology optimization into the equivalent stress state acting on the interface that can be judged by the Tsai-Wu criterion. Then, the Tsai-Wu criterion's judgment result on the interface strength (compression / tension) is transformed into the constraint conditions and sensitivity information in the iterative calculation of structural topology optimization. The introduction of the Tsai-Wu criterion makes this method advantageous for achieving full compression of multi-material interfaces in the optimization results.

[0058] Building upon the foregoing, this application further proposes a constraint threshold. The adaptive constraint method makes the constraint threshold During the optimization process, the constraint values ​​are dynamically adjusted according to the stress state of the interface elements. The constraints are gradually tightened during the optimization process to effectively eliminate the tension interface. At the same time, it avoids the structural design from being too conservative or the iterative calculation from being difficult to converge due to excessive constraints. This makes the topology optimization results of multi-material structures better, and then realizes full compression of multi-material interfaces at the algorithm level. Attached Figure Description

[0059] To more clearly illustrate the technical solutions of the embodiments of this application, the accompanying drawings used are briefly described below.

[0060] Figure 1 This is a flowchart of the structural topology optimization method in this embodiment.

[0061] Figure 2 This is a schematic diagram of the design domain in this embodiment.

[0062] Figure 3 This is a schematic diagram of the global coordinate axis and the interface coordinate axis in this embodiment.

[0063] Figure 4 This is a schematic diagram illustrating the optimization results of the existing method.

[0064] Figure 5 This is a schematic diagram of the optimization results of the structural topology optimization method in this embodiment. Detailed Implementation

[0065] Unless otherwise expressly defined, the use of terms such as "first," "second," or "third" in the claims, description, and accompanying drawings of this application is for the purpose of distinguishing different objects, not for describing a specific order.

[0066] Unless otherwise specified, in the claims and description, the terms "comprising," "having," and variations thereof mean "including but not limited to."

[0067] In the claims and description, unless otherwise specified, the term "have" means that a technical feature that follows is part of a technical feature that precedes it.

[0068] Unless otherwise specified, the terms “fixed connection” or “relatively fixed” in the claims and description shall be interpreted broadly to mean any connection in which there is no displacement or relative rotation relationship between the two parties, including non-removable fixed connection, detachable fixed connection, fixed connection by fasteners, integral connection, direct fixed connection or indirect fixed connection.

[0069] Unless otherwise specified, in the claims and description, the terms “upper,” “lower,” “front,” “rear,” “left,” “right,” etc., indicate directions or positional relationships based on the directions explicitly shown in the drawings.

[0070] Unless otherwise specified in the claims and description, the terms "up and down direction", "front and back direction", and "left and right direction" essentially mean that the up and down direction, the front and back direction, and the left and right direction are relatively perpendicular to each other, that is, the up and down direction is perpendicular to the front and back direction, the up and down direction is perpendicular to the left and right direction, and the front and back direction is perpendicular to the left and right direction.

[0071] Unless otherwise specified in the claims and description, "When not used as a subscript, it represents a natural constant."

[0072] Unless otherwise specified in the claims and description, "Indicates two adjacent materials" The interface phase field, " represents the entire interface field, " "Referring to the first Material subdomains of a type of material and The intersection is All the interface phase fields are combined into .

[0073] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings.

[0074] Example This embodiment provides a method for topology optimization of structures under full compression at multi-material interfaces, and its flowchart is as follows: Figure 1 As shown.

[0075] like Figure 2 As shown, the design objective in this embodiment is to minimize the flexibility of the three-dimensional cantilever beam structure and ensure full-interface compression. The design domain is a three-dimensional cantilever beam of 80 mm × 40 mm × 20 mm. The constraint condition is that a uniformly distributed downward load with a resultant force F = 500 kN acts within a 5 mm × 5 mm area (red area in the figure) at the center of the right side of the beam. The left side is a fully fixed boundary. The upper limit of the volume fraction of each material sub-region is set to 10% of the total volume of the design domain. Two materials are selected (M = 2, The total number of materials (M1 in the attached diagram refers to the first material, and M2 refers to the second material) is shown in the table below:

[0076] The specific optimization method is as follows: First, the structure within the design domain is discretized. The design domain is discretized using 1mm linear hexahedral elements. Design unit.

[0077] Secondly, establish the material subdomain. Construct... Each material subdomain includes a corresponding design unit. Each material subdomain unit corresponds to a design unit. Each material subdomain element. Two materials are used in this example, therefore… =2, therefore one design unit corresponds to two material subdomains. Each material subdomain corresponds to one design variable. ,in When used as a subscript, it represents the design unit number. =1…… , Number the corresponding material subdomain. =1…… Therefore Read as the first In the material subdomain, the first The design variables of each element. The properties of the design domain elements are determined by their corresponding... The interpolation of design variables for each material subdomain element is obtained, for example, when =1, =0, which indicates that the material selected for this design unit is the first type of material; when =0, =1, which indicates that the material selected for this design unit is the second type of material; when =0, =0 indicates that the design unit is not filled with material; does not exist. =1, The case where =1. Therefore, when the material subdomain element design variables... Once the values ​​of are determined, the structure of the design domain is also determined, meaning that the design variable set for the topology optimization method in this application is . A collection of [items / items].

[0078] Furthermore, initialize the design variable set. The initial values ​​are all set to 0.5.

[0079] Finally, iterative design is performed until the optimal set of design variables is obtained. The specific process of iterative design is as follows: Step S1: Based on the initial design variable set or the design variable set formed in the previous iteration, adjust the design variables... Filtration is performed to obtain the unit filtration density. , forming a filter density field The filter radius is set to 3.5 unit dimensions, which is 3.5 mm in this example. The filtration method can adopt existing technologies, such as Bruns TE, Tortorelli DA. Topology optimization of non-linear elastic structures and compliant mechanisms. Comput Methods Appl Mech Eng 2001;190(26–27):3443–59.

[0080] Step S2: Use the Heaviside function to filter the cell density. Mapping is performed to obtain the unit physical density. Forming a physical density field Each material subdomain corresponds to a physical density field, and subsequent density fields follow the same principle. The expression for the Heaviside function in this step is as follows:

[0081] in, The parameter controlling the steepness of the Heaviside function is initially set to 2, and incremented every 25 iterations for a total of five increments, eventually reaching a value of 32. =0.5.

[0082] Step S3: Obtain the element stiffness matrix of the design domain The details are as follows:

[0083]

[0084]

[0085] in, The material penalty factor has a value of 4. Indicates except the first Other than this type of material Such materials, For the first The elastic modulus of this material can be obtained from a material handbook. For the first The unit stiffness matrix of the material is given by the existing formula, which can be found in textbooks (Zeng Pan. Finite Element Analysis and Applications [M]. Tsinghua University Press, 2004).

[0086] Step S4: Perform finite element analysis. Based on the element stiffness matrix of the design domain. Assemble the global stiffness matrix according to the number. Then, according to the static equilibrium equation Calculate structural displacement ,in This is the equivalent load vector.

[0087] Step S5.1: Assess the physical density of the unit cell. Filtration is performed to obtain the secondary unit filtration density. This forms a secondary filtration density field. The filter radius is 6mm, and the filtration method is the same as in step S1.

[0088] Step S5.2: Use the Heaviside function to filter the density of the secondary unit. Mapping is performed to obtain the cell expansion density. , forming an expansion density field The Heaviside function used in this step is as follows:

[0089] in, The initial value is 2, and it is incremented every 25 iterations, for a total of five increments, eventually reaching 32. =0.2.

[0090] Step S5.3: Calculate the values ​​of two adjacent materials Interface phase field and interface density The details are as follows: use and Characterizing the two adjacent materials at the interface, Indicates the first Such materials, Indicates the first Such materials, =1…… , =1…… , The subscript int indicates the interface. According to step S5.2, the unit expansion density is obtained. and expansion density field After that, you can obtain the first Unit expansion density of a material subdomain and expansion density field and the Unit expansion density of a material subdomain and expansion density field Then we have:

[0091]

[0092]

[0093] Where, assuming This is to avoid interface phase field Repeated calculations.

[0094] Step S5.4: Calculate the values ​​of two adjacent materials Element stress in the interface domain The details are as follows: In practice, interfaces often contain two materials simultaneously, and their mechanical properties are determined by the mixing ratio of these two materials. Therefore, the elastic coefficient matrix of the interface element (the design unit within the interface domain) can be determined based on the material ratio. It is assumed that the elastic coefficient matrix of the interface element is obtained by linear interpolation of the elastic coefficients of adjacent materials. Elastic coefficient matrix of material interface It can be represented as:

[0095] in, For the first The elastic modulus of the material, For the first The elastic modulus of the material, assumed here. The value of 0.01 is to ensure that the interface domain unit still has a large elastic coefficient matrix when the physical density of both materials is very small. and The value of the unit physical density has been obtained in step S2. Obtained at that time. Element stress in material interface domain It can be represented as:

[0096] in, The element displacement vectors are obtained in step S4, where the structural displacements have already been acquired. According to the unit number The corresponding element displacement vector can be extracted. , For the geometric matrix (see: Zeng Pan. Finite Element Analysis and Applications [M]. Tsinghua University Press, 2004.). , , , , , The stress components of the interface domain element stress. , , The global rectangular coordinate system in which the interface domain unit is located ( Coordinate axes, such as Figure 3 As shown.

[0097] Step S5.5: Obtain the interfacial tensile and compressive coefficients based on the Tsai-Wu criterion. .

[0098] The Tsai-Wu criterion assesses failure risk by comprehensively evaluating the stress state of material points in different directions, and this assessment depends on the coordinate system referenced by the stress. However, in multi-material topologies, interfaces have arbitrary spatial orientations, making the Tsai-Wu criterion unsuitable for direct application. The inventive point of this application lies in solving this problem, enabling the Tsai-Wu criterion to be applied to topology optimization of multi-material structures, as detailed below: First, calculate the rotation matrix between the global coordinate system and the interface coordinate system. .

[0099] The relationship between the global coordinate system and the interface coordinate system is as follows: Figure 3 As shown, the interface coordinate system ( In the coordinate system, axis 3 is along the interface normal, while axes 1 and 2 lie within the interface tangent plane. The transformation between them is achieved through two rotations around the axes: first around the global coordinate system... Axis rotation The result after rotating around Axis rotation The sign of the rotation angle follows the rules of the right-hand coordinate system.

[0100] The correspondence of rotations is as follows: , , The expressions for the direction cosines between each axis are shown in the table below:

[0101] The corresponding coordinate axis transformation relationship can be expressed as:

[0102] in, This represents the coordinate values ​​on their respective coordinate axes. =0, =0, When =1, then:

[0103]

[0104]

[0105]

[0106] in, , , These represent the coordinates of the unit vector normal to the interface. , , It can be obtained through the following formula:

[0107] in, For the first The secondary filtration density field of the material subdomain of the material has been obtained in step S5.1. Obtained at that time.

[0108] After obtaining the rotation matrix, the transformation relationship between the stress of the interface domain elements in the global coordinate system and the stress of the interface domain elements in the interface coordinate system is obtained based on the rotation matrix.

[0109]

[0110] in, The possible values ​​are shown in the table above. In coordinate system ( (under) Element stress in the material interface domain.

[0111] Finally, the interfacial tensile and compressive coefficients were calculated using the Tsai-Wu criterion. Determine the tension / compression state of the interface.

[0112]

[0113] in, There are twelve intensity tensor coefficients, whose values ​​can be obtained using the following formula:

[0114]

[0115]

[0116]

[0117] Among them, and These represent the tensile and compressive strengths of the interface element along axis 1, respectively. and These represent the tensile and compressive strengths of the interface element in the axial direction 2, respectively. and These represent the tensile and compressive strengths of the interface element along axis 3, respectively. , , These are the shear strengths in the yz, zx, and xy planes, respectively.

[0118] Given that the core objective of this step is to assess the tensile-compressive asymmetry of the interface rather than to accurately characterize the anisotropic strength of the material, the shear strength parameter is assigned based on widely adopted engineering empirical methods:

[0119]

[0120]

[0121] To simulate the mechanical properties where the normal compressive strength of the interface is much higher than its tensile strength, this study uses a set of virtual strength parameters for characterization:

[0122] in, Strength reduction factor ( <1), used to describe the characteristic that the tensile strength of the interface normal (axis 3) is significantly lower than the compressive strength, and is set to 0.05.

[0123] In this application, , , , , , , , and The settings are as follows: =850MPa =42.5MPa =24.5MPa =490.5MPa Step S5.6: Use the Heaviside function to obtain the interface tension / compression state coefficients. .

[0124] Due to the interfacial tensile and compressive coefficients obtained based on the Tsai-Wu criterion The tensile and compressive states of an interface are determined by positive and negative values. Positive values ​​indicate that the interface element is in a tensile state, and negative values ​​indicate that the interface element is in a compressive state. However, in multi-material structure topology optimization, the tensile and compressive coefficients of each interface need to be summed, which can cause positive and negative values ​​to cancel each other out. Therefore, the tensile and compressive coefficients of the interfaces... It cannot be directly applied to topology optimization of multi-material structures. This step introduces the Heaviside function to... Mapping is performed to obtain the interface tension / compression state coefficient. This is to facilitate its application in structural topology optimization. Furthermore, since the algorithm in this application is a gradient algorithm, the mapping function in this step cannot be a step-type Heaviside function. The specific calculation process is as follows:

[0125] In this step, the Heaviside function has a range of [0,1] and is always non-negative. Since it is a non-step function, it has a transition interval; its range is not strictly 0 or 1 but includes continuous values ​​between 0 and 1. Values ​​above 0.5 indicate that the interface element is under tension, values ​​below 0.5 indicate that the interface element is under compression, and a value of 0.5 indicates that the interface element is not under stress (treated as compression in this method). This allows us to characterize the tension and compression state of each interface element, thus enabling its application in subsequent structural topology optimization. The parameter used to control how closely the Heaviside function approximates 0 and 1 is set to 8 in this example. Left side of the middle The natural constant is represented by the subscript on the right. Number the design unit.

[0126] Step S5.7: Establish interface tension / compression state constraints. Details are as follows: because This represents the tension / compression state of a single unit within the interface domain; therefore, in this step, the entire interface domain is considered. Upper interface tensile and compressive state coefficient Perform integration to construct global interface pull and pressure indicators. ,

[0127] Global interface tension and compression indicators This can characterize the number of tensile elements, although it was already defined in step S5.6. Values ​​≤ 0.5 are considered under pressure, but since this value is also > 0, the global interface tension / compression index is calculated after integration. He made a contribution, therefore The value will be slightly larger than the actual number of tension elements, so if the absolute value of 0 is used directly... Imposing constraints can lead to difficulties in algorithm convergence, thus requiring the introduction of a non-zero constraint. To reduce the sensitivity of this constraint to specific problems and improve its applicability across different problems, this application... A normalized constraint form is proposed. :

[0128] in, It is the first design domain The volume of each unit is determined when the unit is divided. It is the first The maximum allowable structural volume fraction for a material subdomain is set to 10%. Refers to the first Material subdomains of a type of material This represents the constraint threshold. This expression improves the adaptability of the parameters to problems with different scales and numbers of interfaces by normalizing the integral value to the structural target volume. This can serve as a reference value for other issues. Within this framework, The physical interpretation is approximately the volume fraction of the design domain occupied by the tension interface element.

[0129] Step S5.8: Construct adaptive constraints for interface tension and compression states.

[0130] To determine the appropriate To ensure that all interface units are under pressure while avoiding algorithm convergence difficulties due to excessively small values, this application proposes a method... Adaptive threshold selection strategy. In the first 30 iterations of optimization, The value is fixed at the initial value. The value is 0.05, and it is updated every 10 iterations as follows:

[0131] in, and This represents the constraint threshold that has been updated twice consecutively. When =0, The value of is . for In the interface domain The integral value on the surface represents the integral value of the tension element. for In the interface domain The integral value on the figure represents the integral value of the compressed element. , This indicates that the integral value of the compressed unit is within the total integral value (global interface tension / compression index). The percentage within ) . Ideally, the entire interface is under pressure, The value is This refers to the volume fraction of compressed elements in the structure that are mapped to positive values. Constraint threshold. This adaptive strategy can gradually tighten constraints and effectively eliminate tension interfaces during optimization, while avoiding conservative structural design due to excessive constraints. When all interfaces approach a compressed state, the constraint threshold... This will tend to stabilize, thus promoting algorithm convergence.

[0132] In the above steps, steps S5.1 to S5.3 can be executed after step S2, and step S5.4 can be executed after steps S4 and S5.3. Essentially, each step can begin execution after obtaining its required input.

[0133] Step S6: Construct a mathematical model that minimizes structural flexibility while considering interface constraints and volume constraints.

[0134]

[0135]

[0136]

[0137]

[0138]

[0139]

[0140] in, It is the structural compliance function. It is the interface tension / compression state adaptive constraint method function. It is the first Volume constraints of material subdomains of a given material. Set to 10 -6 This is used to prevent matrix singularities during numerical computation.

[0141] Step S7: Perform sensitivity analysis.

[0142] In this step, besides... Sensitivity analyses other than those described above are all existing techniques. Therefore, the following only introduces sensitivity analyses for... The sensitivity analysis process.

[0143] Relative to unit physical density The derivative expression is as follows:

[0144] in,

[0145]

[0146]

[0147] in, For geometric matrices,

[0148]

[0149] The interaxial direction cosines with respect to the element physical density The derivative is:

[0150]

[0151]

[0152]

[0153] The interface gradient is related to the unit physical density. The derivative is

[0154]

[0155] Finally, according to the chain rule, we can obtain... Relative to unit design variables Derivative.

[0156] Step S8: Use the moving asymptote algorithm to solve the optimization problem, update the design variables of each unit, and obtain and store the set of design variables formed in this iteration.

[0157] Step S9: Determine if the current iteration meets the exit condition. If the current iteration meets the exit condition, it exits the iteration, and the set of design variables formed in this iteration is recorded as the optimal design variable set; if the exit condition is not met, the next iteration begins. Specifically, the method for determining whether the current iteration meets the exit condition is to compare the set of design variables formed in this iteration with the set of design variables formed in the previous iteration, or to compare the objective function value obtained in this iteration with the objective function value obtained in the previous iteration, and both the interface tension / compression constraints and volume constraints are satisfied. For example, the exit condition is that the absolute value of the difference between all design variables in the design variable set formed in this iteration and the corresponding design variables in the design variable set formed in the previous iteration is less than 0.05. In this embodiment, the optimization convergence criterion is that the maximum change in design variables is less than 0.001, or the total number of iterations reaches 300.

[0158] The optimal structure can then be obtained based on the optimal set of design variables. The effectiveness of this method is verified below.

[0159] Figure 4 This paper demonstrates the optimization configuration of a three-dimensional cantilever beam using existing methods (Yoon GH, Park YC, Kim YY. Element stacking method for topology optimization with material-dependent boundary and loading conditions. Mechanics of Materials and Structures 2007;2:883-95.), as well as the interfacial tensile and compressive coefficients from different perspectives. and interface tensile and compressive state coefficient The distribution Figure 5 The corresponding optimization results of this method are displayed. The comparison results of key performance indicators for the two optimization schemes are shown in the table below.

[0160]

[0161] from Figure 4 and Figure 5 It can be seen from this: ① In existing methods, the material interface portion is located in the tension region. Positive values ​​exist; however, in this method, the material interface clearly migrates towards the compression region, indicating that this method can reduce the generation of tension interfaces.

[0162] ②In existing methods, The positive and negative value regions are symmetrical, indicating the simultaneous existence of interfaces under compression and tension; however, in this method, the optimized configuration exhibits asymmetry, and the global domain... All values ​​are less than 0, indicating that the material interfaces are all under pressure, proving that this method can achieve full pressure on the material interfaces of three-dimensional multi-material structures.

[0163] As can be seen from the table above, the structural compliance obtained by this method is 5.91. Slightly higher than the existing method's 5.85. This demonstrates that the proposed method achieves full compression at multi-material interfaces with minimal impact on structural stiffness. Therefore, although errors may exist during actual product manufacturing, the algorithmic approach effectively achieves full compression at multi-material interfaces.

[0164] The description of the above specification and embodiments is used to explain the scope of protection of this application, but does not constitute a limitation on the scope of protection of this application.

Claims

1. A topology optimization method for multi-material structures, characterized in that, This includes steps such as discretizing the structure within the design domain, establishing material subdomains, initializing the design variable set, and iterative design until the optimal design variable set is obtained. In the step of establishing material subdomains, the following steps are constructed: A material subdomain, The total number of materials is multi-material, and the design unit for each design domain corresponds to... Each material subdomain element corresponds to a design variable. ,in, Number the design unit. Number the corresponding material subdomain. =1…… All design variables Should meet In the iterative design process until the optimal set of design variables is obtained, each iterative design performs the following steps: S1: Based on the initial design variable set or the design variable set formed in the previous iteration, for Filtration is performed to obtain the unit filtration density. , forming a filter density field ; S2: Using the Heaviside function to filter cell density Mapping is performed to obtain the unit physical density. Forming a physical density field : in, The parameter controlling the steepness of the Heaviside function is initially set to 2, and incremented every 25 iterations for a total of five increments, eventually reaching a value of 32. =0.5; S3: Obtain the element stiffness matrix of the design domain ; S4: Based on the element stiffness matrix of the design domain Assemble the global stiffness matrix according to the number. Then, according to the static equilibrium equation Calculate structural displacement ,in This is the equivalent load vector; S5.1: For unit physical density Filtration is performed to obtain the secondary unit filtration density. This forms a secondary filtration density field. ; S5.2: Using the Heaviside function to filter the density of the quadratic unit Mapping is performed to obtain the cell expansion density. , forming an expansion density field : in, The initial value is 2, and it is incremented every 25 iterations, for a total of five increments, eventually reaching 32. =0.2; S5.3: Calculate the values ​​of two adjacent materials Interface phase field and interface density : in, and Characterizing the two adjacent materials at the interface, Indicates the first Such materials, Indicates the first Such materials, =1…… , =1…… ; S5.4: Calculate the values ​​of two adjacent materials Element stress in the interface domain : in, for The elastic coefficient matrix of the material interface, For the first The elastic modulus of the material, For the first The elastic modulus of the material, For element displacement vectors, For geometric matrices, , , , , , The stress components of the interface domain element stress. , , The global rectangular coordinate system in which the interface domain unit is located ( coordinate axes; S5.5: Obtain the interfacial tensile / compressive coefficient : Among them, and These represent the tensile and compressive strengths of the interface element along axis 1, respectively. and These represent the tensile and compressive strengths of the interface element in the axial direction 2, respectively. and These represent the tensile and compressive strengths of the interface element along axis 3, respectively. , , These represent the shear strengths in the yz, zx, and xy planes, respectively, with axes 1, 2, and 3 forming the interface coordinate system. The three-dimensional axes of the interface coordinate system and the global rectangular coordinate system are transformed by two rotations around the axes. First, the rotations around the global coordinate system are transformed by two rotations around the axes. Axis rotation The result after rotating around Axis rotation The sign of the rotation angle follows the right-hand coordinate system rule, and the corresponding relationship of the rotation is as follows: , , The expressions for the direction cosines between each axis are shown in the table below: S5.6: Use the Heaviside function to obtain the interface tensile / compressive state coefficients. : in, The parameters that control how closely the Heaviside function approximates 0 and 1 are... Left side of the middle The natural constant is represented by the subscript on the right. Number the design units; S5.7: Establish interface tension / compression state constraints: in, Represents the entire interface field. It is the first design domain Volume of each unit It is the first The maximum allowable structural volume fraction of a material subdomain. Refers to the first Material subdomains of a type of material Indicates the constraint threshold; S6: Construct a mathematical model that minimizes structural flexibility while considering interface and volume constraints: in, It is the structural compliance function. It is the interface tension / compression state adaptive constraint method function. It is the first Volume constraints on the material subdomains of a given material; S7: Perform sensitivity analysis; S8: The moving asymptote algorithm is used to solve the optimization problem, update the design variables of each unit, and obtain and store the set of design variables formed in this iteration; S9: Determine whether the current iteration meets the exit condition. If the current iteration meets the exit condition, exit the iteration and record the set of design variables formed in this iteration as the optimal design variable set; if the exit condition is not met, proceed to the next iteration. In the above iterative design steps, each step can begin execution after obtaining the required inputs.

2. The multi-material structure topology optimization method as described in claim 1, characterized in that, Between S5.7 and S6, there is also S5.8: Constructing adaptive constraints for interface tension and compression states: In the first 30 iterations, The value is fixed at the initial value. Then, every 10 iterations, it is updated in the following way: in, and This represents the constraint threshold that has been updated twice consecutively. When =0, The value of is , for In the interface domain The integral value on the surface represents the integral value of the tension element. for In the interface domain The integral value on the surface represents the integral value of the compressed element.

3. The multi-material structure topology optimization method as described in claim 2, characterized in that, In S7 Relative to unit physical density The derivative expression is as follows: in, in, For geometric matrices, The interaxial direction cosines with respect to the element physical density The derivative is: The interface gradient is related to the unit physical density. The derivative is According to the chain rule, we can obtain Relative to unit design variables Derivative.

4. The multi-material structure topology optimization method as described in claim 3, characterized in that, In S5.5: in, Strength reduction factor ( <1).

5. The multi-material structure topology optimization method as described in claim 4, characterized in that, In S3, in, The material penalty factor has a value of 4. Indicates except the first Other than this type of material Such materials, For the first The elastic modulus of the material, For the first The unit stiffness matrix of a material.

6. The multi-material structure topology optimization method as described in claim 4, characterized in that, In S9, the method to determine whether the current iteration meets the exit conditions is to compare the design variable set formed in the current iteration with the design variable set formed in the previous iteration, or to compare the objective function value obtained in the current iteration with the objective function value obtained in the previous iteration, and the interface tension and compression constraints and volume constraints are both satisfied, or the total number of iterations reaches 400.

7. The multi-material structure topology optimization method as described in claim 6, characterized in that, The exit condition is that the absolute value of the difference between all design variables in the design variable set formed in the current iteration and the corresponding design variables in the design variable set formed in the previous iteration is less than 0.001, or the total number of iterations reaches 300.

8. The multi-material structure topology optimization method as described in claim 1, characterized in that, in S6: in, Set to 10 -6 .

Images