BESO-based viscoelasticity isogeometric topology optimization method under transient load

Abstract

The invention provides a viscoelasticity isogeometric topology optimization method under transient load based on BESO, and belongs to the technical field of isogeometric topology optimization, and the method comprises the following steps: carrying out isogeometric gridding on a model to be designed, and obtaining a control point and a Gaussian point coordinate corresponding to each unit; constructing a long-term global stiffness matrix and a global quality matrix; determining an optimized target volume; on the basis of constructing a long-term global stiffness matrix, a tangent stiffness matrix and a historical stiffness matrix are obtained, and geometric analysis such as viscoelasticity is carried out. Calculating a target value of the structure according to the optimization problem definition, and calculating the sensitivity of a target function relative to a design variable by adopting an adjoint method; after the sensitivity value of each unit is filtered through a distance filter, the sensitivity value is expanded into sensitivity at a control point through an NURBS filter, and a global sensitivity field is established; and updating the design variables according to the set evolution ratio parameters. According to the invention, the response of the structure under complex dynamic loading can be accurately simulated.

Description

Technical Field

[0001] The present invention relates to the field of geometric analysis technology, and in particular to a geometric topology optimization method such as viscoelasticity under transient load based on BESO. Background Art

[0002] With the rapid development of computer technology, topology optimization has been widely used as a powerful tool for achieving innovative design. Among them, the bidirectional evolutionary structural optimization (BESO) method is an algorithm for optimizing the material distribution in the structure. Its goal is to find the optimal material distribution that meets specific performance and constraints in a given space, while minimizing the use of materials through the selection of binary materials.

[0003] Different from traditional finite element analysis (FEA), isogeometric analysis (IGA) uses non-uniform rational B-spline (NURBS) as a geometric description tool, directly combining geometric modeling with numerical analysis, avoiding the errors caused by geometric transformation in traditional finite element method, thus realizing the seamless connection between computer-aided design (CAD) and computer-aided engineering (CAE). The application of IGA has significantly improved the efficiency and accuracy of numerical simulation and design processes, so the isogeometric topology optimization method based on BESO has received widespread attention. However, most of the current topology optimization methods based on IGA rely on linear elastic assumptions and are mainly used for static analysis, ignoring the common nonlinear characteristics of materials and complex load conditions in engineering, especially the limitations of viscoelastic materials. Viscoelastic materials exhibit time-dependent behavior when subjected to loads, including both transient elastic response and time-varying viscous flow. Under transient loads, viscoelastic materials exhibit stress relaxation or creep phenomena. This time-dependent nonlinear mechanical behavior is particularly critical in many engineering applications.

[0004] Existing isogeometric topology optimization methods usually fail to fully consider these time-dependent characteristics of viscoelastic materials, resulting in certain limitations in engineering applications involving structural stress-strain evolution under transient loads. Therefore, in order to better solve this problem, a viscoelastic isogeometric topology optimization method under transient loads based on BESO is proposed. It is very necessary to analyze the time-dependent response of the structure through the viscoelastic constitutive relationship, which can more accurately simulate the viscoelastic behavior of the structure and optimize the material distribution to meet the complex mechanical requirements in actual engineering. Summary of the invention

[0005] In view of this, the present invention proposes a geometric topology optimization method such as viscoelasticity under transient loads based on BESO, which can accurately simulate the response of the structure under complex dynamic loading and can also consider the viscoelasticity and time-varying characteristics of the material during the optimization process.

[0006] The present invention provides a geometric topology optimization method based on BESO based on viscoelasticity under transient load, comprising the following steps:

[0007] S1: According to the initially set grid size, the model to be designed iso-geometrically meshed to obtain an iso-geometric grid model, the iso-geometric grid model is divided, and the coordinates and numbers of the control points and Gaussian points corresponding to each unit are obtained; the convergence condition and the evolution ratio parameter are determined;

[0008] S2: Assemble the unit stiffness matrix and unit mass matrix using the coordinates and numbers of the control points and Gaussian points, summarize the stiffness matrix and mass matrix of each unit according to the degree of freedom number, and construct the long-term global stiffness matrix and global mass matrix under the isogeometric mesh model;

[0009] S3: Determine the optimized target volume according to the evolution ratio parameter and the current design domain volume;

[0010] S4: Based on the construction of the long-term global stiffness matrix, the tangent stiffness matrix and the historical stiffness matrix are obtained to perform geometric analysis such as viscoelasticity. The structural response at each time step is obtained by solving the equilibrium equation of the viscoelastic material.

[0011] S5; according to the optimization problem definition, the target value of the structure is calculated, and the sensitivity of the target function with respect to the design variables is calculated by the adjoint method;

[0012] S6: After filtering the sensitivity value of each unit through the distance filter, it is expanded into the sensitivity at the control point through the NURBS filter to establish a global sensitivity field;

[0013] S7: Update the design variables according to the set evolution ratio parameters;

[0014] S8: Check whether the updated design variables meet the convergence conditions. If they do, output the optimization results; if not, return to step S3.

[0015] On the basis of the above technical solution, preferably, the content of step S1 is to use the H-refinement method to refine the design domain according to the preset initial grid size, generate an isogeometric grid model of a specified grid size, divide the isogeometric grid model, and determine the coordinates and numbers of the corresponding control points and Gaussian points in each unit according to the characteristics of the isogeometric grid model, and determine the convergence condition, filter radius and evolution ratio parameter according to the requirements, wherein the relationship of the convergence condition is: τ and N0 represent the allowable convergence error and a positive integer respectively. D Represents the target value under a certain iteration, i0 is the current iteration number, J∈N0.

[0016] Preferably, the content of step S2 is that the long-term element stiffness matrix k ∞ And the global mass matrix M is solved according to the following relationship: Where N is the NURBS basis function, B is the strain-displacement matrix, N1, N2, ...N n represents the NURBS basis function at different control points, (x, y) is the physical space coordinate, Ω is the physical domain, is the integration domain, D ∞ is the long-term elasticity matrix, v is the Poisson's ratio of the material, E ∞ represents the long-term elastic modulus, J1 and J2 represent the transformation relationship from NURBS parameter space to physical space and from integral parameter space to NURBS parameter space, respectively. (ξ,η) are the coordinates in parameter space, is the integral space coordinate, ξ i0+1 and i0 is the coordinate of the adjacent control point along the ξ direction in the parameter coordinate system, η i0+1 and η i0 are the coordinates of the adjacent control points along the η direction in the parameter coordinate system.

[0017] Further preferably, the content of step S3 is to set the volume fraction V of the j0th optimization iteration j0 The calculation formula is: V j0 =max[V req , V j0-1 (1-er)], er represents the evolution ratio, V req represents the target volume, V j0-1 represents the volume fraction of the j0-1th optimization iteration.

[0018] More preferably, the specific solution steps of the viscoelastic and other geometric analysis in step S4 are as follows: first, the tangent stiffness matrix k is obtained respectively through the long-term unit stiffness matrixT and the historical stiffness matrix k hist : represents the number of terms in the Prony series, γ jp =E jp / E ∞ , Δt is the interval of each time step, τ jp and E jp are the relaxation time and relaxation modulus respectively; for linear viscoelastic materials, the total stress σ n+1 Defined as C ∞ According to the long-term elastic modulus E ∞ The constitutive matrix solved, time step iT = 1, 2, ..., n, n + 1, ..., N t , They represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, respectively. n+1 ,t n-1 ,t iT It represents the time required to reach the n+1th, n-1th, and iTth time steps. According to the definition of internal force, the internal force is solved by the following formula: Rewrite the above formula as The variables The internal force at time step n+1 is expressed as Each time step ensures the balance of forces, so the internal force is equal to the external force, satisfying is the external force, and the external force at time step n+1 is expressed as is the historical force vector at time step n; is the historical parameter at time step n, indicating the influence of n-1 time steps before the current time step on time step n;

[0019] Assemble the above-mentioned unit stiffness matrix, historical stiffness matrix, and historical parameters into a global form, and let the global stiffness matrix Global history stiffness matrix Global History Vector and The external force vector represents the global form of the external force at time step iT and n+1, where e∈N e , N e is the number of units, and the subscript e is the unit number; the equilibrium equation of the system is expressed as: According to the smoothing equation, the overall displacement vector at each time step is obtained u n and u n+1 represents the overall displacement vector corresponding to time steps n and n+1, are the global forms corresponding to the historical force vector and historical parameters at time step n respectively.

[0020] Further preferably, the content of step S5 is that, according to the definition of the objective function, the target value E is solved by the following formula: where u iT and u iT-1 is the overall displacement vector corresponding to time steps iT and iT-1 obtained according to the overall displacement vector formula of step S4, is the external force vector in global form corresponding to the external force at time step iT-1; the residual R at each time step iT It is expressed by the following formula: is the global form corresponding to the historical parameters at time step iT-1; introduce the adjoint vector multiplier ψ iT , the objective function is expressed as In order to obtain the sensitivity of the design variable ρ, the objective function is differentiated with respect to the design variable, and the sensitivity expression is obtained as follows: in are the external force vectors in global form corresponding to the external forces at the last two time steps; the adjoint matrix multiplier of the last residual equation is

[0021] Solve 1 to N in descending order t The adjoint matrix multiplier of -1 time step is expressed as follows: where δ k-1，iT is the Cronel function, when k-1 = iT, the function value is 1, otherwise it is 0; the residual term R iT The partial derivative with respect to displacement is Will get and the adjoint vector multiplier ψ iT Substituting into the sensitivity expression, we can obtain the sensitivity of the objective function relative to the design variables.

[0022] More preferably, step S6 includes the following contents:

[0023] S61: filtering and averaging the sensitivity values ​​so that the updated sensitivity values ​​contain all sensitivity information of previous iterations, and obtaining the sensitivity field of the global control point;

[0024] S62: Utilizing the threshold parameters of material addition and removal obtained through the iterative algorithm and implementing the updating of topological variables.

[0025] Further preferably, the content of step S61 is to filter the unit sensitivity values ​​according to the following relationship: Where nel is based on unit e as the filter center, r min is all cells in the neighborhood of radius, rmin is the sensitivity filter radius, the linear weight ω(r ej )=max[0,r min -dis(e,j)], where dis(e,j) is the distance between the centers of cells e and j; α e and α j are the sensitivity values ​​of units e and j respectively; the sensitivity is filtered using a NURBS filter to obtain the sensitivity value based on the control point. The specific relationship is as follows: where ρ i , ei and eij represent the variable of the ith control point, the unit set affected by the ith control point and the jth unit in ei, respectively. ij is a NURBS basis function; the control point sensitivity values ​​are averaged so that the updated sensitivity values ​​contain all sensitivity information of the previous iterations: Where m represents the number of optimizations, thereby establishing a global control point sensitivity field.

[0026] More preferably, the specific steps of the iterative algorithm in step S62 are as follows:

[0027] 1) Threshold parameter for adding materials and the threshold parameter for material removal Satisfies the following relationship: α th Optimize the threshold value for the design variable according to the volume fraction V of the i0+1th iteration i0+1 Sure;

[0028] 2) Calculate the allowable volume ratio AR, which is defined as the number of added units divided by the total number of units in the current design; if AR < er, skip step 3); if AR ≥ er, perform step 3) and recalculate and

[0029] 3) Calculate by first sorting the sensitivity values ​​of empty cells 0 The number of cells that switch from 0 to 1 will be equal to the evolution ratio er multiplied by the total number of cells in the current design, is the sensitivity of the element immediately following the last added unit; then determine Make the volume removed equal to (V i0+1 -V i0 + volume of the added unit), V i0 is the volume fraction of the i0th iteration;

[0030] The updating of topological variables is performed according to the following scheme:

[0031] Further preferably, the value of the allowable convergence error τ is 0.001; the value of the positive integer N0 is 5;

[0032] The present invention provides a geometric topology optimization method based on BESO under transient loads, which has the following beneficial effects compared with the prior art:

[0033] (1) The present invention uses viscoelastic isogeometric analysis technology under transient loads: it adopts advanced viscoelastic modeling technology, takes into account the time dependence and hysteresis effect of materials under transient loads, and can capture the complex deformation behavior of materials during loading and unloading; by introducing isogeometric analysis technology, NURBS basis functions are directly used to represent geometric shapes and responses, avoiding geometric conversion errors in traditional finite element methods, and more accurately simulating the viscoelastic response of structures under transient loads. It can not only accurately simulate the response of structures under complex dynamic loading, but also consider the viscoelasticity and time-varying characteristics of materials during the optimization process. This has broad application prospects for engineering applications involving long-term loads or periodic loads, such as bridges, building structures, and mechanical equipment.

[0034] (2) Different from the traditional topology optimization method based on linear elastic or elastoplastic analysis, the present invention comprehensively considers the viscoelastic behavior of the material and adjusts the topological layout of the structure in combination with the influence of transient loads. This can more accurately reflect the time-varying characteristics and energy dissipation of the material during the optimization process, ensuring that the final structure has the best performance under complex load conditions, especially in shock absorption, energy dissipation and long-term loading environments.

[0035] (3) A multi-layer filter is used to smooth the sensitivity information, which not only avoids the common grid dependency and checkerboard phenomenon, but also ensures that each iteration fully considers the previous optimization information through the average processing of historical data. Especially for the topology optimization of viscoelastic materials, due to the complexity of transient response, this processing method effectively improves the stability and convergence speed of sensitivity calculation, and ensures the reliability of the final optimization results. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0037] Figure 1 It is a flow chart of a geometric topology optimization method based on BESO under transient loads;

[0038] Figure 2 A schematic diagram of a Kelvin-Voigt constitutive model of a viscoelastic material based on a BESO-based viscoelastic isogeometric topology optimization method under transient loads according to the present invention;

[0039] Figure 3 A schematic diagram of control points of units in an isogeometric model mesh topology optimization structure of a viscoelastic isogeometric topology optimization method under transient load based on BESO of the present invention;

[0040] Figure 4 It is a schematic diagram of density filtering and NURBS filtering constructed in an embodiment of a geometric topology optimization method based on BESO under transient loads;

[0041] Figure 5 It is a schematic diagram of an optimization example obtained by a BESO basic method in a BESO-based viscoelastic and other geometric topology optimization method under transient loads of the present invention;

[0042] Figure 6 It is a schematic diagram of optimization examples under different simulation times obtained in an embodiment of a geometric topology optimization method based on viscoelasticity under transient loads of BESO in the present invention. DETAILED DESCRIPTION

[0043] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0044] Existing isogeometric topology optimization methods usually fail to fully consider these time-dependent properties of viscoelastic materials, resulting in certain limitations in engineering applications involving structural stress-strain evolution under transient loads. Figure 1 As shown, the present invention provides a geometric topology optimization method based on BESO based on viscoelasticity under transient load, comprising the following steps:

[0045] S1: According to the initially set grid size, the model to be designed iso-geometrically meshed to obtain an iso-geometric grid model, the iso-geometric grid model is divided, and the coordinates and numbers of the control points and Gaussian points corresponding to each unit are obtained; the convergence condition and the evolution ratio parameter are determined;

[0046] This method adopts the H-refinement method, which generates an isogeometric mesh model of a specified mesh size by refining the design domain according to the preset initial mesh size, divides the isogeometric mesh model, and determines the coordinates and numbers of the corresponding control points and Gauss points in each unit according to the characteristics of the isogeometric mesh model. Figure 3 As shown. Determine the convergence condition, filter radius, evolution ratio and other parameters according to actual needs, where the convergence condition relationship is: τ and N0 represent the allowable convergence error and a positive integer respectively. D represents the target value of a certain iteration, i0 is the current iteration number, and J∈N0. In this embodiment, the value of the allowed convergence error τ is 0.001; the value of the positive integer N0 is 5.

[0047] like Figure 2 Schematic diagram of the Kelvin-Voigt constitutive model of the selected viscoelastic material of the present invention, including an elastic part and a damping part.

[0048] S2: Use the coordinates and numbers of the control points and Gaussian points to assemble the unit stiffness matrix and unit mass matrix, summarize the stiffness matrix and mass matrix of each unit according to the degree of freedom number, and construct the long-term global stiffness matrix and global mass matrix under the isogeometric grid model.

[0049] Long-term element stiffness matrix k ∞ And the global mass matrix M is solved according to the following relationship: Where N is the NURBS basis function, B is the strain-displacement matrix, N1, N2, ...N n represents the NURBS basis function at different control points, (x, y) is the physical space coordinate, Ω is the physical domain, is the integration domain, D ∞ is the long-term elasticity matrix, v is the Poisson's ratio of the material, E ∞ represents the long-term elastic modulus, J1 and J2 represent the transformation relationship from NURBS parameter space to physical space and from integral parameter space to NURBS parameter space, respectively. (ξ,η) are the coordinates in parameter space, is the integral space coordinate, ξ i0+1 and i0 is the coordinate of the adjacent control point along the ξ direction in the parameter coordinate system, η i0+1 and η i0 are the coordinates of the adjacent control points along the η direction in the parameter coordinate system.

[0050] S3: Determine the optimized target volume according to the evolution ratio parameter and the current design domain volume.

[0051] The target volume of the local iteration is determined by the evolution ratio and the current design domain volume.

[0052] Let the volume fraction V of the j0th optimization iteration be j0 The calculation formula is: V j0 =max[V req , V j0 -1(e=1-er)], er represents the evolution ratio, V req represents the target volume, V j0-1 It is important to note that the actual volume obtained during the iterations usually deviates slightly from the specified target value, but once the desired material volume fraction is reached, the optimization algorithm will only change the topology and the volume fraction will always remain the same.

[0053] S4: Based on the construction of the long-term global stiffness matrix, the tangent stiffness matrix and the historical stiffness matrix are obtained to perform geometric analysis such as viscoelasticity; the structural response at each time step is obtained by solving the equilibrium equations of the viscoelastic material.

[0054] The specific solution steps for geometric analysis such as viscoelasticity are as follows: First, the tangent stiffness matrix k is obtained through the long-term element stiffness matrix T and the historical stiffness matrix k hist : represents the number of terms in the Prony series, γ jp =E jp / E ∞ , Δt is the interval of each time step, τ jp and E jp They are relaxation time and relaxation modulus respectively, which can be found in the table below.

[0055] Prony series number of terms jp <![CDATA[Relaxation time τ jp Unit: second]]> <![CDATA[Relaxation modulus E jp in GPa]]> 1 0.1 0.5 2 1 0.3 3 10 0.13 4 100 0.07 ∞ ∞ 0.30 0 - 1.30

[0056] The constitutive matrix D of the two-dimensional plane state jp It can be expressed in terms of relaxation modulus and Poisson's ratio:

[0057] For a linear viscoelastic material, the total stress σ n+1 Defined as C ∞ According to the long-term elastic modulus E ∞ The constitutive matrix solved, time step iT = 1, 2, ..., n, n+1, ..., N t , They represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, respectively. n+1 ,t n-1 ,t iT It represents the time required to reach the n+1th, n-1th, and iTth time steps. According to the definition of internal force, the internal force is solved by the following formula: Rewrite the above formula as The variables The internal force at time step n+1 is expressed as

[0058] During the analysis, each time step must ensure the balance of forces, so the internal force is equal to the external force, satisfying is the external force, and further deduction shows that the external force at time step n+1 is expressed as is the historical force vector at time step n; It is the historical parameter at time step n, which indicates the influence of n-1 time steps before the current time step on time step n.

[0059] Assemble the above-mentioned unit stiffness matrix, historical stiffness matrix, and historical parameters into a global form, and let the global stiffness matrix Global history stiffness matrix Global History Vector and The external force vector represents the global form of the external force at time step iT and n+1, where e∈N e , N e is the number of units, and the subscript e is the unit number; the equilibrium equation of the system is expressed as: According to the smoothing equation, the overall displacement vector at each time step is obtained u n and u n+1 represents the overall displacement vector corresponding to time steps n and n+1, are the global forms corresponding to the historical force vector and historical parameters at time step b respectively.

[0060] Reviewing the above content, the design goal of the method of the present invention is to maximize the dissipated work of the structure, which is equivalent to maximizing the mechanical work consumed by the model during the deformation process. The following conditions need to be met:

[0061]

[0062] ρ i =0 or 1; u iT and u iT-1 is the overall displacement vector corresponding to time steps iT and iT-1 obtained according to the overall displacement vector formula of step S4, is the external force vector in global form corresponding to the external force at time step iT-1; ρ i is the density of the i-th unit, parameter v i Represents how much of the total volume the i-th unit occupies; nel is all the units within the specified filter radius with the filter unit as the filter center. Figure 3 and Figure 4 As shown, the filter center is ρ i5 , and its 8 adjacent units together serve as a nel.

[0063] S5; According to the optimization problem definition, the target value of the structure is calculated, and the sensitivity of the target function with respect to the design variables is calculated using the adjoint method.

[0064] According to the definition of the objective function, the target value E is solved by the following formula: The residual R at each time step is iT It is expressed by the following formula: is the global form corresponding to the historical parameters at time step iT-1; introduce the adjoint vector multiplier ψ iT , the objective function is expressed as In order to obtain the sensitivity of the design variable ρ, the objective function is differentiated with respect to the design variable, and the sensitivity expression is obtained as follows: in are the external force vectors in global form corresponding to the external forces at the last two time steps; the adjoint matrix multiplier of the last residual equation is

[0065] Solve 1 to N in descending order t The adjoint matrix multiplier of -1 time step is expressed as follows: where δ k-1，iT is the Cronel function, when k-1 = iT, the function value is 1, otherwise it is 0; the residual term R iT The partial derivative with respect to displacement is Will get and the adjoint vector multiplier ψ iT Substituting into the sensitivity expression, we can obtain the sensitivity of the objective function relative to the design variables.

[0066] S6: After filtering the sensitivity value of each unit through the distance filter, it is expanded into the sensitivity at the control point through the NURBS filter to establish the global sensitivity field. The filtering operation diagram is as follows Figure 4 shown.

[0067] Step S6 includes the following contents:

[0068] S61: Filter and average the sensitivity values ​​so that the updated sensitivity values ​​contain all sensitivity information of the previous iterations, and obtain the global control point sensitivity field.

[0069] Sensitivity filtering of unit sensitivity values ​​can avoid checkerboard phenomenon and grid dependency problems. Filter unit sensitivity values ​​according to the following relationship: Here, nel is based on unit e as the filter center, r min is all cells within the neighborhood of radius r. min is the sensitivity filter radius, the linear weight ω(r ej )=max(0,r min -dis(e,j)], where dis(e,j) is the distance between the centers of cells e and j; α e and α j are the sensitivity values ​​of units e and j respectively; the sensitivity is filtered using a NURBS filter to obtain the sensitivity value α based on the control point i , the specific relationship is as follows: where ρ i , ei and eij represent the variable of the ith control point, the unit set affected by the ith control point and the jth unit in ei, respectively. ij is the NURBS basis function.

[0070] The control point sensitivity values ​​are averaged so that the updated sensitivity values ​​contain all sensitivity information from previous iterations: Where m represents the number of optimizations, thereby establishing a global control point sensitivity field.

[0071] S62: Utilizing the threshold parameters of material addition and removal obtained through the iterative algorithm and implementing the updating of topological variables.

[0072] The specific steps are as follows:

[0073] 1) Threshold parameter for adding materials and the threshold parameter for material removal Satisfies the following relationship: α th Optimize the threshold value for the design variable according to the volume fraction V of the i0+1th iteration i0+1 Determine; for example, if there are 1000 elements in the design domain and α1>α2>...>α 1000 , V i0+1 For a design with 725 elements, α th =α 726 .

[0074] 2) Calculate the allowable volume ratio AR, which is defined as the number of added units divided by the total number of units in the current design; if AR < er, skip step 3); if AR ≥ er, perform step 3) and recalculate and

[0075] 3) Calculate by first sorting the sensitivity values ​​of empty cells 0 The number of cells switched from 0 to 1 will be equal to the evolution ratio er multiplied by the total number of cells in the current design. The evolution ratio er is the maximum volume addition ratio specified. is the sensitivity of the element immediately following the last added unit; then determine Make the volume removed equal to (V i0+1 -V i0 + volume of the added unit), V i0 is the volume fraction at iteration i0. For example, if there are 1000 elements in the design and α1>α2>…>α 1000 , er is 0.02, then the number of units switched from 0 to 1 is 1000×0.02=20. According to the sensitivity order, the sensitivity of the element immediately after the 20th added unit is obtained, and the

[0076] The updating of topological variables is performed according to the following scheme: This scheme shows that when the sensitivity of the entity control point is less than Then delete the control point and restore the sensitivity greater than After the update algorithm converges, the cell density is obtained by interpolation. Finally, it is determined whether the target volume is reached and the iteration converges. If so, the iteration is stopped and the visualization program is entered. Otherwise, the above steps are repeated. In the actual operation, the threshold value is taken in this embodiment. This avoids the complicated calculation process of sensitivity sorting. Figure 5 As shown in the figure, a schematic diagram of an optimization example is provided.

[0077] S7: Update the design variables according to the set evolution ratio parameters;

[0078] S8: Check whether the updated design variables meet the convergence conditions. If they do, output the optimization results; if not, return to step S3.

[0079] like Figure 6 As shown in FIG. 1 , in an implementation example, a clamping beam with a length L to width H ratio of 8:1 and a short-medium force at the bottom is selected, and topology optimization result diagrams are obtained at different times, such as 0.5 seconds, 5 seconds, and 50 seconds.

[0080] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A viscoelastic isogeometric topology optimization method under transient load based on BESO, characterized in that: The steps include: S1: According to the initially set grid size, the model to be designed iso-geometrically meshed to obtain an iso-geometric grid model, the iso-geometric grid model is divided, and the coordinates and numbers of the control points and Gaussian points corresponding to each unit are obtained; the convergence condition and the evolution ratio parameter are determined; S2: Assemble the unit stiffness matrix and unit mass matrix using the coordinates and numbers of the control points and Gaussian points, summarize the stiffness matrix and mass matrix of each unit according to the degree of freedom number, and construct the long-term global stiffness matrix and global mass matrix under the isogeometric mesh model; S3: Determine the optimized target volume according to the evolution ratio parameter and the current design domain volume; S4: Based on the construction of the long-term global stiffness matrix, the tangent stiffness matrix and the historical stiffness matrix are obtained to perform geometric analysis such as viscoelasticity. The structural response at each time step is obtained by solving the equilibrium equation of the viscoelastic material. S5； According to the optimization problem definition, the target value of the structure is calculated, and the sensitivity of the target function with respect to the design variables is calculated using the adjoint method; S6: After filtering the sensitivity value of each unit through the distance filter, it is expanded into the sensitivity at the control point through the NURBS filter to establish a global sensitivity field; S7: Update the design variables according to the set evolution ratio parameters; S8: Check whether the updated design variables meet the convergence conditions. If they do, output the optimization results. If not satisfied, return to step S3.

2. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 1 is characterized in that: The content of step S1 is to use the H-refinement method to refine the design domain according to the preset initial grid size, generate an isogeometric grid model of the specified grid size, divide the isogeometric grid model, determine the coordinates and numbers of the corresponding control points and Gaussian points in each unit according to the characteristics of the isogeometric grid model, and determine the convergence condition, filter radius and evolution ratio parameter according to the requirements, where the relationship of the convergence condition is: τ and N0 represent the allowable convergence error and a positive integer respectively. D Represents the target value under a certain iteration, i0 is the current iteration number, and J∈NO.

3. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 2 is characterized in that: The content of step S2 is that the long-term element stiffness matrix L and the global mass matrix M are solved according to the following relationship: Where N is the NURBS basis function, B is the strain-displacement matrix, N1, N2, ...N n represents the NURBS basis function at different control points, (x, y) is the physical space coordinate, Ω is the physical domain, is the integration domain, D∞ is the long-term elasticity matrix, v is the Poisson's ratio of the material, E ∞ represents the long-term elastic modulus, J1 and J2 represent the transformation relationship from NURBS parameter space to physical space and from integral parameter space to NURBS parameter space, respectively. (ξ,η) are the coordinates in parameter space, is the integral space coordinate, ξ i0+1 and i0 is the coordinate of the adjacent control point along the ξ direction in the parameter coordinate system, η i0+1 and η i0 are the coordinates of the adjacent control points along the η direction in the parameter coordinate system.

4. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 3 is characterized in that: The content of step S3 is to set the volume fraction V of the j0th optimization iteration j0 The calculation formula is: V j0 =max[V req , V j0-1 (1-er)], er represents the evolution ratio, V req represents the target volume, V j0-1 represents the volume fraction of the j0-1th optimization iteration.

5. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 4 is characterized in that: The specific solution steps for the viscoelastic and other geometric analysis described in step S4 are as follows: First, the tangent stiffness matrix k is obtained by using the long-term element stiffness matrix T and the historical stiffness matrix k hist : jp=1,2,...,N p represents the number of terms in the Prony series, γ jp =E jp / E ∞ , Δt is the interval of each time step, τ jp and E jp are the relaxation time and relaxation modulus respectively; for linear viscoelastic materials, the total stress σ n+1 Defined as C ∞ According to the long-term elastic modulus E ∞ The constitutive matrix obtained, time step iT = 1, 2, ..., n, n+1, ..., N t , They represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, respectively. n+1 ,t n-1 ,t iT Indicates the time required to reach the n+1th, n-1th, and iTth time steps; According to the definition of internal force, the internal force is solved by the following formula: Rewrite the above formula as The variables The internal force at time step n+1 is expressed as Each time step ensures the balance of forces, so the internal force is equal to the external force, satisfying is the external force, and the external force at time step n+1 is expressed as is the historical force vector at time step n; is the historical parameter at time step n, indicating the influence of n-1 time steps before the current time step on time step n; Assemble the above-mentioned unit stiffness matrix, historical stiffness matrix, and historical parameters into a global form, and let the global stiffness matrix Global history stiffness matrix Global History Vector and The external force vector represents the global form of the external force at time step iT and n+1, where e∈N e , N e is the number of units, and the subscript e is the unit number; the equilibrium equation of the system is expressed as: According to the smoothing equation, the overall displacement vector at each time step is obtained u n and u n+1 represents the overall displacement vector corresponding to time steps n and n+1, are the global forms corresponding to the historical force vector and historical parameters at time step n respectively.

6. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 5, characterized in that: The content of step S5 is that according to the definition of the objective function, the target value E is solved by the following formula: where u iT and u iT-1 is the overall displacement vector corresponding to time steps iT and iT-1 obtained according to the overall displacement vector formula of step S4, is the external force vector in global form corresponding to the external force at time step iT-1; the residual R at each time step iT It is expressed by the following formula: is the global form corresponding to the historical parameters at time step iT-1; introduce the adjoint vector multiplier ψ iT , the objective function is expressed as In order to obtain the sensitivity of the design variable ρ, the objective function is differentiated with respect to the design variable, and the sensitivity expression is obtained as follows: in are the external force vectors in global form corresponding to the external forces at the last two time steps; the adjoint matrix multiplier of the last residual equation is Solve 1 to N in descending order t The adjoint matrix multiplier of -1 time step is expressed as follows: where δ k-1，iT is the Cronel function, when k-1 = iT, the function value is 1, otherwise it is 0; the residual term R iT The partial derivative with respect to displacement is Will get and the adjoint vector multiplier ψ iT Substituting into the sensitivity expression, we can obtain the sensitivity of the objective function relative to the design variables.

7. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 6, characterized in that: Step S6 includes the following contents: S61: filtering and averaging the sensitivity values ​​so that the updated sensitivity values ​​contain all sensitivity information of previous iterations, and obtaining the sensitivity field of the global control point; S62: Utilizing the threshold parameters of material addition and removal obtained through the iterative algorithm and implementing the updating of topological variables.

8. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 7, characterized in that: The content of step S61 is to filter the unit sensitivity values ​​according to the following relationship: Where nel is based on unit e as the filter center, r min is all cells in the neighborhood of radius, r min is the sensitivity filter radius, the linear weight ω(r ej )=max[0,r min -dis(e,j)], where dis(e,j) is the distance between the centers of cells e and j; α e and α j are the sensitivity values ​​of units e and j respectively; the sensitivity is filtered using a NURBS filter to obtain the sensitivity value based on the control point. The specific relationship is as follows: where ρ i , ei and eij represent the variable of the ith control point, the unit set affected by the ith control point and the ith unit in ei, respectively. ij is a NURBS basis function; the control point sensitivity values ​​are averaged so that the updated sensitivity values ​​contain all sensitivity information of the previous iterations: Where m represents the number of optimizations, thereby establishing a global control point sensitivity field.

9. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 8, characterized in that: The specific steps of the iterative algorithm in step S62 are as follows: 1) Threshold parameter for adding materials and the threshold parameter for material removal Satisfies the following relationship: α th Optimize the threshold value for the design variable according to the volume fraction V of the i0+1th iteration i0+1 Sure; 2) Calculate the allowable volume ratio AR, which is defined as the number of added units divided by the total number of units in the current design; if AR < er, skip step 3); if AR ≥ er, perform step 3) and recalculate and 3) Calculate by first sorting the sensitivity values ​​of empty cells 0 The number of cells that switch from 0 to 1 will be equal to the evolution ratio er multiplied by the total number of cells in the current design, is the sensitivity of the element immediately following the last added unit; then determine Make the volume removed equal to (V i0+1 -V i0 + volume of the added unit), V i0 is the volume fraction of the i0th iteration; The updating of topological variables is performed according to the following scheme:

10. The viscoelastic and other geometric topology optimization method under transient load based on BESO according to claim 9, characterized in that: The value of the allowable convergence error τ is 0.001; the value of the positive integer N0 is 5;

Images