A linear elastic structure topology optimization method and system

Abstract

The application discloses a linear elastic structure topology optimization method and system, and relates to the field of structure design.The method comprises the following steps: constructing a three-dimensional model of a target linear elastic structure according to size information and boundary information; discretizing the three-dimensional model by using a numerical flow method to obtain a gridded model; determining a stiffness matrix and a load matrix of the gridded model; constructing a level set interpolation point according to node information of the gridded model by using a parameterized level set method; initializing a level set value of the level set interpolation point, and constructing an initial topology structure of the target linear elastic structure according to the initialized level set value; constructing a topology optimization model based on the stiffness matrix and the load matrix, and taking the minimum flexibility as the target under the set volume constraint; and optimizing the initial topology structure by using the topology optimization model based on an iterative method to obtain an optimal topology structure of the target linear elastic structure.The application can obtain a high-precision topology structure of a linear elastic structure.

Description

Technical Field

[0001] This invention relates to the field of structural design, and in particular to a method and system for topology optimization of linear elastic structures. Background Technology

[0002] Topology optimization, as a structural design method, enables a structure to find the optimal structural shape that satisfies target performance within a given design space, while meeting certain constraints. Compared to other topology optimization methods, the parametric level set method can generate clear and smooth structural boundaries. However, the traditional parametric level set method is limited by the low continuity of finite elements and the high mesh requirements, resulting in significant limitations in computational accuracy and precise representation of complex boundaries. Therefore, achieving high-precision structural topology optimization has become a pressing problem. Summary of the Invention

[0003] Based on this, embodiments of the present invention provide a method and system for optimizing the topology of linear elastic structures, so as to improve the accuracy of calculation and the precision of boundary representation, thereby obtaining a high-precision topology of the linear elastic structure.

[0004] To achieve the above objectives, embodiments of the present invention provide the following solutions:

[0005] A method for topology optimization of linear elastic structures, comprising:

[0006] Obtain the dimensional and boundary information of the target linear elastic structure;

[0007] A three-dimensional model of the target linear elastic structure is constructed based on the size information and the boundary information.

[0008] The three-dimensional model is discretized using the numerical manifold method to obtain a meshed model; the meshed model comprises multiple numerical manifold elements.

[0009] Displacement boundary conditions are applied to the meshed model to obtain the stiffness matrix and load matrix of the meshed model;

[0010] The parametric level set method is used to construct level set interpolation points based on the node information of the gridded model; the node information includes the mathematical grid nodes, the intersection points of the mathematical grid and the physical grid, and the information of the physical domain vertices in the gridded model.

[0011] The horizontal set values ​​of the interpolation points of the horizontal set are initialized to obtain the initialized horizontal set values, and the initial topology of the target linear elastic structure is constructed based on the initialized horizontal set values.

[0012] Based on the stiffness matrix and load matrix, a topology optimization model is constructed with the goal of minimizing flexibility under the set volume constraints.

[0013] Based on an iterative method, the initial topology is optimized using the topology optimization model to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology is within a set flexibility range.

[0014] The present invention also provides a linear elastic structure topology optimization system, comprising:

[0015] The information acquisition module is used to acquire the size and boundary information of the target linear elastic structure;

[0016] The model building module is used to construct a three-dimensional model of the target linear elastic structure based on the size information and the boundary information.

[0017] The model discretization module is used to discretize the three-dimensional model using the numerical manifold method to obtain a meshed model; the meshed model includes multiple numerical manifold elements.

[0018] The stiffness and load matrix determination module is used to apply displacement boundary conditions to the meshed model to obtain the stiffness matrix and load matrix of the meshed model.

[0019] The interpolation point selection module is used to construct level set interpolation points based on the node information of the gridded model using the parametric level set method; the node information includes the mathematical grid nodes, the intersection points of the mathematical grid and the physical grid, and the information of the physical domain vertices in the gridded model.

[0020] The interpolation point initialization module is used to initialize the horizontal set values ​​of the horizontal set interpolation points, obtain the initialized horizontal set values, and construct the initial topology of the target linear elastic structure based on the initialized horizontal set values.

[0021] The optimization model building module is used to construct a topology optimization model based on the stiffness matrix and load matrix, with the goal of minimizing flexibility under set volume constraints.

[0022] The topology optimization module is used to optimize the initial topology using the topology optimization model based on an iterative method to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology is within a set flexibility range.

[0023] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0024] This invention proposes a method and system for topology optimization of linear elastic structures. It combines numerical manifold methods and parameterized level set methods to construct an initial topology. Based on the stiffness matrix and load matrix, and under set volume constraints, a topology optimization model is constructed with the goal of minimizing flexibility. Using an iterative method, the initial topology is optimized using this model to obtain the optimal topology of the target linear elastic structure. This invention achieves topology optimization of linear elastic structures based on numerical manifold methods, effectively improving the accuracy of computation and boundary representation, thereby obtaining a high-precision topology for the linear elastic structure. Attached Figure Description

[0025] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0026] Figure 1 A flowchart of a linear elastic structure topology optimization method provided in an embodiment of the present invention;

[0027] Figure 2 A schematic diagram of the level set;

[0028] Figure 3 A schematic diagram of the method for selecting interpolation points for a parameterized level set;

[0029] Figure 4 Comparison images before and after the boundary smoothing process;

[0030] Figure 5 A schematic diagram of a cut manifold unit located within the design domain;

[0031] Figure 6 A schematic diagram of a cut manifold element located at the boundary of the design domain;

[0032] Figure 7 This is a schematic diagram of the iterative process of a two-dimensional cantilever beam.

[0033] Figure 8 A schematic diagram of a complex design domain optimization problem;

[0034] Figure 9 This is a schematic diagram of a three-dimensional cantilever beam optimization problem.

[0035] Figure 10 This is a schematic diagram of a three-dimensional curved beam optimization problem. Detailed Implementation

[0036] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0037] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] Example 1

[0039] Numerical manifold methods, based on two sets of covering systems, are used for model discretization and computation, effectively improving computational accuracy and the precision of boundary representation. Therefore, this invention proposes numerical manifold elements and establishes a parameterized horizontally bundled elastic structure topology optimization method based on numerical manifold elements.

[0040] See Figure 1 The linear elastic structure topology optimization method of this embodiment includes:

[0041] Step 101: Obtain the size and boundary information of the target linear elastic structure.

[0042] Step 102: Construct a three-dimensional model of the target linear elastic structure based on the size information and the boundary information.

[0043] Step 103: Discretize the three-dimensional model using the numerical manifold method to obtain a meshed model; the meshed model includes multiple numerical manifold elements.

[0044] The global displacement function of the numerical manifold element is: U = Td; where U is the global displacement function, T is the covering system, and d is the generalized degree of freedom vector; the numerical manifold element is a two-dimensional numerical manifold element or a three-dimensional numerical manifold element.

[0045] Step 103 specifically includes: using a mathematical mesh to cut the physical domain of the three-dimensional model to obtain a meshed model.

[0046] Step 104: Apply displacement boundary conditions to the meshed model to obtain the stiffness matrix and load matrix of the meshed model.

[0047] Step 105: Using the parametric level set method, construct level set interpolation points based on the node information of the gridded model; the node information includes the information of mathematical grid nodes, the intersection of the mathematical grid and the physical grid, and the vertices of the physical domain in the gridded model.

[0048] Step 106: Initialize the level set values ​​of the level set interpolation points to obtain initialized level set values, and construct the initial topology of the target linear elastic structure based on the initialized level set values. Specifically, the method for determining the initialized level set values ​​is as follows:

[0049] (1) Assign initial values ​​to the horizontal set values ​​of the interpolation points of the horizontal set to obtain the original horizontal set values; a horizontal set value greater than 0 indicates that the corresponding interpolation point is a solid region, and a horizontal set value less than 0 indicates that the corresponding interpolation point is a void region.

[0050] (2) The original level set values ​​of the level set interpolation points in the region near the boundary are smoothed to obtain smoothed level set values; the region near the boundary includes the physical domain boundary and the region within a set range from the physical domain boundary. The smoothing process is as follows:

[0051] For the first type of level set interpolation point, it is processed according to the first condition to obtain the first smoothed level set value; the first type of level set interpolation point is the level set interpolation point on the physical domain boundary; the first condition includes: if the original level set value of the first type of level set interpolation point is greater than 0, then the original level set value of the first type of level set interpolation point is fixed to 0; if the original level set value of the first type of level set interpolation point is less than 0, then it remains unchanged.

[0052] For the second type of level set interpolation points, the second condition is applied to obtain the second smoothed level set value; the second type of level set interpolation points are level set interpolation points located outside the design domain within the region near the boundary; the second condition includes: determining the original level set value of the second type of level set interpolation points as the negative signed distance function of the nearest reference point; the nearest reference point is the level set interpolation point of the first type of level set interpolation points that is closest to the second type of level set interpolation point.

[0053] For the third type of level set interpolation points, the third condition is applied to obtain the third smoothed level set value. The third type of level set interpolation points are level set interpolation points located within the design domain in the region near the boundary. The third condition includes: if the original level set value of the third type of level set interpolation point is negative, it remains unchanged; if the original level set value of the third type of level set interpolation point is positive and the original level set value is greater than a positive signed distance function, it remains unchanged; if the original level set value of the third type of level set interpolation point is positive and the original level set value is less than a positive signed distance function, the level set value of the interpolation point is set to a positive signed distance function.

[0054] The first smoothed level set, the second smoothed level set, and the second smoothed level set constitute the smoothed level set.

[0055] (3) The original level set value is updated using the smoothed level set value to obtain the initialized level set value.

[0056] Step 107: Based on the stiffness matrix and load matrix, construct a topology optimization model with the objective of minimizing flexibility under the set volume constraints.

[0057] The topology optimization model is as follows:

[0058]

[0059] N is the number of level set interpolation points, α is the level set value matrix, α1 is the level set value of the first level set interpolation point, α2 is the level set value of the second level set interpolation point, α3 is the level set value of the third level set interpolation point, and α... i Let α be the level set value of the i-th level set interpolation point. N Let J(u, Φ) be the level set value at the Nth level set interpolation point, T denote the transpose, J(u, Φ) be the objective function, u be the allowable displacement (i.e., the real displacement field) in the displacement field U, Φ be the level set value, Ω be the design domain, ε be the strain field, and f, j, k, l represent different dimensional directions. fj (u) is used to calculate the displacement component u in the f direction. f The strain tensor, E, obtained by the partial derivative with respect to the j-direction. fjkl For the elastic modulus, ε kl (u) is used to calculate the displacement component u in the k direction. k The strain tensor obtained by the partial derivative with respect to the l-direction, H(Φ) represents the Heaviside function, G(Φ) is the volume constraint function, and V... m醸x Let v be the maximum allowable volume fraction, v be the virtual displacement, u0 be the displacement on the Dirichlet boundary (as a boundary condition), φ(u, v, Φ) = l(u, v) be the weak form of the elastic equilibrium equation, φ(u, v, Φ) be the bilinear energy form, and l(u, v) be the linear load form. For the structural boundary, α i，min For α i The upper limit, α i，m醸x For α i The lower limit.

[0060] Step 108: Based on the iterative method, the initial topology is optimized using the topology optimization model to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology is within the set flexibility range.

[0061] Step 108 specifically includes:

[0062] For the t-th iteration, the stiffness, volume fraction, load, and level set value of the t-th iteration are substituted into the topology optimization model to calculate whether the compliance of the t-th iteration meets the set convergence condition. The level set value of the 1st iteration is the initial level set value in the initial topology. The set convergence condition includes: the relative error of compliance between five adjacent iterations does not exceed 10. -3 .

[0063] If so, the level set value under the t-th iteration is determined as the optimal level set value, and the optimal topology of the target linear elastic structure is constructed from the optimal level set value.

[0064] If not, the moving asymptotic algorithm is used to update the level set value under the t-th iteration based on the sensitivity information to obtain the level set value under the (t+1)-th iteration; considering the influence of the shape of the actual integration region of the cut numerical manifold element in the meshed model on the stiffness under the t-th iteration, the stiffness under the t-th iteration is updated to obtain the stiffness under the (t+1)-th iteration; the area of ​​the cut numerical manifold element in the meshed model under the t-th iteration is compared with the area at the first iteration to obtain the volume fraction under the (t+1)-th iteration, and then the process proceeds to the (t+1)-th iteration.

[0065] The following section provides a detailed introduction to several aspects involved in the above-mentioned linear elastic structure topology optimization method.

[0066] I. Three-dimensional numerical manifold units

[0067] First, the global displacement function of this element can be constructed in the following form:

[0068]

[0069] Where u(x, y, z), v(x, y, z), and w(x, y, z) are the global displacement vectors defined on the numerical manifold element. u, v, and w are the nodal displacement vectors, u x u y u z Let N(x, y, z) be the partial derivative of the nodal displacement vector with respect to x, y, and z, and let N(x, y, z) be the shape function associated with the nodal displacement vector. x (x, y, z) is a shape function associated with the derivative of the nodal displacement vector with respect to x, N y (x, y, z) is a shape function associated with the derivative of the nodal displacement vector with respect to y, N z (x, y, z) is a shape function associated with the derivative of the nodal displacement vector with respect to z. (x, y, z) represents the coordinate information of the point in three-dimensional space.

[0070] As can be seen, the generalized degrees of freedom of the element at this time are [uu x u y u z At this point, we can combine the geometric equations from elasticity:

[0071]

[0072] Where, ε x Let ε be the strain component in the x-direction. y Let ε be the strain component in the y-direction. z γ is the strain component in the z-direction. xy γ represents the shear strain in the xy plane. xz γ represents the shear strain in the xz plane. yz ω represents the shear strain in the yz plane. z For the rotation angle information about the z-axis, ω y For the rotation angle information about the y-axis, ω x This provides information on the rotation angle around the x-axis.

[0073] By performing a simple transformation on the degrees of freedom, we can obtain the following formula:

[0074]

[0075] At this point, the global displacement function of the unit proposed in this invention will be transformed into the following form:

[0076] U=Td#(4)

[0077] Where U is the global displacement function, T is the covering system, and d is the generalized degree of freedom vector.

[0078] Among them, U=[u(x,y,z)v(x,y,z)w(x,y,z)] T And there are:

[0079]

[0080] d=[uvwε x ε y ε z γ xz γ xz γ yz ω x ω y ω z ] T #(6)

[0081] The expression for the shape function is:

[0082]

[0083] Where k represents the k-th node of the hexahedral element, ξ is the isoparametric coordinate in the x-direction, η is the isoparametric coordinate in the y-direction, ζ is the isoparametric coordinate in the z-direction, and N... k Let represent the shape function at the k-th node with respect to the nodal displacement vector. Let represent the shape function associated with the derivative of the nodal displacement vector with respect to ξ at the k-th node. Let represent the shape function associated with the derivative of the nodal displacement vector with respect to η at the k-th node. ξ represents the shape function associated with the derivative of the nodal displacement vector with respect to ζ at the k-th node. k η represents the value of ξ at the k-th node. k ζ represents the value of η at the k-th node. k Let ξ0 represent the value of ζ at the k-th node. k ξ, η0=η k η, ζ0=ζ0ζ.

[0084] As can be seen, the degree of freedom of the element has changed from the aforementioned displacement and partial derivatives of displacement to displacement and strain components with clear physical meaning. This means that the degree of freedom of the element no longer only serves to improve the accuracy of the element, but also greatly facilitates the calculation of mechanical parameters such as stress and strain.

[0085] II. Two-Dimensional Form of a Novel Hexahedral Numerical Manifold Unit

[0086] The three-dimensional form mentioned above, neglecting the degree of freedom in the z-direction, and utilizing the elasticity geometric equations of the two-dimensional problem, can be easily degenerated into the following form:

[0087]

[0088]

[0089]

[0090] The global displacement function can be written in the following form:

[0091] u=Td#(11)

[0092]

[0093] d=[u vε x ε y γ xy ω] T #(13)

[0094] The expression for the shape function is:

[0095]

[0096] In numerical manifold methods, since the boundaries of the mathematical mesh may not coincide with the boundaries of the physical domain, displacement boundary conditions cannot usually be applied directly. Currently, the Lagrange multiplier method and the penalty function method are commonly used. This embodiment will use the penalty function method to apply boundary conditions. Based on this, the weak form of the governing equations in structural analysis can be written as follows:

[0097]

[0098] Among them, Γ S For stress boundary, Γ d For displacement boundary, For Γ d Given the displacement above, and b as the body force, Given the surface force, k is a freely selectable penalty value, ε is the strain field, σ is the stress field, Ω is the design domain, S represents the design domain boundary, and u represents the actual displacement field.

[0099] Equation (15) can be transformed as follows:

[0100]

[0101]

[0102] ∫ Ω u T bdΩ=d T ∫ Ω T T bdΩ#(18)

[0103]

[0104] Where D is the elasticity matrix, matrix B = LT, and the specific expressions for D and L are:

[0105]

[0106]

[0107] in,

[0108]

[0109] Based on equations (17) to (19), equation (16) can be written in the following form:

[0110]

[0111] According to the variational principle, the functional Π P The condition for a stationary value is that its first variation is 0, that is:

[0112]

[0113] Therefore, the solution equations for the numerical manifold element can be obtained as follows:

[0114] Kd=P#(25)

[0115] Where K is the stiffness matrix and P is the load matrix, K and P are derived from the element stiffness matrix K e With element load matrix P e It was formed from, and is specifically defined as:

[0116]

[0117]

[0118] Unless otherwise specified, the penalty value k is generally taken as 10 in this paper. 4 E, where E is the elastic modulus.

[0119] III. Parameterized Level Set Method Based on New Numerical Manifold Units

[0120] (1) Parametric level set method

[0121] In topology optimization design methods based on the level set method, the core idea is to describe the evolution of the structural topology through the dynamic evolution of a level set function with a dimension one higher than the design domain. To better illustrate the concept of the level set method, a two-dimensional problem is used as an example. The boundary of a two-dimensional design domain D... The evolution will be described through the zero level plane of a three-dimensional level set function. For example... Figure 2 As shown, the projection of regions with level set function values ​​greater than 0 onto the zero level plane will represent the solid material region S in the design domain, while the projection of regions with level set function values ​​less than 0 onto the zero level plane will represent the void region D\S in the design domain. The relationship between level set function values ​​and material distribution is as follows:

[0122]

[0123] Where Φ(x) represents the level set value of the node with coordinate x.

[0124] Topology optimization of a structure is a dynamic evolutionary process. To give the level set function time characteristics, a dynamic virtual time factor t needs to be introduced into the level set function. At this time, the zero level surface can be expressed as Φ(x,t)=0. Taking the partial derivative of this equation with respect to t, we can obtain the following Hamilton-Jacobi equation:

[0125]

[0126] Among them, Vn The normal velocity field representing the design variables is in the following form:

[0127]

[0128]

[0129] Where V is the velocity field of the design variable.

[0130] For the traditional level set method, the main factor affecting its computational efficiency is solving a partial differential equation that couples time and space variables. To overcome this challenge, a structural topology optimization model based on radial basis function interpolation is introduced into the level set method. By interpolating the level set function using radial basis functions, the PDE partial differential equation can be transformed into an ODE partial differential equation with decoupled time and space variables. At this point, the level set function Φ(x, t) can be transformed into the following form:

[0131]

[0132] in, Let represent the radial basis functions, and α(t) be the expansion coefficient. After interpolation using the radial basis functions, the shape of the level set function can be determined by solving for the expansion coefficient α(t) in each iteration, simplifying the numerical solution.

[0133] Considering the high precision characteristics of numerical manifold elements, this invention selects compactly supported radial basis functions with C2 continuity, whose mathematical description is as follows:

[0134]

[0135] Taking a two-dimensional problem as an example, the expression for the support radius r is:

[0136]

[0137] Where d is the radial basis function at the interpolation point (x i y i The radius of influence at point () is the support radius. As an adjustable parameter in the radial basis functions, the support radius affects the final topology result. A larger support radius results in a simpler optimization, while a smaller support radius tends to lead to rougher boundaries in the final optimization result. Generally, the support radius is 2-5 times the mesh size.

[0138] (2) Solving using the parameterized level set method

[0139] Equation (32) can be written in the following matrix form:

[0140] Φ=Aα#(35)

[0141] Among them, Φ=[Φ(x1,t),Φ(x2,t),…,Φ(x N ,t)] T Let A be the level set function matrix, A be the compactly supported radial basis function interpolation matrix, and α be the coefficient matrix.

[0142]

[0143] α=[α1(t) α2(t) … α N (t)] T #(37)

[0144] Substituting equation (35) into equation (29) yields the following:

[0145]

[0146] in,

[0147]

[0148]

[0149]

[0150] Therefore, considering equation (35), equation (38) can be written in the following form:

[0151]

[0152] in,

[0153]

[0154] For equation (42), the solution obtained using the first-order Euler method is:

[0155]

[0156] α i+1 =α i +Δt·A -1 B(α,t)#(45)

[0157] Where Δt is the time step. To avoid unbounded growth of the level set function, the level set function is reinitialized in the following way:

[0158]

[0159] Where, Φ new This represents the level set function after reinitialization. Let the level set value be the value at the r-th interpolation point near the boundary. Let be the magnitude of the gradient of the level set at the interpolation point. Simultaneously, the Dirac function δ(Φ) is introduced into equation (45) for constraint, and equation (45) can then be rewritten as follows:

[0160]

[0161] in,

[0162]

[0163] Here, δ(Φ) is an approximation function used to avoid unbounded growth.

[0164]

[0165] As can be seen from equation (38), in the parameterized level concentration, the velocity field V n In optimization calculations, it is naturally extended to the entire design domain, which enables the MMA algorithm based on the gradient information of the objective function to be applied in the parametric level set method.

[0166] (3) Topology optimization model for minimum compliance problem

[0167] In this invention, all problems are addressed with the objective function of minimizing structural compliance under certain volume constraints. Furthermore, the expansion coefficient α is used as the design variable in the parameterization level set; therefore, the optimization model can be constructed as follows:

[0168]

[0169] Where N is the number of interpolation points, J(u, Φ) is the objective function, ε is the strain field, H(Φ) represents the Heaviside function, and V max Let α be the maximum allowable volume fraction, u represent the allowable displacement in the displacement field U, v represent the virtual displacement, and u0 be the displacement on the Dirichlet boundary. G(Φ) is the volume constraint function, a(u, v, Φ) = l(u, v) is the weak form of the elastic equilibrium equation, and α i，min With α i，max The upper and lower bounds of the design variables are defined. Meanwhile, the expressions for H(Φ), a(u, v, Φ), and l(v, Φ) are as follows:

[0170]

[0171]

[0172]

[0173] In equation (53), Ω represents the entire design domain, Γ is the structural boundary, and b is the body force. For traction force.

[0174] In the minimum compliance problem, the expression for sensitivity analysis is as follows:

[0175]

[0176] Where, the expression for β1 is ξ is the sensitivity factor.

[0177] After obtaining the above sensitivity information, the design variables can be updated using the Moving Asymptotes (MMA) algorithm.

[0178] (4) Interpolation point selection method

[0179] For parametric level sets, the design variables are associated with interpolation points covering the entire design domain. The sign of the level set value at the interpolation point determines which areas in the problem domain are solid regions and which are void regions, and the visualization of the design results also requires information about the interpolation points. In traditional finite element-based parametric level sets, interpolation points are often directly selected as nodes of the finite element. In this invention, numerical manifold elements are used to replace finite elements to complete the discretization and structural analysis of the design domain. A major difference between numerical manifold elements and finite elements is that numerical manifold elements use two non-overlapping coverage systems. In the parametric level set method based on the novel numerical manifold element, such as... Figure 3 As shown, the interpolation point consists of nodes of the mathematical grid, intersections of the mathematical grid and the physical domain, and vertices of the physical domain.

[0180] In some cases, the mathematical grid will not perfectly coincide with the boundary of the design domain, such as when the boundary of the design domain is a curve. In this case, for interpolation points, not all points fall within the design domain; some points will lie outside the boundary. The level set values ​​of these points are fixed as negative numbers and do not participate in the iterative process, but they still serve as reference points for interpolating the level set function using radial basis functions. As mentioned earlier, the basic principle of the level set method is to represent the optimization process of the structure using the evolution of a higher-dimensional level set function, and the visualization of the topology lies in obtaining the zero level plane of the level set function. However, as... Figure 4 As shown, if the horizontal set values ​​at interpolation points near the structural boundary are not smoothed, the resulting structure will have jagged boundaries, making it impossible to accurately capture the boundaries of the design domain. Figure 4 Part (a) represents the horizontal set surface before smoothing. Figure 4 Part (b) represents the model obtained before smoothing. Figure 4 Part (c) represents the horizontal set surface after smoothing. Figure 4The (d) part represents the smoothed model. Therefore, in this invention, the following method is used to smooth the level set values ​​of the interpolation points near the boundary:

[0181] (1) Based on the positional relationship between the interpolation point and the design domain, the interpolation points are divided into two categories: interpolation points located inside the design domain and interpolation points located outside the design domain.

[0182] (2) Obtain the interpolation point information located on the physical domain boundary. If the level set of these interpolation points is greater than 0, then fix the level set value of these interpolation points to 0; if the level set value of the interpolation points is less than 0, then keep the level set value unchanged.

[0183] (3) Taking the interpolation point located on the boundary as the reference point, for the interpolation point located outside the design domain, its level set value is selected as the signed distance function from the interpolation point to its nearest reference point, and fixed as a negative value.

[0184] (4) For interpolation points located inside the design domain, if their level set value is negative, they remain unchanged. If their level set value is positive, they need to be divided into two cases: if their original level set value is greater than the positive sign distance function, the original level set value remains unchanged; if their original level set value is less than the positive sign distance function, the level set value of the interpolation point is set to the sign distance function.

[0185] (5) Stiffness update scheme and volume fraction calculation

[0186] During the dynamic evolution of the level set function, it is inevitable that elements passing through the zero level set surface will be cut. Therefore, the stiffness matrix of the element needs to be recalculated during the iteration process based on the element cutting situation. In level set topology optimization based on the finite element method, the most commonly used stiffness update method is the "stiffness reduction" method. The "stiffness reduction" method does not consider the shape of the integration region of the cut element. Instead, it subdivides the mesh inside the element, determines the level set value of the subdivided mesh node by linear interpolation, and then statistically analyzes the level set value of each subdivided mesh node to roughly obtain the ratio of the area of ​​the cut element to the area of ​​the original element. Then, it uses linear reduction to obtain the stiffness information of the cut element. However, this method does not consider the influence of the actual integration region shape of the cut element on the stiffness, which will lead to serious distortion of the boundary representation of the structure and affect the accuracy of the objective function calculation.

[0187] Considering the limitations of the "stiffness reduction" method, this invention employs a novel element stiffness update scheme. Before performing the stiffness update calculation, each manifold element needs to be categorized into elements located within the design domain and elements containing nonlinear boundaries, based on their position within the design domain.

[0188] For manifold elements located within the design domain, their shape completely overlaps with the mathematical mesh blocks. Therefore, multi-level meshing techniques are not required for stiffness updates of this type of element. Based on the level set values ​​of each node of the manifold element, this type of segmented element can be divided into, for example: Figure 4 The five cases are shown. If the level set of all nodes of the element is greater than 0, then the element can be considered not to be cut by the boundary, and its stiffness does not need to be recalculated (e.g., ...). Figure 5 (As shown in part (a)). Similarly, if the level set values ​​of all nodes are less than 0, it can be assumed that the element is now entirely outside the structure due to boundary evolution. To prevent singularities in the overall stiffness matrix, it is necessary to use an elastic modulus of E. min =10 -4 E's weak material fills the element (e.g., Figure 5 (as shown in part (b)). When the level set values ​​of the nodes in a manifold element are inconsistent in sign, the result depends on the number of nodes with positive level set values, such as... Figure 5 Part (c) to Figure 5 As shown in part (e), there are three cases. For these three cases, it is necessary to first find the zero level set value point on the boundary, and then perform simplex integration on the polygon formed by this point and the points with positive level set values ​​to obtain the stiffness information of the cut element.

[0189] For manifold elements containing nonlinear boundaries, such as Figure 5 As shown, stiffness updates can be divided into the following five cases. If the level set values ​​of all nodes in the manifold element are positive (e.g., ... Figure 6 (as shown in part (a)) or both are negative (as shown in part (a)). Figure 6 As shown in part (b), the method for calculating stiffness is the same as that for elements located inside the design domain. When the horizontal set values ​​of the element nodes are inconsistent in sign, there are three possible cases.

[0190] First, such as Figure 6 As shown in section (c), the level set values ​​of the two nodes located on the nonlinear boundary are both positive, while the nodes with negative level set values ​​are located on the linear element boundary. In this case, it is only necessary to find the zero level set value point on the linear boundary by linear interpolation, and then obtain the stiffness information of the cut element by the simplex integral method based on multi-level subdivision technique.

[0191] Second, such as Figure 6 As shown in part (d), if all the level set values ​​of the nodes on the boundary are negative, then the multi-level partitioned region containing the nonlinear boundary needs to be removed from the integral region of the element. At this time, the remaining integral region is a polygon, and the stiffness information can be obtained directly using the traditional simplex integral.

[0192] Finally, as Figure 6 As shown in part (e), if the level set values ​​of the nodes on the boundary are all different, there will be a zero level set value point on the nonlinear boundary, and the multilevel partitioned region will be cut off by the structural boundary. To obtain a new multilevel partitioned region, it is necessary to first find the projection point of the zero level set value point on the line connecting the two nodes. Draw a perpendicular line through this projection point and intersect it with the nonlinear boundary. This intersection point is the zero level set value point located on the nonlinear boundary. After obtaining this zero level set value point, replace the node with the negative level set value with this point to form a new multilevel partitioned region. Use multilevel partitioning technology to obtain the integral information of the new multilevel partitioned region.

[0193] In the parametric level set method based on finite elements, the calculation of volume fraction is similar to the calculation of stiffness. Based on the refined mesh, the volume fraction of the cut element is determined according to the proportion of refined nodes with positive level set values ​​among the total refined nodes. In this invention, since simplex integrals are used for stiffness updates, and the area calculation of the integration region is already completed in the simplex integral, there is no need to calculate the volume fraction again. The volume fraction of the element is obtained by comparing the area of ​​the cut element calculated during the iteration process with its initial area.

[0194] (6) Verification by example

[0195] Based on the above work, in order to more fully illustrate the features of the present invention and its applicability to practical problems, the following two-dimensional and three-dimensional calculation examples will be used for verification. First, as... Figure 7 The example shown is a two-dimensional cantilever beam. The dimensional information of the example is as follows: Figure 7 As shown in part (a), the volume constraint is set to 50%, the material parameters are E=1, and the elastic modulus of the material in the void region is E min =10 -4 Poisson's ratio is set to v = 0.3. From Figure 7 As can be seen in part (b), the objective function and the volume fraction reach convergence at step 64. Figure 7 Part (c) shows the structures corresponding to the initial iteration, iteration step 3, iteration step 10, iteration step 30, iteration step 40, and final iteration. The final optimized structure obtained by this invention is very similar to the optimized structures obtained in classic literature, which also verifies the effectiveness of the parameterized level set method based on numerical manifold elements. It is also noted that the convergence rate of this method is faster than that of the parameterized level set method based on finite elements.

[0196] Meanwhile, to verify the applicability of the parameterized level set method based on numerical manifold elements to problems in complex design domains, for example... Figure 8 The two problems with complex nonlinear boundaries shown were optimized. From Figure 8As can be seen, the parameterized level set method based on numerical manifold units can efficiently obtain stable topology optimization results, while also obtaining smooth structural boundaries and accurately capturing the geometric features of nonlinear boundaries. Figure 8 Part (a) shows the dimensions of the curved beam. Figure 8 Part (b) shows the optimized curved beam structure. Figure 8 Section (c) shows the wrench size information. Figure 8 Section (d) shows the optimized wrench structure.

[0197] Finally, as Figure 9 and Figure 10 As shown, the application of novel numerical manifold elements in three-dimensional problems was investigated, and the results were validated using a three-dimensional cantilever beam example and a three-dimensional curved beam example. Comparison with results from classical literature reveals that the parameterized level set method based on novel numerical manifold elements can obtain effective topology optimization results, generate smooth boundaries, and has a faster convergence speed compared to the parameterized level set method based on finite elements, demonstrating significant advantages. Figure 9 Part (a) shows the dimensional information of the three-dimensional cantilever beam. Figure 9 Part (b) illustrates the optimization iterative process for the three-dimensional cantilever beam. Figure 9 Section (c) shows the optimization results for the three-dimensional cantilever beam. Figure 10 Part (a) shows the dimensional information of the three-dimensional curved beam. Figure 10 Part (b) illustrates the optimization iterative process of the three-dimensional curved beam. Figure 10 Section (c) shows the optimization results of the three-dimensional curved beam.

[0198] Based on the above, in practical applications, a specific implementation process of the above method is as follows:

[0199] Step 1: Based on the dimensional information of the structure to be analyzed and the mathematical characteristics of its boundaries, establish a model of the problem to be solved, i.e., a three-dimensional model.

[0200] Step 2: Based on the numerical manifold units described above, complete the discretization of the 3D model.

[0201] Step 3: Based on the characteristics of the numerical manifold covering system, read the node information, including mathematical grid nodes, intersection points of the mathematical and physical grids, and physical domain vertices, to obtain the interpolation points of the level set, and initialize the level set values ​​of the interpolation points. This step is specifically implemented using the interpolation point selection method described above.

[0202] Step 4: Based on the above steps, according to the loading conditions and boundary conditions of the model, obtain the corresponding load matrix and stiffness matrix and apply boundary conditions to perform structural analysis and solve (Equations (26) and (27)).

[0203] Step 5: Using the stiffness update scheme and volume fraction calculation method proposed above, the stiffness information and volume fraction are accurately updated during the optimization iteration process.

[0204] Step 6: Based on the results obtained from the structural analysis, perform sensitivity analysis (Equation (54)), and use the MMA algorithm to update the design variables in combination with the sensitivity information. Finally, after meeting the convergence condition, the design variables converge to a stable optimized structure.

[0205] Step 7: Based on the above structure, compare the results with those of traditional finite element-based topology optimization to study the performance and advantages of the optimization method proposed in this invention in the minimum flexibility design problem of linear elastic structures; and study the influence of factors such as mesh size and initial voids on the final optimized structure.

[0206] The method in this embodiment establishes the complete flow of the optimization algorithm, fully leveraging the advantages of numerical manifold methods and parameterized level sets, while overcoming the shortcomings of existing methods. Combining the characteristics of the novel numerical manifold method, a new element stiffness update method is proposed, which can effectively improve the accuracy of element stiffness calculation during iteration and achieve high fidelity for element boundaries. This novel parameterized level set method converges faster and is more efficient than the parameterized level set method based on finite elements. Based on the advantages and characteristics of this method, the optimization algorithm proposed in this invention has great application prospects in dynamic optimization, multi-objective optimization, and microstructure design.

[0207] Specifically, combining the numerical manifold method with the parametric level set method fully leverages the advantages of both methods and overcomes the shortcomings of the finite element-based parametric level set method. Stiffness updates are achieved using multi-level partitioning techniques and simplex numerical integration, which, compared to the traditional "stiffness reduction" method, can accurately capture stiffness information and volume fraction changes during the iteration process. By combining the characteristics of the numerical manifold method covering the system, interpolation point selection is implemented, unifying the interpolation points of the parametric level set with the nodes of the numerical manifold covering system. Two-dimensional and three-dimensional optimization problems based on the parametric level set method using numerical manifold elements are realized. The examples demonstrate that the parametric level set method based on numerical manifold elements effectively utilizes the advantages of high accuracy of numerical manifold elements, efficient covering systems, and accurate boundary description, thus improving the applicability of the parametric level set method.

[0208] Example 2

[0209] In order to implement the method corresponding to Embodiment 1 above and achieve the corresponding functions and technical effects, a linear elastic structure topology optimization system is provided below.

[0210] The system includes:

[0211] The information acquisition module is used to acquire the size and boundary information of the target linear elastic structure.

[0212] The model building module is used to construct a three-dimensional model of the target linear elastic structure based on the size information and the boundary information.

[0213] The model discretization module is used to discretize the three-dimensional model using the numerical manifold method to obtain a meshed model; the meshed model includes multiple numerical manifold elements.

[0214] The stiffness and load matrix determination module is used to apply displacement boundary conditions to the meshed model to obtain the stiffness matrix and load matrix of the meshed model.

[0215] The interpolation point selection module is used to construct level set interpolation points based on the node information of the gridded model using the parametric level set method; the node information includes the information of mathematical grid nodes, the intersection of the mathematical grid and the physical grid, and the vertices of the physical domain in the gridded model.

[0216] The interpolation point initialization module is used to initialize the horizontal set values ​​of the horizontal set interpolation points, obtain the initialized horizontal set values, and construct the initial topology of the target linear elastic structure based on the initialized horizontal set values.

[0217] The optimization model construction module is used to construct a topology optimization model based on the stiffness matrix and load matrix, with the goal of minimizing flexibility under set volume constraints.

[0218] The topology optimization module is used to optimize the initial topology using the topology optimization model based on an iterative method to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology is within a set flexibility range.

[0219] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the descriptions are relatively simple; relevant parts can be referred to the method section.

[0220] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method of linear-elastic structural topology optimization, characterized by, include: Obtain the dimensional and boundary information of the target linear elastic structure; A three-dimensional model of the target linear elastic structure is constructed based on the size information and the boundary information. The three-dimensional model is discretized using the numerical manifold method to obtain a meshed model; the meshed model comprises multiple numerical manifold elements. Displacement boundary conditions are applied to the meshed model to obtain the stiffness matrix and load matrix of the meshed model; The parametric level set method is used to construct level set interpolation points based on the node information of the gridded model; the node information includes the mathematical grid nodes, the intersection points of the mathematical grid and the physical grid, and the information of the physical domain vertices in the gridded model. The horizontal set values ​​of the interpolation points of the horizontal set are initialized to obtain the initialized horizontal set values, and the initial topology of the target linear elastic structure is constructed based on the initialized horizontal set values. Based on the stiffness matrix and load matrix, a topology optimization model is constructed with the goal of minimizing flexibility under the set volume constraints. Based on an iterative method, the initial topology is optimized using the topology optimization model to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology satisfies the set iterative convergence condition.

2. The method of topology optimization of linearly elastic structures according to claim 1, characterized in that, Based on an iterative method, the initial topology is optimized using the aforementioned topology optimization model to obtain the optimal topology of the target linear elastic structure, specifically including: For the t-th iteration, the stiffness, volume fraction, load, and level set value under the t-th iteration are substituted into the topology optimization model to calculate whether the compliance of the t-th iteration satisfies the set convergence condition; the level set value under the 1st iteration is the initial level set value in the initial topology. If so, the level set value under the t-th iteration is determined as the optimal level set value, and the optimal topology of the target linear elastic structure is constructed from the optimal level set value; If not, the moving asymptotic algorithm is used to update the level set value under the t-th iteration based on the sensitivity information to obtain the level set value under the (t+1)-th iteration; considering the influence of the shape of the actual integration region of the cut numerical manifold element in the meshed model on the stiffness under the t-th iteration, the stiffness under the t-th iteration is updated to obtain the stiffness under the (t+1)-th iteration; the area of ​​the cut numerical manifold element in the meshed model under the t-th iteration is compared with the area at the first iteration to obtain the volume fraction under the (t+1)-th iteration, and then the process proceeds to the (t+1)-th iteration.

3. The method of topology optimization of linearly elastic structures according to claim 1, characterized in that, The level set values ​​of the interpolation points of the level set are initialized to obtain initialized level set values, specifically including: Initial values ​​are assigned to the horizontal set values ​​of the interpolation points to obtain the original horizontal set values; a horizontal set value greater than 0 indicates that the corresponding interpolation point is a solid region, and a horizontal set value less than 0 indicates that the corresponding interpolation point is a void region. The original horizontal set values ​​of the horizontal set interpolation points in the region near the boundary are smoothed to obtain smoothed horizontal set values; the region near the boundary includes the physical domain boundary and the region within a set range from the physical domain boundary. The original level set value is updated using the smoothed level set value to obtain the initialized level set value.

4. The method of topology optimization of linearly elastic structures according to claim 3, characterized in that, The original level set values ​​in the region near the boundary are smoothed to obtain smoothed level set values, specifically including: For the first type of level set interpolation point, it is processed according to the first condition to obtain the first smoothed level set value; the first type of level set interpolation point is the level set interpolation point on the physical domain boundary; the first condition includes: if the original level set value of the first type of level set interpolation point is greater than 0, then the original level set value of the first type of level set interpolation point is fixed to 0; if the original level set value of the first type of level set interpolation point is less than 0, then it remains unchanged; For the second type of level set interpolation points, the second condition is applied to obtain the second smoothed level set value; the second type of level set interpolation points are level set interpolation points located outside the design domain within the region near the boundary; the second condition includes: determining the original level set value of the second type of level set interpolation points as the negative signed distance function of the nearest reference point; the nearest reference point is the level set interpolation point in the first type of level set interpolation points that is closest to the second type of level set interpolation point; For the third type of level set interpolation points, the third condition is applied to obtain the third smoothed level set value. The third type of level set interpolation points are level set interpolation points located within the design domain in the region near the boundary. The third condition includes: if the original level set value of the third type of level set interpolation point is negative, it remains unchanged; if the original level set value of the third type of level set interpolation point is positive and the original level set value is greater than a positive signed distance function, it remains unchanged; if the original level set value of the third type of level set interpolation point is positive and the original level set value is less than a positive signed distance function, the level set value of the interpolation point is set to a positive signed distance function. The first smoothed level set, the second smoothed level set, and the third smoothed level set constitute the smoothed level set.

5. The method of topology optimization of linearly elastic structures according to claim 1, characterized in that, The topology optimization model is as follows: N is the number of level set interpolation points, α is the level set value matrix, α1 is the level set value of the first level set interpolation point, α2 is the level set value of the second level set interpolation point, α3 is the level set value of the third level set interpolation point, and α... i Let α be the level set value of the i-th level set interpolation point. N Let J(u, Φ) be the level set value at the Nth level set interpolation point, T denote the transpose, J(u, Φ) be the objective function, u be the allowable displacement in the displacement field U, Φ be the level set value, Ω be the design domain, and ε be the strain field. fj (u) is used to calculate the displacement component u in the f direction. f The strain tensor, E, obtained by the partial derivative with respect to the j-direction. fjkl For the elastic modulus, ε kl (u) is used to calculate the displacement component u in the k direction. k The strain tensor obtained by the partial derivative with respect to the l-direction, H(Φ) represents the Heaviside function, G(Φ) is the volume constraint function, and V... max Let v be the maximum allowable volume fraction, v be the virtual displacement, u0 be the displacement on the Dirichlet boundary, a(u, v, Φ) = l(u, v) be the weak form of the elastic equilibrium equation, a(u, v, Φ) be the bilinear energy form, and l(u, v) be the linear load form. For the structural boundary, α i，min For α i The upper limit, α i，max For α i The lower limit.

6. The linear elastic structure topology optimization method according to claim 1, characterized in that, The global displacement function of the numerical manifold element is: U = Td; U is the global displacement function, T is the covering matrix, and d is the generalized degree of freedom vector; the numerical manifold element is a two-dimensional numerical manifold element or a three-dimensional numerical manifold element.

7. The linear elastic structure topology optimization method according to claim 1, characterized in that, The three-dimensional model is discretized using the numerical manifold method to obtain a meshed model, specifically including: The physical domain of the three-dimensional model is divided using a mathematical mesh to obtain a meshed model.

8. A linear elastic structural topology optimization system, characterized by, include: The information acquisition module is used to acquire the size and boundary information of the target linear elastic structure; The model building module is used to construct a three-dimensional model of the target linear elastic structure based on the size information and the boundary information. The model discretization module is used to discretize the three-dimensional model using the numerical manifold method to obtain a meshed model; the meshed model includes multiple numerical manifold elements. The stiffness and load matrix determination module is used to apply displacement boundary conditions to the meshed model to obtain the stiffness matrix and load matrix of the meshed model. The interpolation point selection module is used to construct level set interpolation points based on the node information of the gridded model using the parametric level set method; the node information includes the mathematical grid nodes, the intersection points of the mathematical grid and the physical grid, and the information of the physical domain vertices in the gridded model. The interpolation point initialization module is used to initialize the horizontal set values ​​of the horizontal set interpolation points, obtain the initialized horizontal set values, and construct the initial topology of the target linear elastic structure based on the initialized horizontal set values. The optimization model building module is used to construct a topology optimization model based on the stiffness matrix and load matrix, with the goal of minimizing flexibility under set volume constraints. The topology optimization module is used to optimize the initial topology using the topology optimization model based on an iterative method to obtain the optimal topology of the target linear elastic structure; the flexibility of the optimal topology is within a set flexibility range.

Images