A thermoelastic topology optimization method for thermal stress singularities

Abstract

A thermal elastic topology optimization method for thermal structure under thermal stress singular problem belongs to the technical field of thermal elastic topology optimization of thermal structure. First, the finite element modeling of the thermal structure is carried out and topology optimization is carried out, and a clear configuration is realized; then, a thermal elastic topology optimization model is established, and a "super spring unit" is used to connect the nodes on the original fixed boundary to provide weak stiffness constraint; the final thermal deformation response is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; sensitivity analysis is carried out according to the thermal deformation response; finally, the design variables are updated and iterated based on the sensitivity information and the gradient optimization algorithm, and the maximum iteration step and the convergence criterion are set; the units below the threshold value after optimization are deleted, and the units above the threshold value are retained, and the innovative topology structure of the optimized thermal structure is obtained. The present application can avoid the hinge phenomenon, the gray problem and the material effect of the design using the fixed boundary condition, which will be a beneficial supplement to the existing thermal elastic topology optimization method.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of thermal-elastic topology optimization method of thermal structure, and relates to a thermal-elastic topology optimization method for thermal stress singularity problem. BACKGROUND

[0002] Thermal structure widely exists in the service environment of industrial equipment structural components. The optimization design method for thermal structure has certain theoretical significance and engineering application value.

[0003] Thermal load depends on the material layout on the design domain and will change in the optimization process, which leads to typical design-dependent thermal load problems such as gray scale and no material effect. Since the essential influencing factors of design-dependent thermal load problem are complex, the research on thermal-elastic topology optimization method for thermal structure has been one of the more difficult research directions for a long time.

[0004] From the perspective of optimization, in order to solve the design-dependent thermal load problem existing in the thermal-elastic topology optimization of thermal structure, the existing literature and research usually adopt two methods: interpolation model method and optimization objective method. The interpolation model method establishes the matching relationship between the unit thermal load and the unit stiffness by adjusting the penalty coefficient, so as to avoid the occurrence of gray phenomenon due to the mismatch of the relationship. The optimization objective method selects or establishes a suitable objective function to limit the direct participation of thermal load in the optimization objective, so as to solve the given optimization problem in a direct or indirect way. Although the interpolation model method and the optimization objective method overcome the design-dependent thermal load problem existing in the thermal-elastic topology optimization of thermal structure to a certain extent, due to the existence of thermal stress singularity problem, it is still extremely challenging to solve the typical design-dependent thermal load problems such as gray scale and no material effect.

[0005] Thermal stress singularity problem refers to that the use of too simplified fixed boundary condition will ignore the actual existing boundary thermal expansion, so that the constraint boundary is too rigid, and then under the action of thermal expansion of thermal structure, the boundary and the interior of thermal structure have high thermal stress, which is obviously much higher than the actual thermal stress, and even the thermal deformation mode does not conform to the actual mode. Under the temperature field of high temperature and large temperature gradient, the thermal stress singularity problem is extremely prominent, which will produce thermal stress much higher than the actual thermal stress, make the already existing design-dependent thermal load problem worse or cause the design-dependent thermal load problem that would not occur, and finally lead to the formation of unreasonable design of thermal-elastic topology optimization of thermal structure.

[0006] From the perspective of analysis, in order to solve the design-dependent thermal load problem existing in the thermal-elastic topology optimization of thermal structure, the existing literatures and researches usually adopt two kinds of strategies: symmetry strategy and spring adding strategy. The symmetry strategy utilizes the symmetric boundary condition to allow the actual boundary thermal expansion, but is restricted by the non-symmetric structure geometry and non-symmetric load, and is only applicable to the symmetric structure and symmetric load, so the applicable range is limited. The spring adding strategy utilizes the spring with weak stiffness coefficient to replace the fixed boundary of the original thermal structure, so as to allow the boundary thermal expansion, but the rigid body displacement is prone to occur under the high mechanical load, so that the thermal-elastic finite element analysis is difficult to converge.

[0007] In summary, the existing literatures and researches have certain deficiencies in facing how to solve the thermal stress singularity problem, and considering the high-precision modeling and the uncertainty of related parameters, therefore, a thermal-elastic topology optimization method for the thermal stress singularity problem is proposed, which can provide effective ideas and method means for solving the design-dependent thermal load problem existing in the thermal-elastic topology optimization of thermal structure by means of simple and effective boundary thermal expansion representation, and has certain theoretical significance and engineering application value. SUMMARY

[0008] In the thermal-elastic topology optimization of thermal structure, the existing literatures and researches usually adopt the fixed boundary condition, but the too simplified fixed boundary condition is difficult to represent the actual existing boundary thermal expansion, so that the calculated thermal stress is much higher than the actual value, and even the deformation mode is different from the actual mode, finally resulting in unreasonable design. In order to solve this problem, the present application proposes a thermal-elastic topology optimization method for the thermal stress singularity problem, which is based on the sequential coupling-based representation for boundary thermal expansion (SC-BTE) to eliminate the limitation of the too simplified fixed boundary condition. The core idea is to divide the thermal-elastic finite element balance equation into mechanical load equation and thermal load balance equation, wherein the mechanical load equation adopts the fixed boundary condition, but the thermal load equation adopts the free thermal expansion boundary condition; on this basis, the Helmholtz type anisotropic filtering method is utilized to construct the stiffening feature, so as to obtain reasonable stiffening layout configuration.

[0009] To achieve the above object, the technical scheme of the present application is as follows:

[0010] A thermal-elastic topology optimization method for the thermal stress singularity problem, the thermal-elastic topology optimization method firstly performs finite element modeling on the thermal structure, including an initial geometric model and a design domain; then, topology optimization is performed on the basis of the established finite element model, and in the process of topology optimization, an interpolation model is adopted to respectively punish the Young's modulus and the thermal conductivity of the element and thermal stress coefficient β to obtain 0-1 distribution of the physical variable of the unit A clear configuration is realized; then, a thermoelastic topology optimization model is established, and the original thermoelastic equilibrium equation is divided into two equations, i.e., a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts a fixed boundary condition, but the thermal load equation adopts a free thermal expansion boundary condition; in the thermal load equation, a new unit called "super spring unit" is established, and the nodes on the original fixed boundary are connected by the "super spring unit" to provide a weak stiffness constraint and avoid rigid body displacement caused by free thermal expansion, so that the analysis of the thermal load equation can stably converge; the final thermal deformation response is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; and sensitivity analysis is performed according to the thermal deformation response; finally, the design variables are updated and iterated based on the sensitivity information and the gradient optimization algorithm; the maximum iteration step and the convergence criterion are set; by setting a threshold, the units below the threshold after optimization are deleted, and the units above the threshold are retained, so that the innovative topology structure of the thermal structure after optimization is obtained. Comprising the following steps:

[0011] Step 1: finite element modeling of the thermal structure, including the initial geometric model and the design domain.

[0012] The design domain refers to the area for topology optimization, and modeling can obtain the geometric file of the thermal structure.

[0013] Further, in step 1, the geometric file format of the thermal structure established is arbitrary, including but not limited to x_t, igs, step, etc.

[0014] Step 2: topology optimization based on the established finite element model, and in the process of topology optimization, the interpolation model is used to respectively punish the Young's modulus and thermal stress coefficient β to obtain 0-1 distribution of the physical variable of the unit A clear configuration is realized.

[0015] Further, in step 2, the interpolation model established is arbitrary, including but not limited to: SIMP, RAMP, MRAMP, etc.

[0016] Step 2.1: topology optimization based on the finite element model established in step 1, and in the process of topology optimization, specifically, MRAMP interpolation model is used to respectively punish the Young's modulus and thermal stress coefficient β, which is expressed as follows:

[0017]

[0018] where S, P are penalty coefficients, S = 16, P = 2; E0is the Young's modulus of the solid element ; a0represents the thermal expansion coefficient of the solid element ; represents the Young's modulus of the penalized element; E min is the minimum Young's modulus required to avoid numerical singularity in thermoelastic finite element analysis; represents the physical variable of the element; represents the thermal stress coefficient of the penalized element; b0represents the thermal stress coefficient of the solid element .

[0019] Further, in the step 2.1, E min is set to 10 -6 E0.

[0020] Step 2.2: mapping the physical variable of the element with 0-1 distribution provided in step 2.1 First, the element design variable p e needs to be filtered to obtain the element filtered variable Then, the element filtered variable is mapped to obtain the element physical variable

[0021] First, the element design variable p e is filtered by using a Helmholtz-type anisotropic filter to obtain the element filtered variable which is calculated as follows:

[0022]

[0023] where is the element filtered vector assembled by the element filtered variable ; p * is the element vector assembled by the element vector ; the matrix H is constructed by the element matrix H e . The transformation matrix L and L * are created by the element vector L e and respectively. The relevant part of equation can be expressed as follows:

[0024]

[0025] where N is the element shape function; Ω e is the element domain. The transformation matrix L maps the element vector p * to the node vector Lp * , and then the transformation matrix L * maps the vector H-1 Lρ * Mapping to unit filter vector Matrix c can be defined as follows:

[0026]

[0027] Here, V is defined as a 3×3 positive definite tensor used to determine the filtering direction. n v t1 and v t2 It is a spatial basis vector in the Cartesian coordinate system, v n Along the thickness direction of the thermal structure (designated as direction 3), v t1 Along direction 1 and v n Vertical, v t2 Along direction 2 and v n Vertical; r n Indicates v n Filter radius in direction; r t1 Indicates v t1 Filter radius in direction; r t2 Indicates v t2 The filtering radius in the direction. r t1 =r t2 r n Set to r t1 100 times, which is usually sufficient to filter all cells in the thickness direction of the reinforcement domain; if not, the factor can be increased further.

[0028] Then, filter variables for the cells. Perform Heaviside mapping to obtain unit physical variables The calculation is as follows:

[0029]

[0030] In the formula Let represent the cell filtering variable, β represent the slope of the Heaviside projection function, and η represent the threshold of the Heaviside projection function.

[0031] Step 3: Establish a thermoelastic topology optimization model, divide the original thermoelastic equilibrium equation into two equations, i.e. a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts fixed boundary conditions, but the thermal load equation adopts free thermal expansion boundary conditions; in the thermal load equation, a new unit called "super spring unit" is established to connect the nodes on the original fixed boundary to provide weak stiffness constraints and avoid rigid body displacement caused by free thermal expansion, so that the analysis of the thermal load equation can stably converge; the final thermal deformation response is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; and sensitivity analysis is performed according to the thermal deformation response.

[0032] Step 3.1: Based on the minimum global strain energy design objective, a thermoelastic topology optimization model is established, and the BTE (Boundary thermal expansion) condition is considered, i.e. the original thermoelastic equilibrium equation is divided into two equations, i.e. a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts fixed boundary conditions, but the thermal load equation adopts free thermal expansion boundary conditions. Therefore, the mathematical formula of the established thermoelastic topology optimization model is as follows:

[0033]

[0034] In the formula, Φ is the global strain energy, D is the elastic matrix, ε is the total strain, ε th is the thermal strain, Ω e is the unit domain, V(ρ) is the volume constraint function, v e is the volume of the e-th unit, f upper is the allowable volume fraction, is the total volume of the design domain, ρ is the design variable vector, n is the number of units in the design domain, ρ e is the design variable of the e-th unit. freedofs_m represents the degrees of freedom DOFs (Degrees of freedom) of the mechanical equilibrium equation, and freedofs_th represents the degrees of freedom DOFs of the thermal load equilibrium equation. U is the node displacement vector, U m is the mechanical load related displacement vector, U th is the temperature equivalent load related displacement vector, F m is the node force vector related to the mechanical load, F th is the node force vector related to the temperature equivalent load. For other variables, K is the global stiffness matrix of the structure, K spring is the stiffness matrix of the "super spring unit", K th is the assembled stiffness matrix of K and K spring . The assembled stiffness matrix K th can be calculated as follows:

[0035]

[0036] where fdofs is the DOFs at the original fixed boundary nodes. Then, U th can be calculated as follows:

[0037] K th (freedofs_th,freedofs_th)U th = F th (9)

[0038] where freedofs_th is the remaining DOFs after fdofs and the virtual n+1 DOF are removed.

[0039] Further, in step 3.1, the design variables, the optimization objective and the constraint functions are arbitrary and can be adjusted according to the specific research problem.

[0040] Step 3.2: In the thermal load equation, a new unit called "super spring unit" is established. The nodes on the original fixed boundary are connected by "super spring unit" to provide weak stiffness constraint, avoid rigid body displacement caused by free thermal expansion, and make the analysis of the thermal load equation stable convergence.

[0041] For two-dimensional and three-dimensional structures, the stiffness matrix K spring of the "super spring unit" can be expressed as follows:

[0042]

[0043] where k represents the spring stiffness; m = 2j + 1 and m = 3j + 1 are used for two-dimensional and three-dimensional structures, respectively. Here, j represents the number of original fixed boundary nodes.

[0044] Step 3.3: By superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively, the thermal elastic finite element balance equation is solved, and the final thermal deformation response is calculated.

[0045] Through step 3.1, K th and U th are assembled for solving the thermal elastic finite element balance equation in step 3.1:

[0046]

[0047] The temperature equivalent load vector F th also needs to be assembled. In the process of topology optimization, the equivalent temperature load of the thermal structure in optimization is updated, the uniform temperature rise is considered, and the temperature equivalent load vector Fi th As follows:

[0048]

[0049] where B i is the strain-displacement matrix of the ith element, D i is the elastic matrix of the ith element, Ω is the design domain, V i is the volume of the ith element, V is the total volume of the design domain, is the thermal strain of the ith element, which can be calculated as follows:

[0050]

[0051] where α is the thermal expansion coefficient, T i is the nodal temperature of the ith element, T ref is the reference temperature of the ith element, which is expressed as follows in the range of uniform temperature rise where the change of Young's modulus with temperature is not significant:

[0052]

[0053] where D0 is the elastic matrix under solid material, is the Young's modulus of the ith element, and the temperature equivalent load vector F th is calculated as follows:

[0054]

[0055] Further, in step 3.3, the value of the uniform temperature rise is not fixed and can be adjusted according to the specific research problem.

[0056] Solve formula (11) and obtain the thermal deformation displacement vector U, which is calculated as follows:

[0057]

[0058] where U m is the mechanical load related displacement vector, U th is the temperature equivalent load related displacement vector, and the thermal deformation displacement vector U is the superposition of the mechanical load related displacement vector U m and the temperature equivalent load related displacement vector U th

[0059] Step 3.4: Sensitivity analysis according to thermal deformation response.

[0060] Solve the sensitivity of step 3.1 in the topology optimization process The sensitivity of the global strain energy physical variable under the thermal-mechanical coupling field can be further described as follows:​

[0061]

[0062] where, is the unit physical variable, is the thermal stress coefficient of the unit, is the Young's modulus of the unit, T e is the unit temperature value, T0 is the unit reference temperature, and the default value is 0℃, and φ is the isotropic constant vector [1 1 0 0 0] in three-dimensional problems, is the elastic matrix, Ω e is the unit domain.

[0063] Finally, based on the chain rule, the sensitivity of the global strain energy Φ to the design variable ρ e can be expressed as follows:

[0064]

[0065] The relationship between the unit physical variable and the unit filtered variable is as follows:

[0066]

[0067] where β represents the slope of the Heaviside projection function, and η represents the threshold value of the Heaviside projection function. The relationship between the unit filtered variable and the design variable ρ e is as follows:

[0068]

[0069] where the matrix H is constructed by the unit matrix H e , and the transformation matrices L and L * are created by the unit vectors L e and , respectively.

[0070] Step 4: Update and iterate the design variable based on the sensitivity information and gradient optimization algorithm; set the maximum iteration step and convergence criterion; by setting a threshold value, delete the units below the threshold value after optimization, and retain the units above the threshold value, i.e., obtain the innovative topology structure after thermal structure optimization.

[0071] Step 4.1: Based on the sensitivity information of the sensitivity of the global strain energy Φ of the thermal structure to the design variable ρ e obtained in step 3.4, update and iterate the design variable ρ e using the gradient optimization algorithm.

[0072] Further, in the step 4.1, the gradient-based algorithm includes: bidirectional asymptotic structural optimization algorithm, global mobile asymptotic optimization algorithm and similar optimization methods, which can be adjusted according to specific research.

[0073] Step 4.2: set the maximum iteration step and the convergence criterion.

[0074] In the topology optimization process, after each update iteration of the design variable p e , return to step 2, and then repeat the above steps; if it does not converge to the maximum iteration step loop max , take the design variable p e of the last iteration step; if the convergence criterion is met, converge in advance, and then retain the current converged design variable p e ;

[0075] Step 4.3: by setting a threshold, deleting the unit below the threshold after optimization, and retaining the unit above the threshold, the innovative topology structure after thermal structure optimization can be obtained.

[0076] Further, in the step 4.3, the threshold range is 0-1.

[0077] The beneficial effects of the present application are:

[0078] (1) The present application proposes a sequential coupling-based boundary thermal expansion representation method (SC-BTE) to replace the overly simplified fixed boundary condition to overcome the thermal stress singularity problem in thermal structure topology optimization, aiming at the hinge phenomenon, gray scale problem and material effect caused by the use of fixed boundary condition in thermal elastic topology optimization.

[0079] (2) The present application is a meaningful supplement to the existing thermal elastic topology optimization method, which eliminates the restriction of using symmetry or adopting spring, and can expand the research scope of thermal elastic topology optimization work, including multi-material, multi-scale, intelligence, uncertainty analysis and large scale.

[0080] (3) The present application does not change the original formula of sensitivity analysis, so it can be easily migrated or applied to other topology method types. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 The present application provides a flow chart of a thermal elastic topology optimization method for thermal stress singularity problem;

[0082] Figure 2The folding rudder structure outer wing design domain and boundary condition schematic diagram provided by the embodiment of the present application; Figure 2 (a) is the boundary condition schematic diagram under the fixed boundary condition; Figure 2 (b) is the boundary condition schematic diagram of the thermal load part under the BTE condition; Figure 2 (c) is the boundary condition schematic diagram of the mechanical load part under the BTE condition; wherein, the length L1 is 618.3mm, the length L2 is 295.7mm, the width W1 is 33.7mm, the width W2 is 15.5mm, and the surface pressure load P is 0.2Mpa.

[0083] Figure 3 The folding rudder structure outer wing topology optimization iteration curve schematic diagram provided by the embodiment of the present application; Figure 3 (a) is the iteration curve under the fixed boundary condition; Figure 3 (b) is the iteration curve under the BTE condition;

[0084] Figure 4 The folding rudder structure outer wing topology optimization configuration result schematic diagram provided by the embodiment of the present application; Figure 4 (a) is the design under the fixed boundary condition; Figure 4 (b) is the design under the BTE condition. DETAILED DESCRIPTION

[0085] In order to make the technical problems solved by the present application, the technical solutions adopted and the technical effects achieved more detailed, the present application will be further described in detail below in conjunction with the drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the present application, but not to limit the present application. In addition, it should be noted that, in order to facilitate the description, only the related parts of the present application are shown in the drawings, but not all the contents.

[0086] Figure 1 The implementation flowchart of the thermal elastic topology optimization method for the thermal stress singular problem provided by the present application. As shown in Figure 1As shown, the thermal elastic topology optimization method for thermal stress singular problem provided by the embodiment of the present application comprises: 1) performing finite element modeling on a thermal structure, including an initial geometric model and a design domain; 2) performing topology optimization on the basis of the established finite element model, in the process of topology optimization, respectively punishing the Young's modulus and the thermal stress coefficient of the unit by using an interpolation model, so as to obtain 0-1 distribution of the physical variable of the unit, and realizing clear configuration; 3) establishing a thermal elastic topology optimization model, dividing the original thermal elastic balance equation into two equations, i.e., a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts a fixed boundary condition, but the thermal load equation adopts a free thermal expansion boundary condition; in the thermal load equation, a new unit called "super spring unit" is established, the nodes on the original fixed boundary are connected by using the "super spring unit" to provide weak stiffness constraint, avoiding rigid body displacement caused by free thermal expansion, so that the analysis of the thermal load equation can stably converge; the final thermal deformation response is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; and sensitivity analysis is performed according to the thermal deformation response; 4) updating and iterating the design variables based on the sensitivity information and the gradient optimization algorithm; and setting the maximum iteration step and the convergence criterion; by setting a threshold value, the units below the threshold value after optimization are deleted, and the units above the threshold value are retained, so that the innovative topology structure of the thermal structure after optimization can be obtained. The specific steps are as follows:

[0087] Embodiment: Thermal force coupling topology optimization and stiffening layout design for outer wing of folding rudder structure

[0088] This embodiment performs thermal force coupling topology optimization and stiffening layout design for the outer wing of the folding rudder structure.

[0089] The outer wing of the folding rudder structure belongs to a thermal thin-walled structure, which is subjected to the combined action of axial pressure, bending moment and shear force, and various types of working conditions need to be considered in the design process. The outer wing of the folding rudder structure is designed to be reinforced to improve its carrying capacity in a thermal force coupling environment. First, the finite element model of the outer wing of the folding rudder structure is modeled, including establishing the initial finite element model of the outer wing of the folding rudder structure and the design domain. Then, on the basis of the established finite element model of the outer wing of the folding rudder structure, topology optimization is performed, and based on the interpolation model, the Young's modulus and the thermal stress coefficient β of the unit of the outer wing of the folding rudder structure are respectively punished to obtain 0-1 distribution of the physical variable of the unit A clear topological configuration of the folding rudder structure outer wing is achieved. Subsequently, a thermoelastic topology optimization stiffening design model for the folding rudder structure outer wing is established. The original thermoelastic equilibrium equation is divided into two equations: a mechanical load equation and a thermal load equation. The mechanical load equation uses fixed boundary conditions, while the thermal load equation uses free thermal expansion boundary conditions. The stiffness matrix and nodal force vectors related to the thermal load are assembled, and the mechanical and thermal load equilibrium equations are solved to obtain the final thermal deformation response of the folding rudder structure outer wing, which is then used for sensitivity analysis. Finally, the design variables of the folding rudder structure outer wing are updated and iterated based on the sensitivity information and gradient optimization algorithm. A maximum number of iterations and a convergence criterion are set. By setting a threshold, elements below the threshold in the optimized folding rudder structure outer wing are deleted, while elements above the threshold are retained, thus obtaining the topological configuration of the folding rudder structure outer wing. The invention will be further described in detail below with reference to an embodiment of the topology optimization design of the folding rudder structure outer wing:

[0090] Step 1: Model the finite element model of the folding rudder structure outer wing, including establishing the initial finite element geometric model and design domain of the folding rudder structure outer wing.

[0091] The entire finite element model is meshed into C3D8 elements; the design domain refers to the region where topology optimization is performed; modeling yields the x_t file of the folding rudder structure outer wing.

[0092] Step 2: Based on the established finite element model of the folding rudder structure outer wing, topology optimization is performed. During the topology optimization process, an interpolation model is used to penalize the Young's modulus of the folding rudder structure outer wing elements. And the thermal stress coefficient β, to obtain the unit physical variables with a 0-1 distribution. Achieving a clear folding rudder structure and outer wing topology configuration.

[0093] Step 2.1: Based on the finite element model of the folding rudder structure outer wing established in Step 1, topology optimization is performed. Specifically, during the topology optimization process, the MRAMP interpolation model is used to penalize the Young's modulus of the folding rudder structure outer wing elements. The thermal stress coefficient β is expressed as follows:

[0094]

[0095] Where S and P are penalty coefficients, S = 16, P = 2; E0 is a solid unit. Young's modulus; α0 represents the solid element. The coefficient of thermal expansion; E represents the Young's modulus of the penalized unit. min It is the minimum Young's modulus required to avoid numerical singularities in thermoelastic finite element analysis; Represents the physical variables of the unit; represents the thermal stress coefficient of the penalized element; β0represents the thermal stress coefficient of the solid element .

[0096] Further, in step 2.1, E min is set to 10 -6 E0.

[0097] Step 2.2: Generating element physical variables with 0-1 distribution for step 2.1 First, the element design variable ρ e of the outer wing of the folding rudder structure needs to be filtered to obtain the element filtered variable Then, the element filtered variable is mapped to obtain the element physical variable

[0098] First, the element design variable ρ e is filtered using a Helmholtz-type anisotropic filter to obtain the element filtered variable is calculated as follows:

[0099]

[0100] where, is the element filtered vector, which is assembled from the element filtered variable ; ρ * is the element vector, which is assembled from the element vector ; the matrix H is constructed from the element matrix H e . The transformation matrices L and L * are created from the element vectors L e and The relevant part of equation can be expressed as follows:

[0101]

[0102] where N is the element shape function of the outer wing of the folding rudder structure, Ω e is the element domain. The transformation matrix L maps the element vector ρ * to the node vector Lρ * , and then the transformation matrix L * maps the vector H -1 Lρ * to the element filtered vector The matrix c can be defined as follows:

[0103]

[0104] where V is defined as a 3x3 positive definite tensor, which is used to determine the filtering direction. vn , v t1 and v t2 are space base vectors in Cartesian coordinate system, v n along the thickness direction of the outer wing structure of the folding rudder structure (set as direction 3), v t1 perpendicular to v n along direction 1, v t2 perpendicular to v n along direction 2; r n represents the filter radius of v n direction; r t1 represents the filter radius of v t1 direction; r t2 represents the filter radius of v t2 direction. r t1 = r t2 = 9, r n = 900, which is usually sufficient to filter all elements in the thickness direction of the reinforcement domain, and if not, the multiple can be increased.

[0105] Then, the element filtering variable v of the outer wing structure of the folding rudder structure is subjected to Heaviside mapping to obtain the element physical variable v The calculation is as follows:

[0106]

[0107] In the formula, v is the element filtering variable, β represents the slope of the Heaviside projection function, and η represents the threshold value of the Heaviside projection function. The initial value of η is set to 0.5, and the initial value of β is set to 1, which is multiplied by 1 every 40 iterations until β max = 32.

[0108] Step 3: Establish a thermal elastic topology optimization model of the outer wing of the folding rudder structure, divide the original thermal elastic equilibrium equation into two equations, a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts fixed boundary conditions, but the thermal load equation adopts free thermal expansion boundary conditions; in the thermal load equation, a new unit called "super spring unit" is established, which is used to connect the nodes on the original fixed boundary to provide weak stiffness constraint and avoid rigid body displacement caused by free thermal expansion, so that the analysis of the thermal load equation can be stably converged; the final thermal deformation response of the outer wing of the folding rudder structure is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; and the sensitivity analysis is carried out according to the thermal deformation response.

[0109] Step 3.1: Based on the minimum global strain energy design objective, a thermoelastic topology optimization model is established, and the BTE (Boundary thermal expansion) condition is considered, i.e. the original thermoelastic equilibrium equation is divided into two equations, a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts a fixed boundary condition, but the thermal load equation adopts a free thermal expansion boundary condition. Therefore, the mathematical formula of the established thermoelastic topology optimization model is as follows:

[0110]

[0111] In the formula, Φ is the global strain energy of the folding rudder structure outer wing, D is an elastic matrix, ε is a total strain, ε th is a thermal strain, Ω e is an element domain, V(ρ) is a volume constraint function, v e is the volume of the e th element, f upper is the allowable volume fraction, which is set to 0.3, is the total volume of the design domain, ρ is a design variable vector, n is the number of elements in the design domain, ρ e is the design variable of the e th element. freedofs_m represents the degrees of freedom DOFs (Degrees of freedom) of the mechanical equilibrium equation, and freedofs_th represents the degrees of freedom DOFs of the thermal load equilibrium equation. U is a node displacement vector, U m is a mechanical load related displacement vector, U th is a temperature equivalent load related displacement vector, F m is a node force vector related to the mechanical load, F th is a node force vector related to the temperature equivalent load. For other variables, K is the global stiffness matrix of the structure, K spring is the stiffness matrix of the “super spring element”, K th is the assembled stiffness matrix of K and K spring , and the assembled stiffness matrix K th can be calculated as follows:

[0112]

[0113] Where fdofs is the degree of freedom DOFs at the original fixed boundary node. Then, U th can be calculated as follows:

[0114] K th (freedofs_m,freedofs_m)U th = F th (9)

[0115] where freedofs_th is the remaining degrees of freedom after removing fdofs and the virtual n+1 degree of freedom.

[0116] Step 3.2: In the thermal load equation, a new element called "super spring element" is established to connect the nodes on the original fixed boundary to provide weak stiffness constraint to avoid rigid body displacement caused by free thermal expansion, so that the thermal load equation analysis of the outer wing of the folding rudder structure can be stable convergence.

[0117] For two-dimensional and three-dimensional structures, the stiffness matrix K spring can be expressed as follows:

[0118]

[0119] where k represents the spring stiffness, which is set to 100; m = 2j + 1 and m = 3j + 1 are used for two-dimensional and three-dimensional structures respectively. Here, j represents the number of original fixed boundary nodes.

[0120] Step 3.3: By superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively, the thermal-elastic finite element equilibrium equation is solved, so as to calculate the final thermal deformation response of the outer wing of the folding rudder structure.

[0121] Through step 3.1, the K th and U th are assembled to solve the thermal-elastic finite element equilibrium equation in step 3.1:

[0122]

[0123] The temperature equivalent load vector F th also needs to be assembled. In the process of topology optimization, the equivalent temperature load of the thermal structure in optimization is updated, considering uniform temperature rise of 200℃, and the temperature equivalent load vector F i th of the i-th element is constructed as follows:

[0124]

[0125] where B i is the strain-displacement matrix of the i-th element of the outer wing of the folding rudder structure, D i is the elastic matrix of the i-th element, Ω is the design domain, V i is the volume of the i-th element, V is the total volume of the design domain, is the thermal strain of the i-th element, which can be calculated as follows:

[0126]

[0127] Wherein, a is the thermal expansion coefficient of the folding rudder structure outer wing, T i is the node temperature of the i-th unit, T ref is the reference temperature of the i-th unit, in the uniform temperature rise range where the change of Young's modulus with temperature is not obvious, it is expressed as follows:

[0128]

[0129] Wherein, D0 is the elastic matrix under the solid material, is the Young's modulus of the i-th unit, the temperature equivalent load vector F th is calculated as follows:

[0130]

[0131] Solve formula (11) and obtain the thermal deformation displacement vector U of the folding rudder structure outer wing, the solution calculation is as follows:

[0132]

[0133] Wherein, U m is the mechanical load related displacement vector, U th is the temperature equivalent load related displacement vector, and the thermal deformation displacement vector U is the superposition of the mechanical load related displacement vector U m and the temperature equivalent load related displacement vector U th .

[0134] Step 3.4: sensitivity analysis according to thermal deformation response.

[0135] Solve the sensitivity of step 3.1 in the topology optimization process Under the thermal-mechanical coupling field, the sensitivity of the global strain energy physical variable of the folding rudder structure outer wing can be further described as follows:

[0136]

[0137] Wherein, is the unit physical variable of the folding rudder structure outer wing, is the thermal stress coefficient of the unit, is the Young's modulus of the unit, T e is the unit temperature value, T0 is the unit reference temperature, the default value is 0℃, and φ is the isotropic constant vector [1 1 1 0 0 0] in three-dimensional problems, is the elastic matrix, Ω e is the unit domain.

[0138] Finally, based on the chain rule, the sensitivity of the global strain energy Φ to the design variable ρ e can be expressed as follows:

[0139]

[0140] Unit physical variable Relationship between the unit filter variable and the design variable p

[0141]

[0142] wherein β represents the slope of the Heaviside projection function, and η represents the threshold value of the Heaviside projection function. The unit filter variable and the design variable p e Relationship between the unit filter variable

[0143]

[0144] wherein the matrix H is constructed by the unit matrix H e , and the transformation matrix L and L * are respectively created by the unit vector L e and .

[0145] Step 4: Based on the sensitivity information and the gradient optimization algorithm, the design variable of the folding rudder structure outer wing is updated and iterated; and the maximum iteration step and the convergence criterion are set; by setting a threshold value, the units below the threshold value after optimization are deleted, and the units above the threshold value are retained, that is, the innovative topology structure of the folding rudder structure outer wing after optimization can be obtained.

[0146] Step 4.1: Based on the sensitivity information of the folding rudder structure outer wing global strain energy Φ to the design variable p e obtained in step 3.4, the moving asymptote-based gradient optimization algorithm is used to update and iterate the design variable p e of the folding rudder structure outer wing.

[0147] Step 4.2: Set the maximum iteration step and the convergence criterion.

[0148] The convergence criterion is that the change of the design variable is less than 0.001 or the maximum iteration step loop max = 240 times. In the topology optimization process of the folding rudder structure outer wing, after completing the update and iteration of the design variable p e , return to step 2, and then repeat the above steps; if it does not converge to the maximum iteration step loop max , take the design variable p e of the last iteration step; if the convergence criterion is met and converges in advance, the current converged design variable p e is retained.

[0149] Step 4.3: By setting a threshold, the units of the folding rudder structure outer wing below the threshold after optimization are deleted, and the units above the threshold are retained, that is, the innovative topology structure of the folding rudder structure outer wing after optimization is obtained, and the threshold value of the embodiment is 0.5.

[0150] Figure 4 The folding rudder structure outer wing topology optimization configuration result comparison schematic diagram provided by the embodiment of the present application shows that, compared with the fixed boundary method, the method of the present application can obtain a reasonable configuration, and the specific performance is that the maximum deformation is reduced by 36.01%, and the maximum Mises stress is reduced by 83.21%. Moreover, the volume fraction can finally converge to the set volume fraction, which is 0.3. Figure 3

[0151] Finally, it should be noted that: the above embodiments are only used to illustrate the method scheme of the present application, but not to limit it; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: the method scheme recorded in the foregoing embodiments is modified, or some or all of the method features are replaced, without making the corresponding method scheme deviate from the scope of the method scheme of the embodiments of the present application.​

Claims

1. A thermo-elastic topology optimization method for thermal stress singularities, characterized in that, The thermoelastic topology optimization method comprises the following steps: Step 1: finite element modeling of a thermal structure, including an initial geometric model and a design domain; the design domain refers to a region for topology optimization, and finite element modeling can obtain a geometric file of the thermal structure; Step 2: On the basis of the established finite element model, topology optimization is carried out, and in the process of topology optimization, the interpolation model is used to punish the Young's modulus and the thermal stress coefficient of the element respectively, the 0-1 distributed physical variables of the element are obtained , and clear configuration is realized; Step 3: establishing a thermoelastic topology optimization model, dividing the original thermoelastic equilibrium equation into two equations, i.e., a mechanical load equation and a thermal load equation, wherein the mechanical load equation adopts a fixed boundary condition, but the thermal load equation adopts a free thermal expansion boundary condition; in the thermal load equation, a new unit called "super spring unit" is established, and the "super spring unit" is used to connect the nodes on the original fixed boundary to provide a weak stiffness constraint, so as to avoid rigid body displacement caused by free thermal expansion, so that the analysis of the thermal load equation can stably converge; the final thermal deformation response is obtained by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively; and sensitivity analysis is performed according to the thermal deformation response; For two-dimensional and three-dimensional structures, the stiffness matrix corresponding to the "super spring element" is represented as follows: （8） wherein, represents the spring stiffness; and are used for two- and three-dimensional structures, respectively; here, represents the number of original fixed boundary nodes; Step 4: updating and iterating the design variables based on the sensitivity information and the gradient optimization algorithm; and setting a maximum iteration step and a convergence criterion; by setting a threshold value, the units below the threshold value after optimization are deleted, and the units above the threshold value are retained, so as to obtain the optimized topology structure of the thermal structure.

2. The method of claim 1, wherein, In the step 1, the geometric file format of the thermal structure established includes x_t, igs, step or other geometric file formats.

3. The thermal elastic topology optimization method for thermal stress singular problem according to claim 1, characterized in that, In the step 2, the interpolation model established includes SIMP, RAMP, MRAMP or other commonly used interpolation models.

4. The method of claim 1, wherein, The step 2 is specifically as follows: Step 2.1: On the basis of the finite element model established in step 1, topology optimization is carried out, and in the process of topology optimization, the interpolation model is used to respectively punish the Young's modulus and the thermal stress coefficient , which are expressed as follows: （1） （2） wherein , is a penalty coefficient, = 16, = 2; is the Young's modulus of the solid element; denotes the thermal expansion coefficient of the solid element; denotes the Young's modulus of the penalized element; is the lowest Young's modulus required to avoid numerical singularities in thermoelastic finite element analysis; denotes the physical variable of the element; denotes the thermal stress coefficient of the penalized element; denotes the thermal stress coefficient of the solid element, said solid element = 1; Step 2.2: Provide 0-1 distributed unit physical variables for step 2.1 First, the unit design variables are filtered to obtain unit filtered variables Then, the unit filtered variables are mapped to obtain unit physical variables ; First, the cell design variables are filtered using a Helmholtz-type anisotropic filter to obtain cell filtered variables which are calculated as follows: （3） wherein, is a cell filter vector assembled from cell filter variables ; is a cell vector assembled from cell vectors ; is a matrix constructed from cell matrices ; is a transformation matrix and is created from cell vectors and Then, the unit filtered variables are obtained by Heaviside mapping of the unit physical variables , computed as follows: （4） wherein denotes a filter variable for the unit, denotes a slope of the Heaviside projection function, denotes a threshold of the Heaviside projection function.

5. The method of claim 4, wherein, In the step 2: Said step 2.1 in, is set to ; In step 2.2, the relevant part of the equation is represented as follows: （5） where, is the element shape function; is the element domain; transformation matrix maps the element vector to the node vector and the transformation matrix maps the vector to the element filter vector ; matrix is defined as follows: （6） where, is defined as a 3x3 positive definite tensor that determines the filter direction; , and are spatial basis vectors in the Cartesian coordinate system, is along the thickness direction of the thermal structure, set as direction 3, is perpendicular to along direction 1, is perpendicular to along direction 2; represents the filter radius along direction; represents the filter radius along direction; represents the filter radius along = , is set to 100 times of , which is usually enough to filter all cells in the reinforcement domain thickness direction, if not, continue to increase the multiple.

6. The method of claim 4, wherein, The step 3 is specifically as follows: Step 3.1: based on the minimum global strain energy design target, a thermoelastic topology optimization model is established, and the BTE condition is considered, i.e., the original thermoelastic equilibrium equation is divided into two equations, i.e., a mechanical load equation and a thermal load equation; wherein the mechanical load equation adopts a fixed boundary condition, but the thermal load equation adopts a free thermal expansion boundary condition; therefore, the mathematical formula of the established thermoelastic topology optimization model is as follows: （7） wherein, is the global strain energy, is the elastic matrix, is the total strain, is the thermal strain, is the element domain, is the volume constraint function, is the volume of the th element, is the allowable volume fraction, is the total volume of the design domain, is the design variable vector, is the number of elements in the design domain, is the design variable of the th element; denotes the degrees of freedom DOFs of the mechanical equilibrium equation, denotes the degrees of freedom DOFs of the thermal load equilibrium equation; is the nodal displacement vector, is the nodal displacement vector related to the mechanical load, is the nodal displacement vector related to the temperature equivalent load, is the nodal force vector related to the mechanical load, is the nodal force vector related to the temperature equivalent load; for other variables, is the global stiffness matrix of the structure, is the stiffness matrix of the "super spring element", is and the assembled stiffness matrix; Step 3.2: in the thermal load equation, a new unit called "super spring unit" is established, and the "super spring unit" is used to connect the nodes on the original fixed boundary to provide a weak stiffness constraint; Step 3.3: by superimposing the mechanical displacement component and the thermal displacement component obtained by the mechanical load equation and the thermal load equation respectively, the thermoelastic finite element balance equation is solved, and the final thermal deformation response is calculated; The calculations assembled by step 3.1 and for solving the thermoelastic finite element equilibrium equations in step 3.1 : （9） In the process of topology optimization, the equivalent temperature load of the updated thermal structure in optimization is considered, the uniform temperature rise is considered, and the temperature equivalent load vector of the first unit is constructed as follows: ​ （10） where is the strain-displacement matrix of the th element, is the elasticity matrix of the th element, is the design domain, is the volume of the th element, is the total volume of the design domain, is the thermal strain of the th element, calculated as follows: （11） wherein, is the coefficient of thermal expansion, is the node temperature of the unit, is the reference temperature of the unit, in the range of uniform temperature rise where the change of Young's modulus with temperature is not significant, is expressed as follows: （12） wherein, is the elastic matrix under the solid material, is the Young's modulus of the unit, and the temperature equivalent load vector is calculated as follows: （13） Solving equation (9) and obtaining the thermal deformation displacement vector The solution calculation is as follows: （14） wherein, is a mechanical load related displacement vector, is a temperature equivalent load related displacement vector, thermal distortion displacement vector is a mechanical load related displacement vector is a temperature equivalent load related displacement vector is a superposition of the temperature equivalent load related displacement vector Step 3.4: sensitivity analysis is performed according to the thermal deformation response; Solving step 3.1 sensitivity in topology optimization process The sensitivity of the global strain energy physical variable under the thermo-mechanical coupling field is further described as follows: （15） wherein, is the unit physical variable, is the thermal stress coefficient of the unit, is the Young's modulus of the unit, is the unit temperature value, is the unit reference temperature, with a default value of 0 °C, is the vector of isotropic constants in three-dimensional problems , is the elasticity matrix, is the unit domain; Based on the chain rule, the global strain energy The sensitivity of the design variable is expressed as follows: （16） wherein, is a unit physical variable; is a unit filter variable; The unit physical variable The relationship between the unit filter variable is as follows: （17） wherein denotes the slope of the Heaviside projection function, denotes the threshold of the Heaviside projection function; The unit filtration variable The relationship between the design variable and the unit filtration variable is as follows: （18） where the matrix is constructed from the unit matrix and the transformation matrix and are created from the unit vectors and respectively.

7. The method of claim 6, wherein, In step 3.1, the stiffness matrix is assembled As follows: （19） where, are the DOFs at the original fixed boundary nodes; then, are calculated as follows: （20） wherein, is the number of degrees of freedom removed and virtual remaining degrees of freedom after the degrees of freedom are removed.

8. The method of claim 6, wherein, The step 4 is specifically as follows: Step 4.1: Obtain the thermal structural global strain energy based on the result of step 3.4 sensitivity information of the design variables , the design variables are updated iteratively using a gradient-based optimization algorithm. Step 4.2: setting a maximum iteration step and a convergence criterion; In the topology optimization process, after each update iteration of the design variables is completed, step 2 is returned, and the above steps are repeated next; if it does not converge to the maximum iteration step , the design variables of the last iteration step are taken; if it converges early, the current converged design variables are retained; Step 4.3: by setting a threshold value, the units below the threshold value after optimization are deleted, and the units above the threshold value are retained, so as to obtain the optimized topology structure of the thermal structure.

9. The method of claim 8, wherein, In the step 4.1, the gradient algorithm includes a bidirectional approximate structure optimization algorithm, a global moving asymptote optimization algorithm or other optimization methods.

10. The method of claim 8, wherein, In the step 4.3, the threshold value is a fixed value in the range of 0-1.

Images