Bimaterial structure topological optimization method considering interface stress constraint

Abstract

The invention discloses a bi-material structure topological optimization method considering interface stress constraint, which introduces the concept of interface mark density, maps the filtering density of a unit through a Heaviside function to obtain the interface mark density, establishes an interface identification function based on the interface mark density, and identifies the interface mark density according to the interface identification function. Whether the unit contains two materials at the same time or not can be recognized through the interface recognition function, so that a method capable of recognizing interface domains of different materials in real time in the optimization process is constructed, and a foundation is laid for follow-up interface performance control. According to the optimization method, a method for carrying out constraint control on interface performance in an interface domain is further provided based on recognition of the material interface domain, during constraint control, an interface domain unit elastic coefficient matrix interpolation mode is adopted, the property of an interface unit can be reflected more reasonably, interface stress can be calculated based on an interpolation model, and the interface performance is improved. Constraint control over the interface stress is better achieved.

Description

Technical Field

[0001] The invention relates to the technical field of engineering structure optimization, and in particular to a dual-material structure topology optimization method considering interface stress constraints. Background Art

[0002] As a key load-bearing component in an aircraft engine, the performance of the turbine disc directly affects the overall performance of the engine. However, with the advancement of the technology of the new generation of aircraft engines, traditional single-material turbine discs have become increasingly difficult to meet the growing performance requirements. To this end, researchers have proposed the concept of dual-material turbine discs, which aims to optimize the performance of turbine discs by combining the characteristics of two or more materials, reducing weight while ensuring sufficient strength and durability.

[0003] Multi-material topology optimization design is one of the effective means to achieve this goal. This method can cleverly arrange the distribution of materials with different characteristics based on the actual stress conditions under specific conditions, thereby maximizing the advantages of each material and creating an innovative dual-material structure. In recent years, with the rapid development of additive manufacturing technology (i.e. 3D printing), it has become possible to manufacture these complex-shaped dual-material combinations, further promoting the practical application of dual-material turbine disks.

[0004] However, in the process of multi-material topology optimization design, the interface characteristics at the junction of different materials have an important influence on the overall structural performance. However, current optimization methods generally assume that the material interface is in an ideal state, which leads to deviations between the optimization results and the performance of the actual structure. Summary of the invention

[0005] The purpose of the present invention is to overcome the above-mentioned defects or problems existing in the background technology and to provide a dual-material structure topology optimization method taking into account interface stress constraints. This method reduces the deviation in performance between the optimization result and the actual structure by adding consideration of material interface stress constraints during the dual-material structure topology optimization process.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] Technical Solution 1: A dual-material structure topology optimization method considering interface stress constraints, whose optimization goal is to minimize the structural mass under stress constraints. It is characterized by including the steps of discretizing the structure in the design domain, establishing the material subdomain, initializing the design variables x of each unit after discretization, and e,m The steps of forming an initial variable set, performing iterative design until an optimal variable set is obtained, and obtaining an optimal structure according to the optimal variable set, wherein the design variables x of all units in all design variable sets are e,m Should satisfy 0≤x e,m ≤1;

[0008] In the step of establishing the material subdomain, each design domain unit corresponds to N m material subdomain units, each of which corresponds to a design variable x e,m ; Where, e is the serial number of the design domain unit, m is the serial number of the corresponding material subdomain unit, and the value of m is not greater than 2;

[0009] In the step of iterative design until the optimal variable set is obtained, each iterative design performs the following steps:

[0010] S1: Based on the initial design variable set or the design variable set formed by the last iteration, the design variables of each unit are filtered using the density filtering method to obtain the filtering density of the e-th unit in the material subdomain m;

[0011] S2: Use the Heaviside function to map the filter density in each material subdomain to obtain the physical density:

[0012]

[0013] Among them, β k It is used to control the steepness of the Heaviside function. k is the number of this iteration. When k=1, it is updated every 20 iterations starting from the first iteration. Each updated value is 1.19 times the previous value, and the maximum value is 8. η∈(0,1) is used to control the inflection point of the Heaviside function.

[0014] S3: According to the different Poisson's ratios of different materials, the unit stiffness matrix is ​​calculated by unit interpolation:

[0015] Among them, K e is the stiffness matrix of the design domain element, K m is the unit stiffness matrix of material m, E m is the elastic modulus of material m; q c is the material penalty factor, its value is 4; z e,m =1 means that the e-th unit in the design domain contains the m-th material, z e,m =0 means that the design domain does not contain the mth material at the eth unit;

[0016] S4: Assemble the global stiffness matrix and calculate the structural displacement u according to the static equilibrium equation e ;

[0017] S5: Calculate the von Mises stress of each material subdomain element

[0018] S6: RAMP method is used to penalize the von Mises stress of each material subdomain unit, and the penalized von Mises stress of each material subdomain unit of each unit is obtained:

[0019] Among them, E min,m A small amount added to avoid singularity; q s is the stress penalty factor, and its value is -0.95;

[0020] S7: Aggregate the penalized von Mises stresses of each material subdomain to obtain the aggregated stress: e=1,2,...,N e ;

[0021] S8: Relax the stress constraints of the material subdomain:

[0022] Among them, c m is the stress relaxation coefficient of material subdomain m in this iteration, c m The value of is updated every 5 iterations; is the maximum stress value of material subdomain m, It can be expressed as: Where α is the damping factor, and its value is 0.4; in this iteration, the maximum stress of the material subdomain m Obtained by the following formula:

[0023] S9: Establishing stress constraints in the interface domain between different materials, including the following steps:

[0024] S9.1: Obtaining Interface Marking Density

[0025] Among them, β k It is used to control the steepness of the Heaviside function, and η∈(0,1) is used to control the inflection point of the Heaviside function;

[0026] S9.2: Establish interface recognition function: S e =s e,1 s e,2 ;

[0027] Among them, s e,1 Represents the interface marker density of material subdomain 1 unit, s e,2 represents the interface marker density of the material subdomain 2 unit. e =1, it means that the design domain unit contains both material 1 and material 2; when S e =0, it means that the design domain unit does not contain one of the materials;

[0028] S9.3: Get the elastic coefficient matrix of the interface domain elements:

[0029]

[0030] Among them, D max and Represent the elastic coefficient matrix and unit physical density of the two materials with larger elastic modulus; D min and They represent the elastic coefficient matrix and unit physical density of the material with the smaller elastic modulus in the two materials respectively;

[0031] S9.4: Computation of Element Von Mises Stresses in Interface Domains

[0032] S9.5: Von Mises stress in the aggregate interface domain yields aggregate stress

[0033] S9.6: Relax the interface stress constraints:

[0034] Among them, c in is the interface stress relaxation coefficient, is the iterative maximum interface stress value, and its update method is as follows: Where α is the damping factor, and its value is 0.4;

[0035] S10: Perform sensitivity analysis on the objective function, where the objective function is a function of the design domain unit density, and the design domain unit density is:

[0036] S11: using a moving asymptote algorithm to solve the optimization problem, updating the design variables of each unit, and obtaining and storing the design variable set formed by this iteration;

[0037] S12: judging whether the current iteration satisfies the exit condition of the iteration. If the current iteration satisfies the exit condition of the iteration, the iteration is exited, and the design variable set formed by the current iteration is recorded as the optimal design variable set; if the exit condition of the iteration is not satisfied, the next iteration is performed;

[0038] The process of the above steps S9.1 to S9.6 allows parallel calculation with the process of steps S2 to S8.

[0039] Technical Solution 2 based on Technical Solution 1: The objective function is to minimize the mass fraction of the structure, and the optimization target is expressed as:

[0040]

[0041] stg S,m =σ a,m -σ L,m ≤0m=1,2,...,Nm ;

[0042] g S,in =σ a,in -σ L,in ≤0

[0043] Ku=fe=1,2,...,N e

[0044] Among them, M f is the mass fraction of the structure, V e is the volume of the e-th unit in the design domain, ρ 0,max is the maximum density among the candidate materials, g s,m is the stress constraint of material m, g S,in represents the interface stress constraint, σ a,m is the maximum stress of material m, σ L,m is the allowable stress of material m, σ a,in is the maximum stress in the interface domain, σ L,in is the allowable stress in the interface domain; K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

[0045] Technical solution three based on technical solution one: in step S12, the method for judging whether the current iteration meets the conditions for exiting the iteration is to compare the design variable set formed by the current iteration with the design variable set formed by the previous iteration, or to compare the objective function value obtained by the current iteration with the objective function value obtained by the previous iteration, and the stresses of the material subdomain and the interface domain both meet the stress constraint conditions.

[0046] Technical solution four based on technical solution three: the exit iteration condition is that the absolute value of the difference between all design variables in the design variable set formed by this iteration and the corresponding design variables in the design variable set formed by the previous iteration is less than a preset exit threshold, and the exit threshold is 0.05.

[0047] Technical Solution 5 based on Technical Solution 1: In step S9.4, the von Mises stress of the interface domain unit is: Among them, σ e,in =D e,in Bu e ; V is a constant matrix, u e is the structural displacement in step S4.

[0048] Technical solution 6 based on technical solution 1: In step S9.5, the p-norm function is used to aggregate the von Mises stress of the interface domain to obtain the aggregated stress:

[0049] Technical solution 7 based on technical solution 1: In step S10, the sensitivity of the objective function is The sensitivity of stress constraints in the interface domain is calculated using the adjoint method.

[0050] Technical Solution 8 based on Technical Solution 1: In the step of performing iterative design until the optimal set of variables is obtained, the minimum number of iterations is greater than 200 times and the maximum number of iterations is less than 400 times.

[0051] It can be seen from the above description of the present invention that, compared with the prior art, the present invention has the following beneficial effects:

[0052] The present invention provides a dual-material structure topology optimization method considering interface stress constraints. The main innovation is that the concept of interface marker density is introduced, the filter density of the unit is mapped by the Heaviside function, the interface marker density is obtained, and an interface identification function is established based on the interface marker density. The interface identification function can identify whether the unit contains two materials at the same time, thereby constructing a method that can identify different material interface domains in real time during the optimization process, which lays a foundation for the control of subsequent interface performance. In addition, the optimization method also proposes a method for constraining interface performance in the interface domain based on the identification of the material interface domain. When performing constraint control, the interface domain unit elastic coefficient matrix interpolation method used can more reasonably reflect the properties of the interface unit, so that the interface stress can be calculated based on the interpolation model, and the constraint control of the interface stress can be better realized. Therefore, the dual-material structure topology optimization method considering interface stress constraints provided by the present invention reduces the deviation in performance between the optimization result and the actual structure by adding consideration of the material interface stress constraints in the dual-material structure topology optimization process. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0054] Figure 1 A flow chart of a dual-material structure topology optimization method considering interface stress constraints provided by an embodiment of the present invention;

[0055] Figure 2 The turbine disk body structure to be optimized in the embodiment of the present invention;

[0056] Figure 3 Schematic diagram of the axisymmetric model of the expanded turbine disk and the load conditions in the embodiment of the present invention;

[0057] Figure 4 An interface recognition function established in an embodiment of the present invention;

[0058] Figure 5 The maximum interface stresses in the embodiments of the present invention are 600Pa, 700Pa, 800Pa, 850Pa, 900Pa and 1000Pa, respectively, and the optimization results and interface areas of the mass minimization problem. DETAILED DESCRIPTION

[0059] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are preferred embodiments of the present invention and should not be regarded as excluding other embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0060] In the claims, description and drawings of the present invention, if the terms "include", "have" and their variations are used, they are intended to mean "including but not limited to".

[0061] In the claims and the specification, various terms have different names in the relevant documents and their true meanings should be determined based on the understanding of the process by those skilled in the art combined with the content of the full specification.

[0062] The technical solutions in the embodiments will be described clearly and completely below in conjunction with the accompanying drawings.

[0063] The embodiment of the present invention relates to the optimization design of the turbine disk, and the overall optimization method thereof is referred to Figure 1 , the disc structure of the turbine disc refers to Figure 2 . When the aircraft engine is working, the turbine disk is in a state of high-speed rotation. The load on the turbine disk mainly comes from two aspects: on the one hand, the centrifugal load generated by the high-speed rotation of the turbine disk itself; on the other hand, the tension on the turbine disk generated by the high-speed rotation of the turbine blades. The centrifugal load generated by the turbine itself is symmetrical with the central rotation axis; the tension on the turbine disk caused by the high-speed rotation of the blades can be approximately evenly distributed on the disk edge surface, and the tension is symmetrical to the central rotation axis. The only constraint on the turbine disk is the axial displacement constraint. Therefore, the turbine disk is an axisymmetric model, and the original three-dimensional model can be simplified into a two-dimensional axisymmetric model, which can greatly save computing resources and improve the efficiency of optimization design. At the same time, considering that topology optimization is a conceptual design, the existing design domain can be expanded.

[0064] Figure 3 The schematic diagram of the axisymmetric model of the expanded turbine disk and the load conditions. The shaded area is the non-design area, and the blank area is the initial design area. The geometric parameters of the turbine disk are shown in Table 1:

[0065] Table 1

[0066]

[0067] The design material of the dual-material turbine disk is a virtual material, and the specific parameters are shown in Table 2:

[0068] Table 2

[0069]

[0070] Reference Figure 1 , Figure 1 The specific steps of the dual-material structure topology optimization method considering interface stress constraints provided by an embodiment of the present invention are shown, specifically, including the steps of discretizing the structure in the design domain, establishing the material subdomain, initializing the design variables x of each unit after discretization, e,m The steps include forming an initial variable set, performing iterative design until an optimal variable set is obtained, and obtaining an optimal structure based on the optimal variable set.

[0071] The above steps are described in detail below.

[0072] Firstly, the finite element method is used to discretize the structure in the design domain using a four-node rectangular axisymmetric unit, and the structural design domain is discretized into N e units.

[0073] Next, the material subdomain is established according to the idea of ​​unit superposition method. Each design domain unit corresponds to N m material subdomain units, each of which corresponds to a design variable x e,m ; Among them, the dual material design variable x e,m Should satisfy 0≤x e,m ≤1, e is the serial number of the design domain unit, m is the serial number of the corresponding material subdomain unit, and the value of m is not greater than 2. Considering the dual material stress constraint and interface stress constraint, the minimum mass fraction of the structure is taken as the objective function, and the optimization objective can be expressed as:

[0074] Among them, M f is the mass fraction of the structure, V e is the volume of the e-th unit in the design domain, ρ 0,max is the maximum density among the candidate materials, g S,m is the stress constraint of material m, g S,in represents the interface stress constraint, σ a,m is the maximum stress of material m, σ L,m is the allowable stress of material m, σ a,in is the maximum stress in the interface domain, σ L,in is the allowable stress in the interface domain; K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

[0075] Among them, the objective function is a function of the design domain unit density, and the design domain unit density is:

[0076] Next, initialize the design variables x of each discretized unit e,m , forming an initial design variable set. In this embodiment, the design variables x of each unit in the design domain are e,m are all initialized to 0.5, and the design variables of each unit in the non-design domain are

[0077] Next, the dual-material topology optimization iterative design is performed until the optimal design variable set is obtained; each iterative design performs the following steps:

[0078] S1: Based on the initial design variable set or the design variable set formed by the last iteration, the design variables of each unit are filtered using the density filtering method to obtain the filtering density of the e-th unit in the material subdomain m.

[0079] In this embodiment:

[0080]

[0081] where v i,m is the volume of the ith element of material subdomain m and is calculated as follows:

[0082] w(x i,m )=max(0,r-||c i,m -c e,m || 2min );

[0083] where r min is the given filter radius, c i,m is the center position of the ith unit in the material subdomain m, c e,m is the center position of the e-th unit in the material subdomain m.

[0084] S2: Use the Heaviside function to map the filter density in each material subdomain to obtain the physical density:

[0085]

[0086] Among them, β k It is used to control the steepness of the Heaviside function. k is the sequence number of this iteration. When k=1, starting from the first iteration, it is updated every 20 iterations. Each updated value is 1.19 times the previous value, and the maximum value is 8. η∈(0,1) is used to control the inflection point of the Heaviside function. The value of η is 0.5.

[0087] S3: According to the different Poisson's ratios of different materials, the unit stiffness matrix is ​​calculated by unit interpolation:

[0088] Among them, K e is the stiffness matrix of the design domain element, K m is the unit stiffness matrix of material m, E m is the elastic modulus of material m; q c is the material penalty factor, its value is 4; z e,m =1 means that the e-th unit in the design domain contains the m-th material, z e,m =0 means that the design domain does not contain the mth material at the eth unit;

[0089] S4: Assemble the global stiffness matrix and calculate the structural displacement u according to the static equilibrium equation e ;

[0090] S5: Calculate the von Mises stress of each material subdomain element

[0091] σ e,m =D 0,m Bu e ; Among them, σ e,m is the stress of the e-th element in the material subdomain m, B is the geometric function matrix of the element, u e is the nodal displacement vector of the e-th element in the design domain, D 0,m is the unit elastic coefficient matrix of material m, which can be expressed as:

[0092]

[0093] where μ m is the Poisson's ratio of material m.

[0094] Von Mises stress for element m in material subdomain It can be expressed as:

[0095] Where V is a constant matrix, expressed as:

[0096]

[0097] S6: RAMP method is used to penalize the von Mises stress of each material subdomain unit, and the penalized von Mises stress of each material subdomain unit of each unit is obtained:

[0098] Among them, E min,m A small amount added to avoid singularity; q s is the stress penalty factor, and its value is -0.95.

[0099] S7: Aggregate the penalized von Mises stresses of each material subdomain to obtain the aggregated stress: e=1,2,...,N e ; Wherein, P is the aggregation coefficient, in this embodiment, P = 8.

[0100] S8: Relax the stress constraints of the material subdomain:

[0101] Among them, c m is the stress relaxation coefficient of material subdomain m in this iteration, c m The value of is updated every 5 iterations; is the maximum stress value of material subdomain m, It can be expressed as: Where α is the damping factor, and its value is 0.4; in this iteration, the maximum stress of the material subdomain m Obtained by the following formula:

[0102] S9: Establishing stress constraints in the interface domain between different materials, including the following steps:

[0103] S9.1: Obtaining Interface Marking Density

[0104] Among them, β k It is used to control the steepness of the Heaviside function, and η∈(0,1) is used to control the inflection point of the Heaviside function;

[0105] S9.2: If Figure 4 As shown, establish the interface recognition function: S e =s e,1 s e,2 ;

[0106] Among them, s e,1 Represents the interface marker density of material subdomain 1 unit, s e,2 represents the interface marker density of the material subdomain 2 unit. e =1, it means that the design domain unit contains both material 1 and material 2; when S e =0, it means that the design domain unit does not contain one of the materials;

[0107] S9.3: Get the elastic coefficient matrix of the interface domain elements:

[0108]

[0109] Among them, D max and Represent the elastic coefficient matrix and unit physical density of the two materials with larger elastic modulus; D min and They represent the elastic coefficient matrix and unit physical density of the material with the smaller elastic modulus in the two materials respectively;

[0110] S9.4: Computation of Element Von Mises Stresses in Interface Domains

[0111] The von Mises stress calculation expression of the interface domain element can be expressed as: The stress calculation expression of the interface domain unit can be expressed as: e,in =D e,in Bu e ; V is a constant matrix u e From step S4.

[0112] S9.5: Von Mises stress in the aggregate interface domain yields aggregate stress

[0113] The p-norm function is used to aggregate the von Mises stress in the interface domain to obtain the aggregate stress Expressed as

[0114] S9.6: Relax the interface stress constraints:

[0115] Among them, c in is the interface stress relaxation coefficient, is the iterative maximum interface stress value, and its update method is as follows: Where α is the damping factor, and its value is 0.4.

[0116] S10: Perform sensitivity analysis on the objective function:

[0117] Sensitivity of the objective function

[0118] The sensitivity of the stress constraint in the interface domain is calculated using the adjoint method:

[0119] Augmented form of interface domain stress It can be expressed as: where λ in Any adjoint vector.

[0120] Physical density of the unit The derivative of :

[0121]

[0122] Augmented form of interface domain stress Physical density of the element on material subdomain 1 The derivative of is:

[0123]

[0124] The expression is in

[0125]

[0126] Adjoint vector λ in The expression is as follows:

[0127]

[0128] Augmented form of interface domain stress Physical density of the element on material subdomain 2 The derivative of is:

[0129] Interface labeling density e,m Physical density The derivative of :

[0130]

[0131] The approximate maximum stress at the interface is Physical density of the unit The derivative is:

[0132]

[0133] Calculating Physical Density Using the Chain Rule For the design variable x e,m The derivative of :

[0134]

[0135] S11: using a moving asymptote algorithm to solve the optimization problem, updating the design variables of each unit, and obtaining and storing the design variable set formed by this iteration;

[0136] S12: judging whether the current iteration satisfies the exit condition of the iteration. If the current iteration satisfies the exit condition of the iteration, the iteration is exited, and the design variable set formed by the current iteration is recorded as the optimal design variable set; if the exit condition of the iteration is not satisfied, the next iteration is performed;

[0137] The process of the above steps S9.1 to S9.6 allows parallel calculation with the process of steps S2 to S8.

[0138] Among them, in step S12, the method for judging whether this iteration meets the exit iteration condition is to compare the design variable set formed by this iteration with the design variable set formed by the previous iteration, or to compare the objective function value obtained by this iteration with the objective function value obtained by the previous iteration, and the stresses of the material subdomain and the interface domain all meet the stress constraint conditions. The exit iteration condition is that the absolute value of the difference between all the design variables in the design variable set formed by this iteration and the corresponding design variables in the design variable set formed by the previous iteration is less than a preset exit threshold, and the exit threshold is 0.05.

[0139] In this embodiment, in the step of performing iterative design until the optimal variable set is obtained, the minimum number of iterations is greater than 200 times and the maximum number of iterations is less than 400 times.

[0140] After the above optimization process, the topology optimization results of the dual-material turbine disk structure provided in this embodiment are shown in Table 3:

[0141] Table 3

[0142]

[0143] Reference Figure 5 , which shows the optimization results and interface areas of the mass minimization problem when the maximum interface stress is 600Pa, 700Pa, 800Pa, 850Pa, 900Pa, and 1000Pa. Table III gives the mass scores of the optimization results of the mass minimization problem when the interface stress takes different values.

[0144] Through Table 3 and Figure 5 It can be seen that (1) when the maximum interface stress changes, the topological optimization results of the structure change significantly. (2) When the maximum interface stress is 600Pa, the shape of the structure is similar to a small single inner ring cavity structure, and the structure is asymmetric up and down, and the optimization result is unreasonable. (3) When the maximum interface stress value increases from 700Pa to 1000Pa, the structure changes from a single inner ring cavity structure formed by two materials to a double inner ring cavity structure in which each material forms a cavity. (4) When the maximum interface stress value increases from 700Pa to 850Pa, the score of the optimization result gradually decreases; when the maximum interface stress value increases from 900Pa to 1000Pa, the score of the optimization result does not change significantly.

[0145] Based on the above description, the main innovation of the dual-material structure topology optimization method considering interface stress constraints involved in the embodiment of the present invention is that the concept of interface marker density is introduced, the filter density of the unit is mapped by the Heaviside function, the interface marker density is obtained, and the interface identification function is established based on the interface marker density. The interface identification function can identify whether the unit contains two materials at the same time, thereby constructing a method that can identify different material interface domains in real time during the optimization process, thereby laying the foundation for the control of subsequent interface performance. In addition, the optimization method also proposes a method for constraining interface performance in the interface domain based on the identification of the material interface domain. When performing constraint control, the interface domain unit elastic coefficient matrix interpolation method used can more reasonably reflect the properties of the interface unit, so that the interface stress can be calculated based on the interpolation model, and the constraint control of the interface stress can be better realized. Therefore, the dual-material structure topology optimization method considering interface stress constraints provided by the present invention reduces the deviation in performance between the optimization result and the actual structure by adding consideration of the material interface stress constraints in the dual-material structure topology optimization process.

[0146] The description of the above specification and embodiments is used to explain the protection scope of the present invention, but does not constitute a limitation on the protection scope of the present invention. Through the enlightenment of the present invention or the above embodiments, ordinary technicians in this field can obtain modifications, equivalent substitutions or other improvements to the embodiments of the present invention or part of the technical features thereof through logical analysis, reasoning or limited experiments, which should be included in the protection scope of the present invention.

Claims

1. A dual-material structure topology optimization method considering interface stress constraints, whose optimization goal is to minimize the structural mass under stress constraints, and its characteristics are: It includes the steps of discretizing the structure in the design domain, establishing the material subdomain, and initializing the design variables x of each unit after discretization. e,m The steps of forming an initial variable set, performing iterative design until an optimal variable set is obtained, and obtaining an optimal structure according to the optimal variable set, wherein the design variables x of all units in all design variable sets are e,m Should satisfy 0≤x e,m ≤1; In the step of establishing the material subdomain, each design domain unit corresponds to N m material subdomain units, each of which corresponds to a design variable x e,m ; Where, e is the serial number of the design domain unit, m is the serial number of the corresponding material subdomain unit, and the value of m is not greater than 2; In the step of iterative design until the optimal variable set is obtained, each iterative design performs the following steps: S1: Based on the initial design variable set or the design variable set formed by the last iteration, the design variables of each unit are filtered using the density filtering method to obtain the filtering density of the e-th unit in the material subdomain m; S2: Use the Heaviside function to map the filter density in each material subdomain to obtain the physical density: Among them, β k It is used to control the steepness of the Heaviside function. k is the number of this iteration. When k=1, it is updated every 20 iterations starting from the first iteration. Each updated value is 1.19 times the previous value, and the maximum value is 8. η∈(0,1) is used to control the inflection point of the Heaviside function. S3: According to the different Poisson's ratios of different materials, the unit stiffness matrix is ​​calculated by unit interpolation: Among them, K e is the stiffness matrix of the design domain element, K m is the unit stiffness matrix of material m, E m is the elastic modulus of material m; q c is the material penalty factor, its value is 4; z e,m =1 means that the e-th unit in the design domain contains the m-th material, z e,m =0 means that the design domain does not contain the mth material at the eth unit; S4: Assemble the global stiffness matrix and calculate the structural displacement u according to the static equilibrium equation e ; S5: Calculate the von Mises stress of each material subdomain element S6: RAMP method is used to penalize the von Mises stress of each material subdomain unit, and the penalized von Mises stress of each material subdomain unit of each unit is obtained: Among them, E min,m A small amount added to avoid singularity; q s is the stress penalty factor, and its value is -0.95; S7: Aggregate the penalized von Mises stresses of each material subdomain to obtain the aggregated stress: S8: Relax the stress constraints of the material subdomain: Among them, c m is the stress relaxation coefficient of material subdomain m in this iteration, c m The value of is updated every 5 iterations; is the maximum stress value of material subdomain m, It can be expressed as: Where α is the damping factor, and its value is 0.4; in this iteration, the maximum stress of the material subdomain m Obtained by the following formula: S9: Establishing stress constraints in the interface domain between different materials, including the following steps: S9.1: Obtaining Interface Marking Density Among them, β k It is used to control the steepness of the Heaviside function, and η∈(0,1) is used to control the inflection point of the Heaviside function; S9.2: Establish interface recognition function: S e =s e,1 s e,2 ; Among them, s e,1 Represents the interface marker density of material subdomain 1 unit, s e,2 represents the interface marker density of the material subdomain 2 unit. e =1, it means that the design domain unit contains both material 1 and material 2; when S e =0, it means that the design domain unit does not contain one of the materials; S9.3: Get the elastic coefficient matrix of the interface domain elements: Among them, D max and Represent the elastic coefficient matrix and unit physical density of the two materials with larger elastic modulus; D min and They represent the elastic coefficient matrix and unit physical density of the material with the smaller elastic modulus in the two materials respectively; S9.4: Computation of Element Von Mises Stresses in Interface Domains S9.5: Von Mises stress in the aggregate interface domain yields aggregate stress S9.6: Relax the interface stress constraints: Among them, c in is the interface stress relaxation coefficient, is the iterative maximum interface stress value, and its update method is as follows: Where α is the damping factor, and its value is 0.4; S10: Perform sensitivity analysis on the objective function, where the objective function is a function of the design domain unit density, and the design domain unit density is: S11: using a moving asymptote algorithm to solve the optimization problem, updating the design variables of each unit, and obtaining and storing the design variable set formed by this iteration; S12: judging whether the current iteration satisfies the exit condition of the iteration. If the current iteration satisfies the exit condition of the iteration, the iteration is exited, and the design variable set formed by the current iteration is recorded as the optimal design variable set; if the exit condition of the iteration is not satisfied, the next iteration is performed; The process of the above steps S9.1 to S9.6 allows parallel calculation with the process of steps S2 to S8.

2. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 1, characterized in that: The objective function is to minimize the mass fraction of the structure, and the optimization objective is expressed as: stg S,m =s a,m -s L,m ≤0m=1,2,...,N m ; g S,in =s a,in -s L,in ≤0 Ku=fe=1,2,...,N e Among them, M f is the mass fraction of the structure, V e is the volume of the e-th unit in the design domain, ρ 0,max is the maximum density among the candidate materials, g s,m is the stress constraint of material m, g S,in represents the interface stress constraint, σ a,m is the maximum stress of material m, σ L,m is the allowable stress of material m, σ a,in is the maximum stress in the interface domain, σ L,in is the allowable stress in the interface domain; K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

3. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 1, characterized in that: In step S12, the method for judging whether this iteration satisfies the conditions for exiting the iteration is to compare the design variable set formed by this iteration with the design variable set formed by the previous iteration, or to compare the objective function value obtained by this iteration with the objective function value obtained by the previous iteration, and the stresses of the material subdomain and the interface domain both satisfy the stress constraint conditions.

4. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 3, characterized in that: The exit iteration condition is that the absolute value of the difference between all design variables in the design variable set formed by this iteration and the corresponding design variables in the design variable set formed by the previous iteration is less than a preset exit threshold, and the exit threshold is 0.

05.

5. The method for topological optimization of a dual-material structure considering interface stress constraints according to claim 1, characterized in that: In step S9.4, the von Mises stress of the interface domain element is: Among them, σ e,in =D e,in Bu e ; V is a constant matrix, u e is the structural displacement in step S4.

6. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 1, characterized in that: In step S9.5, the von Mises stress of the interface domain is aggregated using the p-norm function to obtain the aggregated stress:

7. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 1, characterized in that: In step S10, the sensitivity of the objective function is The sensitivity of stress constraints in the interface domain is calculated using the adjoint method.

8. A dual-material structure topology optimization method considering interface stress constraints as claimed in claim 1, characterized in that: In the step of performing iterative design until the optimal variable set is obtained, the minimum number of iterations is greater than 200 times and the maximum number of iterations is less than 400 times.

Images