Force-heat coupling simulation and structure optimization method and system for ultra-complex structure

Abstract

The invention discloses a force-heat coupling simulation and structure optimization method and system for an ultra-complex structure, and the method comprises the steps: generating parameters and data structures required by subsequent simulation and optimization for a given model and a design domain; comprising the following steps: constructing a background coarse grid, constructing a seed point and a radius, modeling a porous structure based on a voronoi structure, and pre-calculating a stress field of the porous structure; calculating a density field, performing thermodynamic simulation on the generated porous structure based on the background coarse grid and the shape function according to the density field, querying a stress matrix corresponding to each seed point according to the obtained stress field, and performing characteristic decomposition and extraction on the stress matrix to obtain a weight matrix corresponding to the seed point; a topological optimization method is used, control parameters are updated through a GCMMA method according to the sensitivity, and therefore the porous structure is updated, material distribution is optimized, and optimization is conducted continuously till the change of the control parameters is smaller than a threshold value. According to the method, the topological optimization problem of the anisotropic porous structure under the multi-physical-field condition can be quickly and effectively solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of topology optimization, and relates to a method and system for mechanical-thermal coupling simulation and structural optimization of an ultra-complex structure. Background Art

[0002] Porous foam materials have received extensive attention in various fields due to their unique properties and potential applications. The porous structure of these materials enables them to have a wide range of functions, such as enhanced strength-to-weight ratio, thermal insulation, and sound absorption. As a result, porous foams have become an important component in many industries, including aerospace, automotive, and construction.

[0003] Topology optimization is able to fully exploit the potential of porous structures and thus develop structures with superior performance. By systematically redistributing materials within a given design space, topology optimization allows the identification of optimal configurations that maximize the desired properties while satisfying specified constraints [1,2]. Through this process, porous structures can be optimized to achieve improved mechanical, thermal, acoustic performance, etc.

[0004] However, traditional voxel-based topology optimization methods have both advantages and limitations when applied to porous foam materials. Advantages include the ease of modeling structures and the ability to capture complex geometries. However, these methods have limitations in terms of computational efficiency, especially for large-scale optimization problems [3]. Therefore, alternative methods are needed to overcome these limitations and improve the performance of porous foams through topology optimization.

[0005] Anisotropic structures, where materials exhibit different properties in different directions, are of particular interest in a variety of applications. Multiphysics topology optimization of anisotropic porous structures has the potential to advance the development of materials and structures. This optimization can be achieved through a variety of methods, such as thermodynamic orientation optimization [4], augmented hypothetical strain elements [5], and strength-based topology optimization [6]. The role of anisotropy is also important in the topology optimization of additively manufactured load-bearing structures.

[0006] The topology optimization problem of anisotropic porous structures under multi-physics conditions is currently a hot topic in research. [7] A method for designing the optimal arrangement of porous open-cell foams for acoustic and vibration performance using finite element solutions of the Biot equations is discussed. Harston [8] A topology optimization method is proposed to simultaneously optimize the topology and material properties of anisotropic materials to achieve performance improvement. [9] A new method is proposed to optimize the design of porous structures for additive manufacturing and improve their structural properties.

[0007] However, the realization of topology optimization of anisotropic porous structures under multi-physics conditions remains a major challenge. Its complexity lies in the need to consider various physical factors as well as the anisotropic design of porous structures. The realization of optimal design considering the anisotropic behavior of materials in different directions adds additional complexity to the problem. The interaction between different physical fields further complicates the optimization process, and innovative methods need to be developed to effectively address these challenges.

[0008] [1]J.Zhu,J.Sun,H.Tang,J.Wang,Q.Ao,T.Bao,W.Song,Gradient-structuraloptimization of metal fiber porous materials for sound absorption,PowderTechnology 301(2016)1235–1241.doi:10.1016 / j.powtec.2016.08.006.

[0009] [2] J. Yoon, H. Kim, T. Koh, S. Pyo, Microstructural characteristics of soundabsorbable porous cement-based materials by incorporating natural fibers and aluminum powder, Construction and Building Materials 243 (2020) 118167.

[0010] doi:10.1016 / j.conbuildmat.2020.118167.

[0011] [3] S.Dou, Aprojection approach for topology optimization of porousstructures through implicit local volume control, Structural and Multidisciplinary Optimization 62(2)(2020)835–850.doi:10.1007 / s00158-020-02539-x.

[0012] [4]D.R.Jantos,K.Hackl,P.Junker,Topology optimization with anisotropicmaterials,including a filter to smooth fiber pathways,Structural andMultidisciplinary Optimization 61(5)(2020)2135–2154.doi:10.1007 /

[0013] s00158-019-02461-x.

[0014] [5]G.Zhang,L.Li,K.Khandelwal,Topology optimization of structures withanisotropic plastic materials using enhanced assumed strain elements,Structural and Multidisciplinary Optimization 55(6)(2016)1965–1988.

[0015] doi:10.1007 / s00158-016-1612-1.

[0016] [6]A.M.Mirzendehdel,B.Rankouhi,K.Suresh,Strength-based topologyoptimization for anisotropic parts,Additive Manufacturing 19(2018)104–113.

[0017] doi:10.1016 / j.addma.2017.11.007.

[0018] [7]E.Lind-Nordgren,P. Optimising open porous foam foracoustical and vibrational performance,Journal of Sound and Vibration 329(7)(2010)753–767.

[0019] doi:10.1016 / j.jsv.2009.10.009.

[0020] [8] S.Harston, C.Mattson, M.Koecher, Atopology optimization method with anisotropic materials, in: 13th AIAA / ISSMO Multidisciplinary Analysis Optimization Conference, American Institute of Aeronautics and Astronautics,

[0021] 2010.doi:10.2514 / 6.2010-9176.

[0022] [9] Summary of the invention

[0023] In view of the shortcomings of the prior art, the present invention provides a method and system for mechanical and thermal coupling simulation and structural optimization of ultra-complex structures, aiming to solve the topological optimization problem of anisotropic porous structures under multi-physical field conditions. Specifically, the contributions of the present invention can be summarized as follows:

[0024] 1) A porous structure modeling method with clear design-meaningful control parameters is provided to accurately and robustly simulate porous foams.

[0025] 2) A simulation method based on a fixed background grid and CBN shape function is provided to achieve efficient thermodynamic simulation of the model.

[0026] 3) A comprehensive topology optimization framework based on control parameters is provided, which is suitable for the design optimization of multi-objective functions of porous foams.

[0027] 4) A method for optimizing anisotropic porous foams was implemented. The anisotropic optimization of the porous structure was guided by a fixed density field, and the mechanical properties and regularity of the final porous foam were optimized.

[0028] Based on the above contributions 1) to 4), the present invention provides a set of design solutions that can be used for topological optimization of anisotropic porous structures under multi-physical field conditions.

[0029] The technical solution adopted by the present invention is as follows:

[0030] A method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures, including:

[0031] 1) For the given model and design domain, generate and establish the parameters and data structures required for subsequent simulation and optimization; including constructing the background coarse grid, constructing the seed points and radius, modeling the porous structure based on the voronoi structure, and pre-calculating its stress field;

[0032] 2) Calculate the density field, and perform thermodynamic simulation on the generated porous structure based on the background coarse grid and shape function according to the obtained density field; in the process of thermodynamic simulation, calculate the parameter values ​​and sensitivity of the relevant objective function and constraint conditions;

[0033] According to the obtained stress field, the stress matrix corresponding to each seed point is queried, and the stress matrix is ​​subjected to characteristic decomposition and extraction processing to obtain the weight matrix corresponding to the seed point. The objective function value and sensitivity characterizing the structural anisotropy are calculated for the porous structure.

[0034] Using topology optimization, the control parameters are updated through the GCMMA method according to the sensitivity, thereby updating the porous structure and optimizing the material distribution.

[0035] 3) Repeat step 2) until the control parameter change is less than the threshold. When the above process is completed, the final seed point and radius, as well as the topological information of the porous structure are obtained. The above structure is used to construct the final optimized porous foam based on the signed distance field.

[0036] A mechanical-thermal coupling simulation and structural optimization system for ultra-complex structures, used to implement the method described above.

[0037] The beneficial effects of the present invention are:

[0038] 1) The parameterized model based on parameter control generated by the porous structure modeling method provided by the present invention can accurately simulate porous foam. Compared with the voxel model, the simulation and optimization efficiency is improved under the premise of ensuring the characteristic parameters.

[0039] 2) The simulation method based on fixed background grid and CBN shape function used in the present invention is more efficient than the traditional simulation method based on grid.

[0040] 3) The topology optimization framework proposed in this invention avoids expensive global sensitivity calculations and significantly improves computational efficiency.

[0041] 4) The anisotropic porous foam optimization method implemented by the present invention has better mechanical properties compared with the isotropic structure optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 It is a schematic diagram of the overall framework of the method of the present invention;

[0043] Figure 2 The definition of beams in implicit modeling of porous structures includes multiple beams and beam radius description and single beam description;

[0044] Figure 3 A schematic diagram of a design domain and a boundary in a specific example of the present invention;

[0045] Figure 4 A porous foam generated according to the final optimization result in a specific example of the present invention;

[0046] Figure 5 The stress field distribution obtained by simulating the porous structure finally optimized and generated in a specific example of the present invention;

[0047] Figure 6 The temperature field distribution obtained by simulating the porous structure finally optimized and generated in a specific example of the present invention;

[0048] Figure 7 is a convergence curve of the objective function and the constraints in a specific example of the present invention;

[0049] Figure 8 A comparison of the flexibility (stiffness) under anisotropic and isotropic conditions in a specific example of the present invention;

[0050] Fig. 9 1 is a comparison of the average temperature (heat dissipation) under anisotropic and isotropic conditions in a specific example of the present invention. Specific implementation plan

[0051] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments.

[0052] According to a specific embodiment of the present invention, Figure 1 As shown in the figure, a method for mechanical-thermal coupling simulation and structural optimization of an ultra-complex structure includes the following:

[0053] First, for the given model and design domain, the parameters and data structures required for subsequent simulation and optimization are generated and established, such as constructing a background coarse grid through tetrahedral meshing, generating and initializing seed points for building porous structures, etc. At the same time, the model is pre-calculated using a voxel-based topology optimization method to obtain a stress field that guides anisotropic optimization.

[0054] The voronoi structure is generated in the design domain by constructing the background coarse grid, seed points and radius, thereby modeling a porous structure composed of multiple beams with thickness. Furthermore, the distance from each grid unit node to the nearest beam is calculated, and according to the density function, it is converted into the density of each node, and the average density of all nodes of the grid unit is calculated, thereby calculating the density of all grid units and obtaining the density field.

[0055] According to the obtained density field, the generated porous structure is subjected to thermodynamic simulation in sequence using the simulation method based on the background grid and CBN shape function. According to the temperature field obtained from the thermal simulation, the force load on the porous structure under thermal conditions is calculated. The independent force loads applied under the combined force conditions act together on the porous structure, and then a mechanical simulation is performed. During the thermodynamic simulation process, the parameter values ​​and sensitivities of the relevant objective functions and constraints are calculated.

[0056] According to the obtained stress field, the stress matrix corresponding to each seed point is queried, and the stress matrix is ​​subjected to characteristic decomposition and extraction processing to obtain the weight matrix corresponding to the seed point. The objective function value and sensitivity characterizing the structural anisotropy are calculated for the porous structure.

[0057] Furthermore, the topology optimization method is used to update the control parameters through the GCMMA method according to the sensitivity, thereby updating the porous structure and optimizing the material distribution.

[0058] Repeat the above steps until the control parameters change very little or remain unchanged.

[0059] When the above process is completed, the final seed point and radius, as well as the topological information of the porous structure, will be obtained. The above structure is used to construct the final optimized porous foam based on the signed distance field.

[0060] In the above scheme, in the modeling part, the present invention abstracts the model into a beam structure network controlled by the position and radius of the seed point. In order to achieve clear shape control, the present invention uses the edge of the Voronoi diagram with a specific radius to form Voronoi foam.

[0061] In the simulation part, in order to avoid complex grid division, the present invention introduces a fixed background grid. The present invention explores efficient multi-physics field simulation under force-heat coupling conditions and temperature constraints.

[0062] In the optimization part, the present invention studies the generation of anisotropic porous structures in order to maintain a certain regularity while preserving the physical properties, thereby improving the manufacturability of the model.

[0063] On this basis, the present invention establishes a topology optimization framework, which can realize the anisotropic mechanical and thermal coupling topology optimization of porous structures based on Voronoi generation and implicit representation. It is solved according to the classical topology optimization framework and is as follows: Figure 1 The method shown in the figure is carried out. The method first uses a coarse background grid to embed the domain for downstream simulation, and then pre-calculates the stress field through the voxel topology optimization method. After initializing the design variables of seed X and radius r, the variable values ​​are repeatedly updated according to the subsequent steps until convergence: density function calculation, property simulation, sensitivity calculation and design update. Finally, the optimized Voronoi porous structure is constructed by X, r, Ω and output.

[0064] In density function calculation, the present invention constructs an implicit function to represent the Voronoi foam structure and produces smooth foam under explicit geometric control variables. The implicit representation avoids unreliable and very time-consuming geometric operations when explicitly constructing Voronoi foam and generating simulation grids, and allows high-precision continuous finite element integration.

[0065] In the property simulation, the present invention introduces material-aware CBN shape functions to reduce computational cost. The initially constructed embedded coarse mesh is simulated, and a high simulation accuracy is still achieved by tailoring a set of polyhedral material-aware shape functions for each coarse element. In addition, the process operates directly on the implicit foam representation, eliminating the need for error-prone boundary conformal fine mesh generation. This approach significantly reduces computational cost.

[0066] In addition, in order to achieve anisotropic Voronoi foam optimization, the present invention calculates the weight matrix Mx in the objective function based on the pre-calculated stress field to achieve field-guided anisotropic optimization. Compared with the traditional isotropic regular optimization (when Mx is the unit matrix), the method of the present invention can obtain anisotropic Voronoi cells, thereby improving mechanical properties while maintaining regularity.

[0067] The following is a detailed description of the technical details used in the solution.

[0068] 1. Implicit Modeling of Porous Foam

[0069] First, an open Voronoi foam is defined as the edges of the Voronoi cells (with a certain radius) bounded within the model interior. Mixing together the implicit functions of all beams introduced by the Voronoi edges, or together with the shell, an implicit representation of the Voronoi foam defined as a density field is obtained. Details are given below.

[0070] Each Voronoi bubble Ω V(X, r) consists of two parts: the internal part and the boundary part. For the internal part, the present invention first starts from the seed X = {X i ,i=1,…,N s} Compute the Voronoi diagram. In a sufficiently large bounding box, all its edges are then clipped using the outer free shape Ω. The remaining j-th edge becomes a beam with a radius equal to the average of the neighboring seed points:

[0071]

[0072] where χ j is the set of adjacent seeds of edge j, |χ j | is the number of seeds. The radius of the seeds is collected in the form of vector r, such as Figure 2 As shown in the left figure.

[0073] The Voronoi edges corresponding to vertices v1 and v2 define a beam φ, which consists of a beam with a radius of A cylinder with height ||v1-v2|| and two radii The hemispherical end composition: see Figure 2 The implicit form of the beam φ is shown on the right. for:

[0074]

[0075] Where d(x,v1,v2) represents the minimum distance from point x to edges v1,v2.

[0076]

[0077] in

[0078] a=v2-v1, b=x-v1, e=x-v2,

[0079] The implicit representation of the entire open-cell Voronoi foam is obtained by combining the implicit functions of all beams after smoothing using the Kreisselmeier-Steinhauser (KS) function.

[0080]

[0081] Where φ max =max(φ1,…,φ n ), n is the number of beams, and p = 16. The KS function makes the description function Φ(x) compact and differentiable with respect to x.

[0082] The boundary part associated with the external shape Ω is also implicitly represented by φ Ω (x) indicates that when the shell needs to be closed, φ Ω(x) as additional beams, it is straightforward to expand (5). Otherwise, the present invention includes the intersection between the outer shape and the Voronoi faces, where the thickness of each face is the average of the two seed radii that determine the face. For simplicity, the present invention will not distinguish between Voronoi foams with or without shells.

[0083] Voronoi Foam Omega V (X,r) is finally described as a density field by the regularized Heaviside function H(Φ(x)).

[0084]

[0085] where ∈ controls the size of the regularization, α = 1e -6 The default is a small positive number to avoid global stiffness matrix singularities.

[0086] 2. Finite element formula

[0087] In order to realize the simulation of implicit foam, the present invention immerses the model into a fixed tetrahedral background grid D = {D e ,e=1,2,…,N} in the covering shape Ω, that is, This not only avoids the time-consuming and error-prone model partitioning process, but also fixes the mesh to help achieve convergence in subsequent optimizations.

[0088] In each tetrahedral element D e Medium, density H e Set it as the average value of H(x) on its four nodes,

[0089]

[0090] In the formula, For unit D e The coordinates of the kth node in .

[0091] The linear elastic problem of coupled heat conduction problem and thermal expansion effect is considered. The present invention considers a fixed displacement boundary Γ u and a fixed temperature boundary Γ T The design domain Ω. Figure 3 As shown, the surface force acts on the Neumann boundary Γ N The surface heat flux is located at the boundary Γ h The basic equation of the mechanical and thermal coupling problem is

[0092]

[0093] Where u is the displacement vector, σ is the stress vector, T is the temperature, and T0 is Γ T The temperature specified above, ΔT = T – T ref , T refis the thermal expansion reference temperature, k is the thermal conductivity, f is the load vector, and Q is the internal strength of the heat source.

[0094] The stress vector is expressed as:

[0095] σ=D(ε-ε th ) (9)

[0096] Where D is the elastic matrix, ε is the strain vector, and ε th is the thermal strain vector.

[0097] The strain energy E can be expressed as:

[0098]

[0099] After discretizing the design domain by the finite element method, the displacement vector u and the stress tensor σ are expressed as the relationship between the node displacement U:

[0100]

[0101] Where B is the strain displacement matrix and N is the shape function matrix. E is rewritten as:

[0102]

[0103] The total potential energy can be expressed as

[0104]

[0105] Where f v is the object force vector, f b For boundary traction.

[0106] According to the minimum potential energy principle, the change in the node displacement vector U can be expressed as

[0107]

[0108] Considering the variational principle, when f v = 0, the global equations can be obtained from equation 15:

[0109] K m U=F=F m +F th (16)

[0110] have

[0111]

[0112] Where K m is the stiffness matrix, U is the node displacement vector, F m is the mechanical load vector, F th is the heat load vector.

[0113] The static equilibrium equation of the linear thermoelastic continuum under mechanical and temperature loads can be expressed as Equation 16 and Equation 17. m From the element stiffness matrix K m,e Combination

[0114] K m =∑ e K m,e (18)

[0115] K m,e can be calculated as an integral

[0116]

[0117] Where B e is the unit strain-displacement matrix, Ω e is the unit physical domain. Similarly, the thermal load F th By unit heat load F th.e Combination

[0118] F th =∑ e F th,e (20)

[0119] in

[0120]

[0121] Where ε th , e is the unit thermal strain vector, α is the thermal expansion coefficient, T e is the unit temperature, is the unit node temperature T e The average value of Φ is the vector [1 1 1 0 0 0] in three-dimensional space. Similarly, by interpolating the node displacement U, the displacement u(x) in the background grid D can be obtained.

[0122] u(x)=N e (x)U e (twenty two)

[0123] Where N e (x) is the shape function matrix on the unit node, U e is the unit node displacement.

[0124] In the case of force-heat coupling, the effect of force on heat is small. In order to simplify the analysis, the present invention focuses on the one-way coupling of heat and force. Specifically, the present invention incorporates the thermal stress generated in the material into the mechanical equation.

[0125] The heat balance equation for steady-state heat conduction is

[0126] K T (X,r)T=P (23)

[0127] Where K T K is the thermal conductivity matrix, T is the node temperature vector, and P is the heat flux density vector. T and P is given by the unit thermal conductivity matrix \boldK T,e and the unit heat flow vector P e Combination

[0128]

[0129] The thermal conductivity matrix can be derived using the Galerkin finite element method,

[0130]

[0131] Where λ e is the effective thermal conductivity.

[0132] By interpolating the node temperature T, we can get the temperature T(x) in the background grid D.

[0133] T(x)=N e (x)T e (26)

[0134] Where T e is the unit node temperature, satisfying

[0135]

[0136] 3. Optimize question design

[0137] The Voronoi foam design problem uses the seed and radius X,r as the design variables for performance optimization. Its analytical and discrete expressions are explained as follows.

[0138] The present invention hopes that Voronoi foam Ω V The overall deformation of (X, r) is small while maintaining a certain shape regularity. Therefore, the present invention introduces the terms flexibility and shape. The flexibility C(X, r, U) of the elastic potential of the object widely used in topology optimization is set as the physical target.

[0139] C(X,r,U)=F T U (28)

[0140] The shape adjustment energy S(X) is to regularize the Voronoi cells into approximately regular polyhedrons.

[0141]

[0142] That is, each Voronoi cell Seeds X i With the centroid The Euclidean distance between The weighted sum.

[0143] In order to obtain the anisotropic structure of Voronoi foam, the present invention introduces M x As a weight matrix for any coordinate x, different weights are assigned to different directions and Voronoi cells. x The calculation process is as follows. First, simulation and optimization are performed based on the high-resolution voxel model to obtain the stress field of the optimized model. The Cauchy stress tensor at point x in the stress field is given by

[0144]

[0145] where σ ii is the normal stress, σ ij , i≠j is the shear stress. Due to the symmetry of shear stress, σ ij =σ ji , so σ x is a symmetric matrix. In order to extract σ x The main direction vector of , the present invention decomposes its eigenvalue:

[0146] σ x =VΛV -1 (31)

[0147] Where Λ is a diagonal matrix whose eigenvalues ​​are on the diagonal, and V is a matrix whose eigenvectors are its columns. The present invention extracts the eigenvector corresponding to the maximum eigenvalue as the main direction vector v m .

[0148] v m =[v1 v2 v3] (32)

[0149] M x By v m Components:

[0150]

[0151] For simplicity, for each Voronoi cell Take M x The x in the As shown in formula 29.

[0152] Finally, for the optimization problem with volume and temperature constraints, the optimal design parameter p = {p j ,j=1,2,…,n}, where p j Represents a component of a variable X or r:

[0153]

[0154] The objective function J(p,U) in this example is J(X,r,U) which is the weighted sum of the physical objective and the shape regularization term

[0155] J(p, U) = (1-w) C(p, u) + w S(X), (35)

[0156] V and is the total volume of the structure and the maximum volume constraint, T i and are the node temperatures and the corresponding maximum temperature constraints, p j and is the design variable p j The lower and upper bounds of .

[0157] As above, different physical properties or shape regularization metrics can be introduced for different design purposes, and the target foam can be derived similarly following this process.

[0158] The optimization problem Eq. (34) is solved by using the global convergent moving asymptote method (GCMMA) based on numerical gradient due to its robust convergence in design optimization. The original non-convex problem is approximated by a set of convex subproblems using the gradients of the optimization objective and constraints with respect to the design variable p.

[0159] 4. Sensitivity analysis

[0160] 1.1 Sensitivity of flexibility

[0161] To optimize the objective function, first consider the effect of flexibility on the design parameter p j Sensitivity:

[0162]

[0163] In order to eliminate the The present invention obtains the partial derivative of Eq. (23):

[0164]

[0165] Transposing these terms, we get:

[0166]

[0167] Similarly, find the partial derivative of Eq. (16) and rearrange the terms:

[0168]

[0169] For simplicity, assume that P and pj No matter, Substitute Eq.(38) and Eq.(39) into Eq.(36)

[0170]

[0171] In order to simplify Eq. (40), t Defined as neighboring vectors that satisfy:

[0172]

[0173] λ t It can be obtained by solving the following linear equations:

[0174]

[0175] According to Eq.(16) and Eq.(21),

[0176]

[0177] λ t Substituting into Eq. (40), we finally get the flexibility C versus design parameter p j Sensitivity:

[0178]

[0179] According to Eq. (25), Eq. (19) and Eq. (21), the present invention has

[0180]

[0181]

[0182] and

[0183]

[0184] in

[0185]

[0186] 1.2 Temperature Constraint Sensitivity

[0187] To satisfy the temperature constraints, the present invention will then calculate the sensitivity of temperature to the design parameters. The node temperature is written as

[0188]

[0189] in represents a vector of all zeros except the i-th term which is set to 1. The nodal temperature is related to the design parameter p j The gradient of

[0190]

[0191] Derivative of Eq.(23)

[0192]

[0193] Substituting Eq. (51) into Eq. (50), we obtain

[0194]

[0195] To simplify Eq. (52), similar to λ in Eq. (41) t , the present invention introduces an adjacent vector λ i , by solving the following linear equations we can get:

[0196] K T λ i =-L i (53)

[0197] Then, calculate the sensitivity of the node temperature

[0198]

[0199] 5. Efficient calculation of CBN shape function

[0200] Traditional finite element methods require very fine mesh elements to achieve high-precision simulation, which introduces a large number of degrees of freedom. This poses a challenge to solving linear equations.

[0201] Based on previous research, the present invention introduces a material-aware CBN shape function to replace the original shape function, so as to reduce the degree of freedom and thus speed up the calculation.

[0202] First, the entire virtual domain is coarsely meshed. For each coarse unit e and its virtual domain Ω e , the present invention further divides it into a tetrahedral grid to obtain a refined grid h of e. In order to reduce the calculation degree of freedom while maintaining the calculation accuracy, the present invention uses to calculate simulation values ​​such as the stiffness matrix. Subsequently, the present invention compresses the calculation results by mapping the calculation results to the coarse grid.

[0203] Substituting the unit shape function N in Eq. (26) and Eq. (22) e Through the transformation matrix is replaced by the following element shape function consisting of linear shape functions on the fine mesh:

[0204]

[0205] Where N e,h(x) is Ω e The set of node shape functions on the fine grid. Similarly, the strain displacement matrix B in Eq. (25), Eq. (19), and Eq. (21) e Can be replaced by

[0206]

[0207] Where B e,h Ω e The global strain displacement matrix of the fine mesh. CBN transformation matrix Aims to map CBN node values ​​to Ω e The internal value in. By the boundary interpolation matrix Ψ and the boundary-interior transformation matrix The product of is as follows:

[0208]

[0209] where Ψ and Map the displacement of the CBN to the boundary nodes and then to Ω e The fully refined nodes in .

[0210] The boundary interpolation matrix Ψ maps the CBN node values ​​to the fine mesh boundary node values ​​Ω e It is derived by constructing a bicubic B′ezier interpolation surface on the surface of interest, with the CBN as control points. The matrix Ψ is formed by computing the surface at the nodes of the fine mesh in the plane and collecting them row by row into a matrix for all values ​​of the boundary of the mesh h.

[0211] The transformation matrix maps the boundary node values ​​to Ω e The node values ​​of the fine grid in Ω are calculated by the equilibrium equations. e The fine mesh is used for local simulation.

[0212]

[0213] where k b ,k i ,k bi ,k ib Ω e Upper local stiffness matrix k e The sub-matrix, q b ,q i are the vectors of the boundary and internal nodes respectively, and f b Exposed force vectors on boundary nodes formed for harmonic analysis.

[0214] From the second line, the present invention obtains q i =M e q b The relationship between Therefore, q i and q b The set is q=[q b ,q i ] T The present invention has Form

[0215]

[0216] Among them I 2b is the 2b×2b identity matrix.

[0217] Therefore, the original coarse unit thermal conductivity matrix K T,e Given by Eq. (25), Replace B e ,Right now

[0218]

[0219] Among them, K T,e,h is the fine grid Ω e High-fidelity thermal conductivity matrix, H e,α is the average value of H(x) on the tetrahedral elements in the fine grid h, Ω e,α is the virtual domain of the subelement α in the unit e. Similarly, in the coarsening simulation, the gradient of the unit thermal conductivity matrix in Eq. (45) is calculated as

[0220]

[0221] Similarly, K m,e ,F th,e and its gradient can be replaced by:

[0222]

[0223] and

[0224]

[0225] Where k m,e,h and F th,e,h Ω e High-fidelity stiffness matrix and thermal loads for fine meshes.

[0226] In order to verify the effectiveness of the method proposed in the present invention, the experimental results are presented and explained below.

[0227] Under the conditions of volume constraints and temperature constraints, topology optimization experiments for minimizing flexibility are carried out. Figure 4 , 5Figures 6 and 7 respectively show the optimal design, stress and temperature distribution of the porous structure obtained by optimization. It can be seen that a clear topological structure is obtained through optimization and the temperature requirements of the heat source are met. In addition, with the temperature constraint of the heat source, the maximum temperature of the structure can be controlled within the allowable value of 325°C. Figure 7 The iterative process of the topology optimization problem is shown, indicating that the method of the present invention can achieve stable convergence. The change of the target value compliance and the change trend of the material volume fraction show a synchronous fluctuation pattern (refer to Figure 7 Iteration history comparison) is consistent with the traditional minimum flexibility optimization based on volume constraints. Moreover, the method provided by the present invention is very efficient, and only 200 iterations are required to obtain the optimal design. From the given iteration history, it can be seen that the convergence condition is actually reached at the 50th iteration.

[0228] Subsequently, the optimal topology obtained was compared with the optimal topology obtained under isotropic conditions. The relevant data of the optimization results under isotropic and anisotropic conditions are shown in Figure 2. Figure 8 and 9 shown. Figure 8 The iterative history of flexibility under two conditions is given, which also shows that the structure generated by considering anisotropy has better physical properties. Fig. 9 ), it is found that the structure generated by considering anisotropy has a lower average temperature, which means that the structure optimized by considering anisotropy has better heat dissipation performance.

[0229] In summary, the method proposed in the present invention shows significant advantages in improving design efficiency, ensuring constraint satisfaction, optimizing structural performance, and adapting to complex material properties, providing a powerful tool for design optimization in related fields.

Claims

1. A method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures, characterized in that: include: 1) For a given model and design domain, generate and establish the parameters and data structures required for subsequent simulation and optimization; It includes constructing the background coarse grid, constructing the seed points and radius, modeling the porous structure based on the voronoi structure, and pre-calculating its stress field; 2) Calculating the density field, and performing thermodynamic simulation on the generated porous structure based on the background coarse grid and shape function according to the obtained density field; During the thermodynamic simulation, the parameter values ​​and sensitivities of the relevant objective functions and constraints are calculated; According to the obtained stress field, the stress matrix corresponding to each seed point is queried, and the stress matrix is ​​subjected to characteristic decomposition and extraction processing to obtain the weight matrix corresponding to the seed point. The objective function value and sensitivity characterizing the structural anisotropy are calculated for the porous structure. Using topology optimization, the control parameters are updated through the GCMMA method according to the sensitivity, thereby updating the porous structure and optimizing the material distribution. 3) Repeat step 2) until the control parameter change is less than the threshold. When the above process is completed, the final seed point and radius, as well as the topological information of the porous structure are obtained. The above structure is used to construct the final optimized porous foam based on the signed distance field.

2. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 1, characterized in that: In step 1), a background coarse grid is constructed by tetrahedron-based mesh generation, and the stress field is pre-calculated for the model using a voxel-based topology optimization method to obtain a stress field that guides anisotropic optimization.

3. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 1, characterized in that: In step 1), the model is abstracted into a beam structure network controlled by the position and radius of the seed point, and the seed points for constructing the porous structure are generated and initialized. The Voronoi foam structure is generated in the design domain through the constructed background coarse grid, seed points and radius, thereby modeling a porous structure composed of multiple beams with thickness.

4. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 3 is characterized in that: Construct implicit functions to represent the Voronoi foam structure. Open Voronoi foam is defined as the edges of the Voronoi cells confined inside the model with a certain radius. Mix the implicit functions of all beams introduced by the Voronoi edges together, or together with the shell, to obtain an implicit representation of the Voronoi foam defined as a density field. Specifically: Each Voronoi bubble Ω V (X,r) consists of two parts: the internal part and the boundary part; For the inner part, in the bounding box, all its edges are cut using the outer free shape, and the remaining j-th edge becomes a beam with a radius equal to the average of the adjacent seed points of edge j. The Voronoi edge corresponding to vertices v1 and v2 defines the beam φ, which consists of a radius of A cylinder with height ||v1-v2|| and two radii The implicit form of the beam φ is for: Where d(x,v1,v2) represents the minimum distance from point x to edges v1,v2. in a=v2-v1, b=x-v1, e=x-v2, The implicit representation of the entire open-cell Voronoi foam is obtained by combining the implicit functions of all beams after smoothing using the Kreisselmeier-Steinhauser (KS) function. Where φ max =max(φ1,…,φ n ), n is the number of beams, and p is 16; The boundary part associated with the external shape Ω is also implicitly represented by φ Ω (x) indicates; Voronoi Foam Omega V (x,r) is finally described as a density field by the regularized Heaviside function H(Φ(x)): where ∈ controls the size of the regularization, α = 1e -6 .

5. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 4, characterized in that: In step 2), in order to realize the implicit foam simulation, the model is immersed in a fixed tetrahedral background grid D = {D e ,e=1,2,…,N} in the covering shape Ω, that is, In each tetrahedral element D e Medium, density H e Set it as the average value of H(x) on its four nodes, In the formula, For unit D e The coordinates of the kth node in ; Consider a fixed displacement boundary Γ u and a fixed temperature boundary Γ T The design domain Ω, the surface force acts on the Neumann boundary Γ N The surface heat flux is located at the boundary Γ h superior; Then the static equilibrium equation of the linear thermoelastic continuum under mechanical and temperature loads is expressed as: K m U=F=F m +F th have K m =∫ Ω B T DBdΩ F th H∫ Ω B T Dε th dΩ Where K m is the stiffness matrix, U is the node displacement vector, F m is the mechanical load vector, F th is the heat load vector; B is the strain displacement matrix, N is the shape function matrix; ε th is the thermal strain vector; K m From the element stiffness matrix K m,e Combination K m =∑ e K m,e Where B e is the unit strain-displacement matrix, Ω e is the unit physical domain; Heat load F th By unit heat load F th.e Combination F th =∑ e F th,e e th =α(T e -T ref )F T Where ε th,e is the unit thermal strain vector, α is the thermal expansion coefficient, T e is the unit temperature, is the unit node temperature T e The average value, T ref is the thermal expansion reference temperature, Φ is the vector [1 1 1 0 0 0] in three-dimensional space, and similarly, the displacement u(x) in the background grid D is obtained by interpolating the node displacement U. u(x)=N e (coin e Where N e (x) is the shape function matrix on the unit node, U e is the unit node displacement; Incorporating the thermal stress generated in the material into the mechanical equation, the thermal balance equation for steady-state heat conduction is: K T (x,r)T=P Where K T is the thermal conductivity matrix, T is the node temperature vector, P is the heat flux density vector, K T and P is given by the unit thermal conductivity matrix K T,e and the unit heat flow vector P e Combination of K T (X,r)=∑ e K T,e (X,r) P=∑ e P e Where λ e is the effective thermal conductivity, By interpolating the node temperature T, we get the temperature T(x) in the background grid D. T(x)=N e (x)T e 。 6. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 5, characterized in that: The Voronoi foam design problem uses the seed x and radius r as the design variables for performance optimization, and measures the flexibility C(X,r,U) of the elastic potential of the object as the physical objective. C(X,r,U)=F T U The Voronoi cells are regularized into approximately regular polyhedrons through shape adjustment energy S(X) That is, each Voronoi cell v i Seeds X i With the centroid The Euclidean distance between The weighted sum; For the optimization problem with volume and temperature constraints, find the optimal design parameter p = {p j ,j=1,2,…,n}, where p j Represents a component of a variable x or r: The objective function J(p,U) is set as the weighted sum of the physical objective and the shape regularization term J(p,U)=(1-w)C(p,u)+w S(X), V and is the total volume of the structure and the maximum volume constraint, T i and are the node temperatures and the corresponding maximum temperature constraints, p j and is the design variable p j The lower and upper bounds of The optimization problem is solved by using the gradients of the optimization objective and constraints relative to the design variable p and a global convergent moving asymptote method (GCMMA) based on numerical gradient.

7. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 6, characterized in that: Compliance C versus design parameter p j The gradient of is: have in The temperature gradient with respect to the design parameters is 8. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 6, characterized in that: For each Voronoi cell v i The weight matrix at point x The calculation of is as follows: Based on the stress field obtained by voxel model simulation and optimization, the Cauchy stress tensor at point x in the stress field is where σ ii is the normal stress, σ ij , i≠j is the shear stress. Due to the symmetry of shear stress, σ ij =σ ji , so σ x For a symmetric matrix, decompose its eigenvalues: s x =VΛV -1 Where Λ is a diagonal matrix with eigenvalues ​​on the diagonal, V is a matrix with eigenvectors as its columns, and the eigenvector corresponding to the maximum eigenvalue is extracted as the principal direction vector v m : v m =[v1 v2 v3] M x By v m Components:

9. The method for mechanical-thermal coupling simulation and structural optimization of ultra-complex structures according to claim 6, characterized in that: The material-aware CBN shape function is introduced to replace the original shape function to reduce the degree of freedom and speed up the calculation. Specifically, for the entire virtual domain divided by the coarse grid, each coarsening unit e and its virtual domain Ω e , further divide it into a tetrahedral mesh to obtain the refined mesh h of e; The unit shape function N e Transformation matrix by CBN Replace the element shape functions with those composed of linear shape functions on a fine mesh: Where N e,h (x) is Ω e A collection of node shape functions on a fine grid; Strain displacement matrix B e Replace with Where B e,h Ω e Global strain displacement matrix of fine mesh; CBN transformation matrix Aims to map CBN node values ​​to Ω e The internal value in is obtained by the boundary interpolation matrix ψ and the boundary-interior transformation matrix The product of is as follows: where ψ and Map the displacement of the CBN to the boundary nodes and then to Ω e The fully refined nodes in . Similarly, in the simulation, the gradient of the element thermal conductivity matrix is ​​calculated as Similarly, K m,e ,F th,e and its gradient can be replaced by: and Where K m,e,h and F th,e,h Ω e High-fidelity stiffness matrix and thermal loads for fine meshes.

10. A mechanical-thermal coupling simulation and structural optimization system for ultra-complex structures, characterized in that: Used to implement the method according to any one of claims 1 to 9.

Images