Multi-material structure topological optimization method for accurately predicting and controlling stress

Abstract

The invention discloses a multi-material structure topological optimization method capable of accurately predicting and controlling stress, which comprises the following steps of: introducing multi-material prediction density and a multi-material transformation factor in a step of relaxing stress constraint of a material sub-domain to realize accurate prediction and control on the stress constraint of the multi-material, and solving the stress field when solving the stress field. A Heaviside function is adopted to carry out large-force mapping on the filtering density, the multi-material prediction density is obtained, the number of gray units is smaller than that of physical density, the structural gray change degree in the iteration process is remarkably reduced, a linear or nonlinear stiffness matrix calculation mode can be adaptively selected according to the proportion of the gray units in a design domain, and the multi-material prediction density is obtained. The stress calculation precision is improved; meanwhile, different from a traditional regularization method which adopts a maximum stress value in a design domain in an iteration process to calculate the stress relaxation coefficient, the stress relaxation coefficient is calculated by adopting a maximum predicted stress value of each material sub-domain in the iteration process aiming at basic characteristics of multi-material topological optimization, and accurate control on the maximum stress in the design domain is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of engineering structure optimization, and in particular to a multi-material structure topology optimization method for accurately predicting and controlling stress. Background Art

[0002] As a key load-bearing component in an aircraft engine, the performance of the turbine disk 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 disks have become increasingly difficult to meet the growing performance requirements. To this end, researchers have proposed the concept of multi-material turbine disks, which aims to optimize the performance of turbine disks by combining the characteristics of two or more materials, thereby 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 innovative multi-material structures. In recent years, with the rapid development of additive manufacturing technology (i.e. 3D printing), it has become possible to manufacture these complex shapes of multi-material combinations, further promoting the practical application of multi-material turbine disks.

[0004] However, in the process of multi-material topology optimization design using conventional topology optimization methods considering stress constraints, the following problems still exist: on the one hand, the structural grayscale changes significantly during the iteration process, and the single stiffness matrix interpolation scheme does not match the structural grayscale, resulting in low stress calculation accuracy; on the other hand, the coupling of centrifugal loads and multi-material design variables makes the stress sensitivity highly nonlinear, which makes stress control difficult and is not conducive to optimization convergence. 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 multi-material structure topology optimization method for accurately predicting and controlling stress. The method constructs a method for accurately predicting and controlling multi-material stress by introducing multi-material predicted density and multi-material transformation factor, thereby improving the problems of low stress calculation accuracy and difficult stress control existing in conventional topology optimization methods considering stress constraints.

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

[0007] Technical Solution 1: A multi-material structure topology optimization method that accurately predicts and controls stress, 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, and initializing the design variables x of each unit after discretization. e,mThe 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, and m is the serial number of the corresponding material subdomain unit;

[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 v 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 steps:

[0023] S8.1: Use the Heaviside function to map the filtered density in each material subdomain to obtain the predicted density:

[0024] Among them, β H Controls the jitter of the Heaviside function. Starting from the first iteration, it is updated every 20 iterations. The β value of each update is H The value is the last β H The value is increased by 4, but β H The maximum value is 20;

[0025] S8.2: Calculate multi-material transformation factors:

[0026] Among them, N b The number of elements with predicted density greater than 0.9 for all material subdomains, N w The number of elements with predicted density less than 0.1 for all material subdomains, N m is the number of all material subdomain units;

[0027] S8.3: Calculate the predicted stiffness matrix of the design domain elements:

[0028]

[0029] S8.4: Assemble the global stiffness matrix and calculate the structural displacements based on the static equilibrium equations;

[0030] S8.5: Based on the predicted density Compute the predicted von Mises stresses for material subdomain elements

[0031] S8.6: Penalize the predicted von Mises stress of the element in the material subdomain using a linear penalty to obtain the penalized predicted von Mises stress:

[0032] S8.7: Calculate the maximum predicted stress for this iteration:

[0033] S9: Perform sensitivity analysis on the objective function;

[0034] S10: 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;

[0035] S11: Determine whether this iteration satisfies the exit condition. If so, the iteration is exited and the design variable set formed by this iteration is recorded as the optimal design variable set. If not, the next iteration is performed.

[0036] The process of steps S8.1 to S8.7 above allows parallel calculation with the process of steps S2 to S8.

[0037] Technical Solution 2 based on Technical Solution 1: The objective function is structural flexibility, and the optimization objective is expressed as:

[0038] min C=f T u

[0039] stg S,m =σ a,m -σ L,m ≤0m=1,2,...,N m

[0040]

[0041] Ku=f;

[0042] Where C is the flexibility of the structure, g M,m is the total mass constraint of the structure, M f is the total mass percentage of the structure allowed, Ve 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, σ q,m is the maximum stress of material m, σ L,m is the allowable stress of material m. K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

[0043] Technical Solution 3 based on Technical Solution 1: In step S11, the method for judging whether the current iteration meets the condition for exiting the iteration 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

[0044] 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.

[0045] Technical Solution 5 Based on Technical Solution 1: In step S9, the structural flexibility C is proportional to the physical density The derivative of is: The sensitivity of multi-material stress constraints is calculated using the adjoint method.

[0046] Technical solution 6 based on technical solution 1: 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.

[0047] 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:

[0048] The present invention provides a multi-material structure topology optimization method for accurately predicting and controlling stress. The main innovation point of the method is that in the step of relaxing the stress constraint of the material subdomain, the multi-material prediction density and the multi-material transformation factor are introduced to realize the accurate prediction and control of the multi-material stress constraint, wherein the key is to realize the calculation of the maximum predicted stress of this iteration through steps S8.1 to S8.7; when solving the stress field, the Heaviside function is used to perform large-scale mapping on the filter density to obtain the multi-material prediction density, the number of gray units of which is less than the physical density, the degree of grayscale change of the structure during the iteration is significantly reduced, and the linear or nonlinear stiffness matrix calculation method can be adaptively selected according to the proportion of gray units in the design domain to improve the accuracy of stress calculation; at the same time, unlike the traditional regularization method that uses the maximum stress value in the design domain during the iteration process to calculate the stress relaxation coefficient, in step S8, the maximum predicted stress value of each material subdomain during the iteration process is used to calculate the stress relaxation coefficient c according to the basic characteristics of multi-material topology optimization. m , achieving accurate control of the maximum stress in the design domain. Therefore, the multi-material structure topology optimization method for accurately predicting and controlling stress provided by the present invention constructs a method for accurately predicting and controlling multi-material stress by introducing multi-material prediction density and multi-material transformation factor, thereby improving the problems of low stress calculation accuracy and difficult stress control in conventional topology optimization methods considering stress constraints. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] 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.

[0050] Figure 1 A flowchart of a multi-material structure topology optimization method for accurately predicting and controlling stress provided by an embodiment of the present invention;

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

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

[0053] Figure 4 Schematic diagram of a design domain unit and a material subdomain unit in an embodiment of the present invention;

[0054] Figure 5 A comparison of the single-material optimization result and the multi-material optimization result provided in the embodiments of the present invention. DETAILED DESCRIPTION

[0055] 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.

[0056] 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".

[0057] 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.

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

[0059] 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.

[0060] 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:

[0061] Table 1

[0062]

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

[0064] Table 2

[0065]

[0066] Reference Figure 1 , Figure 1 The specific steps of the multi-material structure topology optimization method for accurately predicting and controlling stress 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, and 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.

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

[0068] 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.

[0069] 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 multi-material design variable x e,m Should satisfy 0≤x e,m ≤1, e is the serial number of the design domain unit, and m is the serial number of the corresponding material subdomain unit. Taking structural flexibility as the objective function, considering the stress constraint of each material and the total mass constraint of the structure, the optimization objective is expressed as follows:

[0070] min C=f T u

[0071] stg S,m =σ a,m -σ L,m ≤0m=1,2,...,N m

[0072]

[0073] Ku=f;

[0074] Where, is the flexibility of the structure, g M,m is the total mass constraint of the structure, M f is the total mass percentage of the structure allowed, 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, σ a,m is the maximum stress of material m, σ L,mis the allowable stress of material m. K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

[0075] 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

[0076] Next, multi-material topology optimization iterative design is performed until the optimal set of design variables is obtained; each iterative design performs the following steps:

[0077] 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.

[0078] In this embodiment:

[0079]

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

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

[0082] 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.

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

[0084]

[0085] 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.

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

[0087] 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;

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

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

[0090] σ 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:

[0091]

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

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

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

[0095]

[0096] 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:

[0097] 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.

[0098] 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.

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

[0100] 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 steps:

[0101] S8.1: Use the Heaviside function to map the filtered density in each material subdomain to obtain the predicted density:

[0102] Among them, β H Controls the jitter of the Heaviside function. Starting from the first iteration, it is updated every 20 iterations. The β value of each update is H The value is the last β H The value is increased by 4, but β H The maximum value is 20; the β of the first iteration H is 4;

[0103] S8.2: Calculate multi-material transformation factors:

[0104] Among them, N b The number of elements with predicted density greater than 0.9 for all material subdomains, N w The number of elements with predicted density less than 0.1 for all material subdomains, N m is the number of all material subdomain units;

[0105] S8.3: Calculate the predicted stiffness matrix of the design domain elements:

[0106]

[0107] S8.4: Assemble the global stiffness matrix and calculate the structural displacements based on the static equilibrium equations;

[0108] S8.5: Based on the predicted density Compute the predicted von Mises stresses for material subdomain elements

[0109] S8.6: Penalize the predicted von Mises stress of the element in the material subdomain using a linear penalty to obtain the penalized predicted von Mises stress:

[0110] S8.7: Calculate the maximum predicted stress for this iteration:

[0111] S9: Perform sensitivity analysis on the objective function:

[0112]

[0113] Structural equivalent load f The derivative of is:

[0114]

[0115] Where G is the transformation matrix from the element stiffness matrix to the global stiffness matrix.

[0116] The global stiffness matrix K of the design domain is The derivative of is:

[0117] Stiffness matrix K of design domain element e Physical density Derivative:

[0118] in The expression is as follows:

[0119] Combined with the above formula, the structural flexibility C is proportional to the physical density The derivative of is:

[0120]

[0121] Compute the sensitivity of multi-material stress constraints using the adjoint method:

[0122] Augmented form of the maximum stress approximation in a material subdomain can be written as: where λ i,m is an arbitrary adjoint vector.

[0123] Augmented form of the maximum stress approximation in a material subdomain The physical density of the cells on this subdomain The derivative of is:

[0124]

[0125] Wherein, the maximum approximate stress in the material subdomain is The physical density of the cells on this subdomain The derivative of is:

[0126]

[0127] By choosing a suitable adjoint vector λ i,m , we can get the adjoint equation:

[0128] Wherein, the maximum approximate stress in the material subdomain is The derivative of the displacement vector u in the design domain is:

[0129] Approximate maximum stress in a material subdomain For the displacement vector u e Derivative:

[0130] Adjoint vector λ m,m The calculation expression is as follows:

[0131] Augmented form of the maximum stress approximation in a material subdomain Physical density of cells on other subdomains Derivative:

[0132] Choose a suitable adjoint vector λ j,m , we can get the adjoint equation: λ m =λ m,m =λ j,m ;

[0133] Due to the maximum approximate stress in the material subdomain m Physical density of cells on other subdomains Not relevant, that is:

[0134] In summary, the material subdomain approximates the augmented form of the maximum stress The derivative of the physical density is:

[0135]

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

[0137] S10: 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;

[0138] S11: Determine whether this iteration satisfies the exit condition. If so, the iteration is exited and the design variable set formed by this iteration is recorded as the optimal design variable set. If not, the next iteration is performed.

[0139] The process of steps S8.1 to S8.7 above allows parallel calculation with the process of steps S2 to S8.

[0140] In step S11, the method for judging whether the current iteration satisfies the condition for exiting the iteration 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 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.

[0141] 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.

[0142] In this embodiment, GH4169 and K418B are both currently used turbine disc manufacturing materials, and there is no obvious difference between the two. GH4169 has a higher yield strength than K418B, showing better strength performance, while K418B has a smaller density and a larger elastic modulus, showing better stiffness performance. In this embodiment, the filter radius is set to 5, the safety factor is 1, and the total mass constraint is not more than 13.0%.

[0143] After the above optimization process, the topology optimization results of the multi-material turbine disk structure provided in this embodiment are as follows: Figure 5 shown. Figure 5The optimization results shown include the optimization results based on GH4169 material, the optimization results based on K418B material, and the optimization results based on two materials. It can be seen that: (1) When a single material is used, the optimized structure based on GH4169 can meet the mass constraint and stress constraint, while the optimized structure based on K418B cannot further meet the mass constraint while meeting the stress constraint. This is because the strength performance of K418B is poor and it cannot meet the set stress requirements when the total mass percentage does not exceed 13.0%. (2) The optimized structure based on K418B shows significantly smaller flexibility, which fully demonstrates the excellent stiffness performance of K418B. (3) The optimized structure based on two materials can meet the corresponding stress constraints and mass constraints, and the structural flexibility is reduced by 6.6% compared with the structure based on GH4169, making full use of the characteristics of the two materials. It can be seen that the multi-material structure topology optimization method provided in this embodiment reasonably utilizes the performance of GH4169 and K418B. After meeting the corresponding stress and lightweight design requirements, the structural stiffness is improved to a certain extent, and a structure with better performance is obtained compared with the single material design.

[0144] Based on the above description, the multi-material structure topology optimization method for accurately predicting and controlling stress involved in the embodiment of the present invention has a main innovation point that, in the step of relaxing the stress constraints of the material subdomain, the multi-material prediction density and the multi-material transformation factor are introduced to achieve accurate prediction and control of the multi-material stress constraints, wherein the key is to achieve the calculation of the maximum predicted stress of this iteration through steps S8.1 to S8.7; when solving the stress field, the Heaviside function is used to perform a large-scale mapping on the filter density to obtain the multi-material prediction density, whose number of grayscale units is less than that of the physical density, and the degree of grayscale change of the structure during the iteration process is significantly reduced, and the linear or nonlinear stiffness matrix calculation method can be adaptively selected according to the proportion of grayscale units in the design domain to improve the accuracy of stress calculation; at the same time, unlike the traditional regularization method that uses the maximum stress value in the design domain during the iteration process to calculate the stress relaxation coefficient, in step S8, the maximum predicted stress value of each material subdomain during the iteration process is used to calculate the stress relaxation coefficient c according to the basic characteristics of multi-material topology optimization. m , achieving accurate control of the maximum stress in the design domain. Therefore, the multi-material structure topology optimization method for accurately predicting and controlling stress provided by the present invention constructs a method for accurately predicting and controlling multi-material stress by introducing multi-material prediction density and multi-material transformation factor, thereby improving the problems of low stress calculation accuracy and difficult stress control in conventional topology optimization methods considering stress constraints.

[0145] 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 multi-material structure topology optimization method for accurately predicting and controlling stress, wherein the optimization goal is to minimize the structural mass under stress constraints, and the method is characterized by: 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, and m is the serial number of the corresponding material subdomain unit; 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 steps: S8.1: Use the Heaviside function to map the filtered density in each material subdomain to obtain the predicted density: Among them, β H Controls the jitter of the Heaviside function. Starting from the first iteration, it is updated every 20 iterations. The β value of each update is H The value is the last β H The value is increased by 4, but β H The maximum value is 20; S8.2: Calculate multi-material transformation factors: Among them, N b The number of elements with predicted density greater than 0.9 for all material subdomains, N w The number of elements with predicted density less than 0.1 for all material subdomains, N m is the number of all material subdomain units; S8.3: Calculate the predicted stiffness matrix of the design domain elements: S8.4: Assemble the global stiffness matrix and calculate the structural displacements based on the static equilibrium equations; S8.5: Based on the predicted density Compute the predicted von Mises stresses for material subdomain elements S8.6: Penalize the predicted von Mises stress of the element in the material subdomain using a linear penalty to obtain the penalized predicted von Mises stress: S8.7: Calculate the maximum predicted stress for this iteration: S9: Perform sensitivity analysis on the objective function; S10: 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; S11: Determine whether this iteration satisfies the exit condition. If so, the iteration is exited and the design variable set formed by this iteration is recorded as the optimal design variable set. If not, the next iteration is performed. The process of steps S8.1 to S8.7 above allows parallel calculation with the process of steps S2 to S8.

2. A multi-material structure topology optimization method for accurately predicting and controlling stress as claimed in claim 1, characterized in that: The objective function is the structural flexibility, and the optimization objective is expressed as: min C=f T u stg S,m =s a,m -s L,m ≤0m=1,2,...,N m Ku=f Where C is the flexibility of the structure, g M,m is the total mass constraint of the structure, M f is the total mass percentage of the structure allowed, 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, σ a,m is the maximum stress of material m, σ L,m is the allowable stress of material m. K is the overall stiffness matrix, u is the node displacement vector, and f is the equivalent node load vector.

3. A multi-material structure topology optimization method for accurately predicting and controlling stress as claimed in claim 1, characterized in that: In step S11, the method for judging whether the current iteration satisfies the exit condition 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 4. A multi-material structure topology optimization method for accurately predicting and controlling stress 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. A multi-material structure topology optimization method for accurately predicting and controlling stress as claimed in claim 1, characterized in that: In step S9, the structural flexibility C is proportional to the physical density The derivative of is: The sensitivity of multi-material stress constraints is calculated using the adjoint method.

6. A multi-material structure topology optimization method for accurately predicting and controlling stress 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