Topological optimization method of elastic buckling multistable energy absorption structure

Abstract

The invention discloses a topological optimization method of an elastic buckling multistable energy absorption structure, and relates to the technical field of topological optimization. The topological optimization method of the elastic buckling multistable energy absorption structure comprises the steps that a topological optimization finite element model is initialized; updating the filtering density, the physical density and the material attribute of each unit in the established finite element model according to the current design variable; finite element calculation of a front buckling stage is carried out on the updated finite element model, a critical buckling load and sensitivity thereof are calculated in a front critical point state, and a first critical point corresponds to loading displacement and sensitivity thereof; carrying out subsequent finite element calculation in a post-buckling stage until the load reaches a critical buckling load, and calculating the loading displacement at the moment and the sensitivity of the loading displacement; calculating an objective function and corresponding sensitivity thereof, and a constraint function and corresponding sensitivity thereof; and constructing an explicit convex approximation optimization sub-problem by adopting a moving asymptote method, and solving through a Lagrange duality method to obtain an updated design variable until optimization convergence is completed.

Description

Technical Field

[0001] The present invention relates to the technical field of topology optimization, and in particular to a topology optimization method for an elastic buckling multi-stable energy-absorbing structure. Background Art

[0002] Energy-absorbing structures dissipate the kinetic energy of the collision body in various ways, thereby preventing the protected body from being damaged. They are widely used in various protective structures, such as vehicle bumpers, packaging for precision instruments, and personal protective equipment. Currently, common energy dissipation methods include metal plastic deformation, fiber or ceramic fracture damage, and material viscosity effect. These energy dissipation methods all have their own limitations. For example, energy dissipation structures based on plastic deformation or fracture damage can only be used once, while structures designed based on viscosity effects can be reused, but their mechanical response is too sensitive to the loading rate, resulting in unsatisfactory energy absorption effects for low-speed / low-frequency problems.

[0003] As an energy-absorbing structure, a multistable structure formed by an array of multiple elastic buckling bistable structures has the advantage of being reusable compared to traditional energy-absorbing structures. However, the energy absorption efficiency of the elastic buckling multistable structure is currently low, which limits its practical application. Therefore, it is of great significance to design an elastic buckling multistable structure with higher energy absorption efficiency.

[0004] Therefore, it is desired to have a technical solution to overcome or at least alleviate at least one of the above-mentioned defects of the prior art. Summary of the invention

[0005] The purpose of the present application is to provide a topology optimization method for an elastic buckling multi-stable energy-absorbing structure to overcome or at least alleviate at least one of the above-mentioned defects of the prior art.

[0006] To achieve the above objectives, the present application provides a topology optimization method for an elastic buckling multi-stable energy absorbing structure, the topology optimization method for an elastic buckling multi-stable energy absorbing structure comprising: Step 1: Initialize the topology optimization finite element model: given the design domain and material properties, discretize the design domain into a finite element mesh, initialize the design variables, couple the degrees of freedom of the loading area, and apply corresponding boundary conditions to the finite element model; Step 2: Update the filter density, physical density, and material properties of each unit in the finite element model established in step 1 according to the current design variables; Step 3: Perform finite element calculation of the pre-buckling stage on the finite element model updated in step 2 according to the internal force vector and each stiffness matrix, and calculate the critical buckling load, critical buckling load sensitivity, the first critical point corresponding loading displacement and the first critical point corresponding sensitivity under the pre-critical point state; Step 4: After obtaining the critical buckling load, the critical buckling load sensitivity, the loading displacement corresponding to the first critical point, and the sensitivity corresponding to the first critical point, perform finite element calculation in the post-buckling stage until the load reaches the critical buckling load, and calculate the loading displacement and loading displacement sensitivity when the critical buckling load is reached; Step 5: Calculate the optimization objective function, the constraint function, the sensitivity corresponding to the optimization objective function, and the sensitivity corresponding to the constraint function according to the information obtained in step 3 and step 4; Step 6: Use the moving asymptote method to construct the optimization subproblem and solve it to obtain the updated design variables; Step 7: Determine whether the topology optimization has converged. If so, the optimization ends. Otherwise, the updated design variables are used as current design variables and steps 2 to 6 are repeated until the topology optimization is determined to have converged.

[0007] Optionally, when performing the finite element calculation, the optimization objective function and the constraint function adopt the following optimization formula: ;in, is the theoretical maximum energy absorption per bistable unit; is the design variable; is the finite element structure displacement vector; is the critical buckling load, is the first critical point corresponding to the loading displacement, is the loading displacement when the load is the critical buckling load in the post-buckling stage; It is an unbalanced force; , The units are Design variables and material volumes; is the total volume of the design domain, is the maximum specified material volume fraction.

[0008] Optionally, in the finite element calculation, the unbalanced force The expression is: ;in, is the node internal force vector when the degree of freedom coupling and displacement boundary conditions are not considered, is the node external force vector when the degree of freedom coupling and displacement boundary conditions are not considered; matrix Defined as: ;in, is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; yes The vector of degrees of freedom that are not constrained by displacement boundary conditions; matrix Defined as: ;in, is the finite element structure displacement vector; is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; The internal force vector considering the degree of freedom coupling and displacement boundary conditions Use the following formula to obtain: ;in, The stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions includes a shear stiffness matrix, a linear stiffness matrix, and a stress stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions. as follows: ; The stress stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the stress stiffness matrix without considering degree of freedom coupling and displacement boundary conditions; The linear stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the linear stiffness matrix without considering degree of freedom coupling and displacement boundary conditions.

[0009] Optionally, the filtering density of each unit in step 2 is obtained by the following method: Through the filtering method, according to the design variables Calculate the filter density ; Through the smooth Heaviside projection function, according to the filtered density Calculate the physical density used for interpolation of element material properties .

[0010] Optionally, the node internal force vector is obtained by assembling the unit internal force vectors of each unit. I The unit internal force vector Use the following formula to obtain: ;in, is the displacement used for the nonlinear part of the strain energy density calculation, is the element internal force vector calculated using the large deformation theory, is the element internal force vector calculated under the small deformation assumption; The tangent stiffness matrix without considering the degree of freedom coupling and displacement boundary conditions is assembled from the element stiffness matrix of each element. I The element stiffness matrix Use the following formula to obtain: ;in, is the element stiffness matrix calculated using the large deformation theory, is the element stiffness matrix calculated under the small deformation assumption.

[0011] Optionally, performing finite element calculation of the finite element model updated in step 2 in the pre-buckling stage according to the internal force vector and each stiffness matrix, and calculating the critical buckling load, the critical buckling load sensitivity, the loading displacement corresponding to the first critical point, and the sensitivity corresponding to the first critical point in the pre-critical point state includes: In step 3, the equilibrium state of the previous incremental step needs to be saved during finite element calculation, and the previous critical point is continuously identified. When it is identified that the current equilibrium state has passed the first critical point, the current equilibrium state is saved, the previous equilibrium state is recorded as the previous critical point, and the finite element structure is restored to the state of the previous critical point. Under the state of the previous critical point, the critical buckling load, critical buckling load sensitivity, loading displacement corresponding to the first critical point, and sensitivity corresponding to the first critical point are calculated through the following steps: (1) The one-point method is used as the calculation method for nonlinear buckling load, and the first critical point is , the generalized eigenvalues ​​are solved by the subspace iteration method of Rayleigh-Ritz projection: ;in, is the stress stiffness matrix considering degree of freedom coupling and displacement boundary conditions, is the linear stiffness matrix considering degree of freedom coupling and displacement boundary conditions; is the modal eigenvalue ( ), is the corresponding eigenmode, represents the order of the buckling mode; (2) After solving the generalized eigenvalue problem, the buckling load is calculated using the following formula: ;in, is the buckling load, is the external force applied on the controlled degree of freedom in the pre-critical state, is the order of the buckling mode chosen to approximate the maximum eigenvalue, is the KS aggregation function; (3) Select the unit from the global degrees of freedom of the structure The degree of freedom matrix is ,calculate , ; If plane strain elements are used, the generalized adjoint force for each element is calculated as follows: ; If axisymmetric elements are used, the generalized adjoint forces for each element are calculated using the following formula: ;in, Representative unit Degrees of freedom of midpoint L The corresponding vector components; is the shape function of node L; , They are the reference configuration coordinates and the current configuration coordinates respectively; is the material tangent stiffness tensor; is the deformation gradient tensor; Each calculated unit The generalized adjoint force vector After assembly, the generalized adjoint force vector of the entire structure is obtained ,calculate , substitute into the adjoint equation and solve: ; Find the first Lagrange multiplier After that, calculate , .

[0012] Calculated by central difference method , , ; The partial derivative of the eigenvalue with respect to the physical density is calculated by the following formula: ; Substitute the calculation results into the following formula to calculate the partial derivative of the KS aggregation function with respect to the physical density: ; Calculation of critical buckling load About the partial derivative of physical density: ; (4) Obtain the current loading displacement as the loading displacement corresponding to the first critical point , solve The adjoint equation corresponding to the sensitivity calculation is: ; In the formula, The corresponding component of the "control degree of freedom" is 1, and the other components are 0. The dimension is the same as Same vector; Find the second Lagrange multiplier After that, calculate , ,calculate About the partial derivative of physical density: .

[0013] Optionally, in step 5, the calculation formula of the sensitivity corresponding to the optimization objective function is as follows: ; In the formula, is the sensitivity of the optimization objective function, is the critical buckling load, is the first critical point corresponding to the loading displacement, is the loading displacement when the load is the critical buckling load in the post-buckling stage, is the design variable, is the filter density, is the physical density.

[0014] The stiffness, buckling load and theoretical maximum specific energy absorption of the multi-stable structure formed by the elastic buckling multi-stable energy absorption structure optimized by the topology optimization method of the elastic buckling multi-stable energy absorption structure of the present application are greatly improved, indicating the effectiveness of the optimization method proposed in the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a flow chart of a topology optimization method of an elastic buckling multi-stable energy absorbing structure according to an embodiment of the present application.

[0016] Figure 2 Schematic diagram of the process of topology optimization method of elastic buckling multi-stable energy absorbing structure according to an embodiment of the present application.

[0017] Figure 3 This is a coupled degree of freedom loading method according to an embodiment of the present application.

[0018] Figure 4 This is a flow chart of the pre-critical point identification algorithm of an embodiment of the present application.

[0019] Figure 5A bistable structure optimized for an embodiment of the present application.

[0020] Figure 6 A multi-stable structure formed by a bistable cell array according to an embodiment of the present application.

[0021] Figure 7 This is a bistable structure design domain of an embodiment of the present application.

[0022] Figure 8 This is a non-optimized cosine-shaped shell structure used for comparison in the embodiments of the present application.

[0023] Fig. 9 This is a bistable curved shell structure obtained by optimization according to an embodiment of the present application.

[0024] Fig.10 This is the cross-sectional shape of the bistable curved shell structure obtained by optimization according to the embodiment of the present application.

[0025] Reference numerals DETAILED DESCRIPTION

[0026] In order to make the purpose, technical scheme and advantages of the implementation of this application clearer, the technical scheme in the embodiment of this application will be described in more detail below in conjunction with the drawings in the embodiment of this application. In the drawings, the same or similar reference numerals throughout represent the same or similar elements or elements with the same or similar functions. The described embodiments are part of the embodiments of this application, not all of them. The embodiments described below with reference to the drawings are exemplary and are intended to be used to explain this application, and should not be construed as limitations on this application. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of this application. The embodiments of this application are described in detail below in conjunction with the drawings.

[0027] In the description of the present application, it should be understood that the terms "center", "longitudinal", "lateral", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", etc., indicating orientations or positional relationships, are based on the orientations or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present application and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be understood as limiting the scope of protection of the present application.

[0028] See also Figure 5 as well as Figure 6 , Figure 5 The elastic buckling multi-stable energy absorbing structure of the present application (such as Figure 6), the present application performs a bistable structure topology optimization design based on the traditional cosine-shaped two-dimensional curved shell, and the cosine-shaped two-dimensional curved shell includes an elastic part 1 and a rigid part 2 (as shown in FIG. Figure 5 As shown), the purpose of this application is to optimize the material distribution characteristics of the elastic part 1 of each basic unit in the elastic buckling multi-stable energy absorption structure of this application.

[0029] like Figure 1 The topology optimization method for the elastic buckling multi-stable energy-absorbing structure shown includes: Step 1: Initialize the topology optimization finite element model: given the design domain and material properties, discretize the design domain into a finite element mesh, initialize the design variables, couple the degrees of freedom of the loading area, and apply corresponding boundary conditions to the finite element model; Step 2: Update the filter density, physical density, and material properties of each unit in the finite element model established in step 1 according to the current design variables; Step 3: Perform finite element calculation of the pre-buckling stage on the finite element model updated in step 2 according to the internal force vector and each stiffness matrix, and calculate the critical buckling load, critical buckling load sensitivity, the first critical point corresponding loading displacement and the first critical point corresponding sensitivity under the pre-critical point state; Step 4: After obtaining the critical buckling load, the critical buckling load sensitivity, the loading displacement corresponding to the first critical point, and the sensitivity corresponding to the first critical point, perform finite element calculation in the post-buckling stage until the load reaches the critical buckling load, and calculate the loading displacement and loading displacement sensitivity when the critical buckling load is reached; Step 5: Calculate the optimization objective function, the constraint function, the sensitivity corresponding to the optimization objective function, and the sensitivity corresponding to the constraint function according to the information obtained in step 3 and step 4; Step 6: Use the moving asymptote method to construct the optimization subproblem and solve it to obtain the updated design variables; Step 7: Determine whether the topology optimization has converged. If so, the optimization ends. Otherwise, the updated design variables are used as current design variables and steps 2 to 6 are repeated until the topology optimization is determined to have converged.

[0030] The stiffness, buckling load and theoretical maximum specific energy absorption of the multi-stable structure formed by the elastic buckling multi-stable energy absorption structure optimized by the topology optimization method of the elastic buckling multi-stable energy absorption structure of the present application are greatly improved, indicating the effectiveness of the optimization method proposed in the present invention.

[0031] When performing the finite element calculation, the optimization objective function and the constraint function adopt the following optimization formula: ;in, is the theoretical maximum energy absorption per bistable unit; is the design variable; is the finite element structure displacement vector; is the critical buckling load, is the first critical point corresponding to the loading displacement, is the loading displacement when the load is the critical buckling load in the post-buckling stage; It is an unbalanced force; , The units are Design variables and material volumes; is the total volume of the design domain, is the maximum specified material volume fraction.

[0032] In finite element calculation, unbalanced force The expression is: ;in, is the node internal force vector when the degree of freedom coupling and displacement boundary conditions are not considered, is the node external force vector when the degree of freedom coupling and displacement boundary conditions are not considered; matrix Defined as: ;in, is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; yes The vector of degrees of freedom that are not constrained by displacement boundary conditions; matrix Defined as: ;in, is the finite element structure displacement vector; is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; The internal force vector considering the degree of freedom coupling and displacement boundary conditions Use the following formula to obtain: ;in, The stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions includes a shear stiffness matrix, a linear stiffness matrix, and a stress stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions. as follows: ; The stress stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the stress stiffness matrix without considering degree of freedom coupling and displacement boundary conditions; The linear stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the linear stiffness matrix without considering degree of freedom coupling and displacement boundary conditions.

[0033] In this embodiment, the filtering density of each unit in step 2 is obtained by the following method: Through the filtering method, according to the design variables Calculate the filter density ; Through the smooth Heaviside projection function, according to the filtered density Calculate the physical density used for interpolation of element material properties .

[0034] In this embodiment, the node internal force vector is obtained by assembling the unit internal force vector of each unit. The unit internal force vector Use the following formula to obtain: ;in, is the displacement used for the nonlinear part of the strain energy density calculation, is the element internal force vector calculated using the large deformation theory, is the element internal force vector calculated under the small deformation assumption; The tangent stiffness matrix without considering the coupling of degrees of freedom and displacement boundary conditions is assembled from the element stiffness matrix of each element. The element stiffness matrix Use the following formula to obtain: ; Among them, the element stiffness matrix The tangent stiffness matrix without considering the coupling of degrees of freedom and displacement boundary conditions is assembled from the element stiffness matrix of each element; is the element stiffness matrix calculated using the large deformation theory, is the element stiffness matrix calculated under the small deformation assumption.

[0035] In this embodiment, the finite element calculation of the pre-buckling stage is performed on the finite element model updated in step 2 according to the internal force vector and each stiffness matrix, and the critical buckling load, the critical buckling load sensitivity, the loading displacement corresponding to the first critical point and the sensitivity corresponding to the first critical point are calculated in the pre-critical point state, including: In step 3, the equilibrium state of the previous incremental step needs to be saved during finite element calculation, and the previous critical point is continuously identified. When it is identified that the current equilibrium state has passed the first critical point, the current equilibrium state is saved, the previous equilibrium state is recorded as the previous critical point, and the finite element structure is restored to the state of the previous critical point. Under the state of the previous critical point, the critical buckling load, critical buckling load sensitivity, loading displacement corresponding to the first critical point, and sensitivity corresponding to the first critical point are calculated through the following steps: The present application is further described in detail below by way of examples. It should be understood that the examples do not constitute any limitation to the present application.

[0036] Step 1: Determine the material properties. In this embodiment, the material constitutive model is the Neo-Hookean hyperelastic model, and the physical property parameter is the elastic modulus. , Poisson's ratio ; Determine the design domain as Figure 7 As shown, , , , The outer dimension of the cross section is The design domain adopts The four-node fully integrated plane strain element is discretized, and the design variable is initialized to 0.5, and the maximum material volume fraction is set to 0.683. A fixed support boundary condition is applied to the left end of the design domain, and a symmetric boundary condition is applied to the right end. After coupling the degrees of freedom of the top "platform" area, a force or displacement load is applied downward (such as Figure 3 ), in order to simulate the compression process of a single bistable structure.

[0037] Step 2: According to the current design variables (in this embodiment, when the 0th round of iteration is performed, the current design variables are the initialization design variables, and the design variables of each subsequent round are the updated design variables obtained in the previous round), the filtering density of each unit is calculated and updated by the filtering method, and the physical density of each unit is calculated and updated by the projection method.

[0038] Step 3: Perform finite element calculations in the pre-buckling stage on the finite element model updated in step 2: switch the nonlinear equation solving algorithm to the arc length method, switch the loading method to coupled degree of freedom displacement control loading, and set the loading displacement value to 32 mm. Perform finite element iterative calculations, and use Figure 4The pre-critical point identification algorithm shown checks whether the current incremental step has passed the pre-critical point. If not, the next incremental step is calculated. Otherwise, the current incremental step is used as the pre-critical point.

[0039] In this embodiment, the finite element structure is restored to the state of the critical point before the increase, the loading mode is switched to force loading, and the buckling load is calculated. and its partial derivative with respect to the physical density , the loading displacement corresponding to the first critical point and its partial derivative with respect to the physical density , which can be calculated by following the steps below: (1) The one-point method is used as the calculation method for nonlinear buckling load, and the first critical point is , the generalized eigenvalue problem is solved by the subspace iteration method of Rayleigh-Ritz projection: (1) In the formula is the stress stiffness matrix considering degree of freedom coupling and displacement boundary conditions, is the linear stiffness matrix considering degree of freedom coupling and displacement boundary conditions; is the modal eigenvalue ( ), is the corresponding eigenmode, subscript Represents the order of the buckling mode.

[0040] (2) After solving the generalized eigenvalue problem, calculate the critical buckling load: (2) In the formula, is the external force applied on the controlled degree of freedom in the pre-critical state, is the order of the buckling mode chosen to approximate the maximum eigenvalue, It is a KS aggregation function, and its expression is: (3) (3) Select the unit from the global degrees of freedom of the structure The degree of freedom matrix is ,calculate , . Compute the generalized adjoint force for each element: (4) In formula (4) Representative unit I Middle Node L Degrees of freedom The corresponding vector components; It is a node L The shape function of , They are the reference configuration coordinates and the current configuration coordinates respectively; is the material tangent stiffness tensor; is the deformation gradient tensor.

[0041] Each calculated unit The generalized adjoint force vector After assembly, the generalized adjoint force vector of the entire structure is obtained ,calculate , substitute into the adjoint equation and solve: (5) Find the first Lagrange multiplier After that, calculate , .

[0042] Calculated by central difference method , , .

[0043] The partial derivative of the eigenvalue with respect to the physical density is calculated by formula (6), and the result is substituted into formula (7) to calculate the partial derivative of the KS aggregation function with respect to the physical density.

[0044] (6) (7) Calculation of buckling loads About the partial derivative of physical density: (8) (4) Obtain the current loading displacement as the loading displacement corresponding to the first critical point , solve The adjoint equation corresponding to the sensitivity calculation is: (9) In the formula, is the dimension in which the corresponding component of the "control degree of freedom" is 1 and the other components are 0. Same vector.

[0045] Find the second Lagrange multiplier After that, calculate , .calculate About the partial derivative of physical density: (10) Step 4: After obtaining the critical buckling load, critical buckling load sensitivity, loading displacement corresponding to the first critical point, and sensitivity corresponding to the first critical point, perform finite element calculations in the post-buckling stage until the load reaches the critical buckling load, and calculate the loading displacement and loading displacement sensitivity when the critical buckling load is reached. Specifically, the following method is used for the subsequent post-buckling stage finite element calculations: Restore the finite element state to the state when the increment step converges after the previous critical point, switch the loading mode to displacement loading, and continue to perform finite element iterative calculations until the displacement of the "controlled degree of freedom" reaches the set loading displacement. Switch the nonlinear equation solving algorithm to Newton's method, switch the loading mode to force-controlled loading, and set the loading force to the previous critical load. Continue to perform finite element iterative calculations until the external force on the "controlled degree of freedom" reaches the previous critical load.

[0046] Read the displacement of the "controlled degree of freedom" at this time as In the finite element state of the current increment, solve The adjoint equation (9) corresponding to the sensitivity calculation is used to obtain the first Lagrange multiplier After that, calculate , .calculate About the partial derivative of physical density: Step 5: Based on the calculation results of steps 3 and 4, calculate the constraint function (structural material volume fraction) and sensitivity, and calculate the objective function (i.e., the theoretical maximum value of energy absorption per bistable unit): (12) Calculate the sensitivity of the objective function; (13) Step 6: Use the moving asymptote method to construct the optimization subproblem and solve it to obtain the updated design variables; Step 7: Determine whether the topology optimization has converged. If so, the optimization ends; otherwise, return to step 2.

[0047] Figure 8 is a non-optimized cosine-shaped shell structure used for comparison, where , and the remaining geometric dimensions are the same as the design domain of the embodiment. The stiffness of the unoptimized cosine-shaped curved shell structure is 80.87 N / mm, the buckling load is 128.37 N, and the theoretical maximum specific energy absorption is 179.58 J / kg.

[0048] The bistable curved shell structure obtained in this topology optimization example is as follows Fig. 9 Its cross-sectional shape is as shown Fig.10As shown in Figure 1, the total material volume is equal to that of the unoptimized cosine-shaped shell structure. The stiffness of the optimized structure is 108.31 N / mm, the buckling load is 331.26 N, and the theoretical maximum specific energy absorption is 444.98 J / kg.

[0049] It can be seen that the stiffness, buckling load and theoretical maximum specific energy absorption of the multi-stable structure formed by the bistable curved shell structure array obtained by topology optimization are greatly improved, indicating the effectiveness of the optimization method proposed in the present invention.

[0050] Although the present application is disclosed as above in terms of a preferred embodiment, it is not intended to limit the present application. Any technical personnel in this field may make possible changes and modifications without departing from the spirit and scope of the present application. Therefore, the scope of protection of the present application shall be based on the scope defined by the claims of the present application.

[0051] Finally, it should be pointed out that the above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application is described in detail with reference to the above embodiments, a person skilled in the art should understand that the technical solutions described in the above embodiments can still be modified, or some of the technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A topology optimization method for an elastic buckling multi-stable energy-absorbing structure, characterized in that: The topology optimization method of the elastic buckling multi-stable energy absorption structure comprises: Step 1: Initialize the topology optimization finite element model: given the design domain and material properties, discretize the design domain into a finite element mesh, initialize the design variables, couple the degrees of freedom of the loading area, and apply corresponding boundary conditions to the finite element model; Step 2: Update the filter density, physical density, and material properties of each unit in the finite element model established in step 1 according to the current design variables; Step 3: Perform finite element calculation of the pre-buckling stage on the finite element model updated in step 2 according to the internal force vector and each stiffness matrix, and calculate the critical buckling load, critical buckling load sensitivity, the first critical point corresponding loading displacement and the first critical point corresponding sensitivity under the pre-critical point state; Step 4: After obtaining the critical buckling load, the critical buckling load sensitivity, the loading displacement corresponding to the first critical point, and the sensitivity corresponding to the first critical point, perform finite element calculation in the post-buckling stage until the load reaches the critical buckling load, and calculate the loading displacement and loading displacement sensitivity when the critical buckling load is reached; Step 5: Calculate the optimization objective function, the constraint function, the sensitivity corresponding to the optimization objective function, and the sensitivity corresponding to the constraint function according to the information obtained in step 3 and step 4; Step 6: Use the moving asymptote method to construct the optimization subproblem and solve it to obtain the updated design variables; Step 7: Determine whether the topology optimization has converged. If so, the optimization ends. Otherwise, the updated design variables are used as current design variables and steps 2 to 6 are repeated until the topology optimization is determined to have converged.

2. The topology optimization method of the elastic buckling multi-stable energy absorption structure according to claim 1 is characterized in that: When performing the finite element calculation, the optimization objective function and constraint function adopt the following optimization formula: ;in, is the theoretical maximum energy absorption per bistable unit; is the design variable; is the finite element structure displacement vector; is the critical buckling load, is the first critical point corresponding to the loading displacement, is the loading displacement when the load is the critical buckling load in the post-buckling stage; It is an unbalanced force; , The units are Design variables and material volumes; is the total volume of the design domain, is the maximum specified material volume fraction.

3. The topology optimization method of the elastic buckling multi-stable energy absorption structure according to claim 2 is characterized in that: In finite element calculation, unbalanced force The expression is: ;in, is the node internal force vector when the degree of freedom coupling and displacement boundary conditions are not considered, is the node external force vector when the degree of freedom coupling and displacement boundary conditions are not considered; matrix Defined as: ;in, is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; yes The vector of degrees of freedom that are not constrained by displacement boundary conditions; matrix Defined as: ;in, is the finite element structure displacement vector; is the displacement vector The vector obtained by removing the coupled degrees of freedom from ; The internal force vector considering the degree of freedom coupling and displacement boundary conditions Use the following formula to obtain: ;in, The stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions includes a shear stiffness matrix, a linear stiffness matrix, and a stress stiffness matrix considering the coupling of degrees of freedom and displacement boundary conditions. as follows: ； The stress stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the stress stiffness matrix without considering degree of freedom coupling and displacement boundary conditions; The linear stiffness matrix considering degree of freedom coupling and displacement boundary conditions is as follows: ;in, is the linear stiffness matrix without considering degree of freedom coupling and displacement boundary conditions.

4. The topology optimization method of the elastic buckling multi-stable energy absorption structure according to claim 3 is characterized in that: The filtering density of each unit in step 2 is obtained by the following method: Through the filtering method, according to the design variables Calculate the filter density ; Through the smooth Heaviside projection function, according to the filtered density Calculate the physical density used for interpolation of element material properties .

5. The topology optimization method of the elastic buckling multi-stable energy absorption structure according to claim 4, characterized in that: The node internal force vector is assembled from the unit internal force vectors of each unit. The unit internal force vector Use the following formula to obtain: ;in, is the displacement used for the nonlinear part of the strain energy density calculation, is the element internal force vector calculated using the large deformation theory, is the element internal force vector calculated under the small deformation assumption; The tangent stiffness matrix without considering the coupling of degrees of freedom and displacement boundary conditions is assembled from the element stiffness matrix of each element. The element stiffness matrix Use the following formula to obtain: ;in, is the element stiffness matrix calculated using the large deformation theory, is the element stiffness matrix calculated under the small deformation assumption.

6. The topology optimization method of the elastic buckling multi-stable energy absorption structure according to claim 5, characterized in that: In step 5, the calculation formula of the sensitivity corresponding to the optimization objective function is as follows: ； In the formula, is the sensitivity of the optimization objective function, is the critical buckling load, is the first critical point corresponding to the loading displacement, is the loading displacement when the load is the critical buckling load in the post-buckling stage, is the design variable, is the filter density, is the physical density.

Images