A method for viscoelastic isogeometric topology optimization under transient load based on BESO

Abstract

The application provides a viscoelastic isometric topology optimization method under transient load based on BESO, and belongs to the technical field of isometric topology optimization, and comprises the following steps: isometrically gridding a to-be-designed model to obtain the coordinates of control points and Gauss points corresponding to each unit; long-term global stiffness matrix and global mass matrix are constructed; the target volume of optimization is determined; on the basis of the long-term global stiffness matrix, tangent stiffness matrix and historical stiffness matrix are obtained, and viscoelastic isometric analysis is carried out. According to the definition of the optimization problem, the target value of the structure is calculated, and the sensitivity of the target function with respect to the design variable is calculated by using the adjoint method; after the sensitivity value of each unit is filtered through a distance filter, the sensitivity at the control points is extended through a NURBS filter, and a global sensitivity field is established; and the design variable is updated according to the set evolution ratio parameter. The application can accurately simulate the response of the structure under complex dynamic loading.

Description

Technical Field

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

[0002] With the rapid development of computer technology, topology optimization has been widely used as a powerful tool for achieving innovative designs. Among them, the Bi-directional Evolutionary Structural Optimization (BESO) method is an algorithm used to optimize the material distribution in a structure. Its goal is to find the optimal material distribution within a given space that meets specific performance and constraints, while minimizing material usage through the selection of binary materials.

[0003] Unlike traditional finite element analysis (FEA), isogeometric analysis (IGA) uses non-uniform rational B-splines (NURBS) as a geometric description tool, directly integrating geometric modeling with numerical analysis. This avoids the errors caused by geometric transformation in traditional finite element methods, thereby achieving a seamless transition between computer-aided design (CAD) and computer-aided engineering (CAE). The application of IGA significantly improves the efficiency and accuracy of numerical simulation and design processes, and as a result, isogeometric topology optimization methods based on BESO have garnered widespread attention. However, most current IGA-based topology optimization methods rely on the linear elasticity assumption and are primarily applied to static analysis, ignoring the nonlinear material properties and complex loading conditions commonly encountered in engineering. These methods are particularly limited when dealing with viscoelastic materials. Viscoelastic materials exhibit time-dependent behavior under load, including both transient elastic response and time-dependent viscous flow. Under transient loads, viscoelastic materials exhibit stress relaxation or creep. This time-dependent nonlinear mechanical behavior is crucial in many engineering applications.

[0004] Existing isogeometric topology optimization methods often fail to fully consider these time-dependent properties of viscoelastic materials, resulting in certain limitations in engineering applications involving the stress-strain evolution of structures under transient loads. Therefore, to better address this issue, a BESO-based viscoelastic isogeometric topology optimization method under transient loads is proposed. By analyzing the time-dependent response of the structure through the viscoelastic constitutive relation, it can more accurately simulate the viscoelastic behavior of the structure and optimize the material distribution to meet the complex mechanical requirements of actual engineering. This is very necessary. Summary of the Invention

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

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

[0007] S1: Perform isogeometric meshing on the model to be designed according to the initially set grid size to obtain an isogeometric grid model, divide the isogeometric grid model, obtain the coordinates and numbers of the control points and Gaussian points corresponding to each unit; determine the convergence condition and evolution ratio parameter;

[0008] 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;

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

[0010] S4: Based on 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.

[0011] S5; According to the optimization problem definition, calculate the target value of the structure and use the adjoint method to calculate the sensitivity of the objective function with respect to the design variables;

[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 the 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 so, 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 to 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 requirements, wherein the relationship formula of the convergence condition is: τ and N0 represent the allowed convergence error and positive integer respectively, E 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. (ξ,η) is the coordinate in parameter space, is the integration space coordinate, ξ i0+1 and ξ i0 is the coordinate of the adjacent control points 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] More preferably, the content of step S3 is to set the volume fraction V of the j0th optimization iteration to 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 relaxation time and relaxation modulus respectively; for linear viscoelastic materials, the total stress σ n+1 Defined as C ∞ Based on the long-term elastic modulus E ∞ The constitutive matrix obtained, time step iT=1,2,...,n,n+1,...,N t , Respectively represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, t n+1 , t n-1 , t iT It represents the time required to reach the n+1, n-1, and iT 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; 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;

[0019] Assemble the above element 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 Represents the external force vector in global form corresponding to 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 of 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 in 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, which has a value of 1 when k-1=iT, and 0 otherwise; the residual term R iT The partial derivative with respect to displacement is Will get and the adjoint vector multiplier ψ iT Substituting it into the sensitivity expression, we can get the sensitivity of the objective function relative to the design variables.

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

[0023] 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;

[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 the filter center of unit e, r min For all cells in the neighborhood of radius rmin is the sensitivity filter radius, linear weight ω(r ej ) = max[0, r min -dis(e, j)], dis(e, j) is the distance between the center of element e and element j; a e and a j are the sensitivity values of elements e and j, respectively; the sensitivity is filtered by a NURBS filter to obtain the control point-based sensitivity value, and the specific relationship is as follows: where p i , e and eij represent the variable of the ith control point, the set of elements affected by the ith control point, and the jth element in e, respectively, N ij is the NURBS basis function; the control point sensitivity value is averaged to make the updated sensitivity value contain all the sensitivity information of the previous iteration: where m represents the number of optimizations, thereby establishing a global control point sensitivity field.

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

[0027] 1) Let the threshold parameter of material addition a and the threshold parameter of material removal a satisfy the following relationship: a th is the design variable optimization threshold, which is determined according to the volume fraction V i0+1 of the ith0+1 iteration;

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

[0029] 3) Calculate a by first sorting the sensitivity values of empty elements 0; the number of elements switched from 0 to 1 will be equal to the evolution ratio er multiplied by the total number of elements in the current design, is the sensitivity of the element immediately after the last added element; then determine a so that the removed volume is equal to (V i0+1 -V i0 + the volume of the added elements), V i0 is the volume fraction of the ith iteration;

[0030] The update of the topological variable is carried out 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 advantages over 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 the geometric conversion errors in traditional finite element methods, and can more accurately simulate 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 in the optimization process. This has broad application prospects for engineering applications involving long-term loads or periodic loading, such as bridges, building structures, and mechanical equipment.

[0034] (2) Unlike traditional topology optimization methods based on linear elasticity 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 optimal performance under complex loading 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 common grid dependencies and checkerboard patterns but also ensures that each iteration fully considers the previous optimization information by averaging historical data. Especially for topology optimization of viscoelastic materials, due to the complexity of transient responses, this processing method effectively improves the stability and convergence speed of sensitivity calculations, ensuring 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 following briefly introduces the drawings required for use in the embodiments or the description of the prior art. 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 any creative work.

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

[0038] Figure 2 Schematic diagram of a Kelvin-Voigt constitutive model of viscoelastic materials based on a BESO-based geometric topology optimization method for viscoelasticity 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 loads based on BESO of the present invention;

[0040] Figure 4 Schematic diagram of density filtering and NURBS filtering constructed in an embodiment of a BESO-based viscoelastic and other geometric topology optimization method under transient loads of the present invention;

[0041] Figure 5 A schematic diagram of an optimization example obtained by a BESO-based method in a BESO-based method for geometric topology optimization of viscoelasticity under transient loads according to the present invention;

[0042] Figure 6 This is a schematic diagram of optimization examples under different simulation times obtained in an embodiment of a geometric topology optimization method based on BESO under transient loads. 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 embodiments described 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 making creative efforts 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 of viscoelasticity under transient load based on BESO, comprising the following steps:

[0045] S1: Perform isogeometric meshing on the model to be designed according to the initially set grid size to obtain an isogeometric grid model, divide the isogeometric grid model, obtain the coordinates and numbers of the control points and Gaussian points corresponding to each unit; determine the convergence condition and evolution ratio parameter;

[0046] This method adopts the H-refinement method, which generates an isogeometric grid model of a specified grid size by refining the design domain according to the preset initial grid size, divides the isogeometric grid 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 grid model. Figure 3 As shown. Determine the convergence condition, filter radius, evolution ratio and other parameters according to actual needs. The relationship between the convergence condition is: τ and N0 represent the allowed convergence error and positive integer respectively, E D represents the target value at a certain iteration, i0 is the current iteration number, and J∈N0. In this embodiment, the allowed convergence error τ is 0.001; the positive integer N0 is 5.

[0047] like Figure 2 FIG. 4 is a schematic diagram of the Kelvin-Voigt constitutive model of the viscoelastic material selected in 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. (ξ,η) is the coordinate in parameter space, is the integration space coordinate, ξ i0+1 and ξ i0 is the coordinate of the adjacent control points 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 based on 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 represents the volume fraction at the j0-1 optimization iteration. It should be noted 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 only changes the topology; the volume fraction will always remain constant.

[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. By solving the equilibrium equations of the viscoelastic material, the structural response at each time step is obtained.

[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[松弛时间τ jp Unit: seconds]]> <![CDATA[松弛模量E jp Unit 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 linear viscoelastic materials, the total stress σ n+1 Defined as C ∞ Based on the long-term elastic modulus E ∞ The constitutive matrix obtained, time step iT = 1, 2, ..., n, n+1, ..., N t , Respectively represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, t n+1 , t n-1 , t iT It represents the time required to reach the n+1, n-1, and iT 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 element 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 Represents the external force vector in global form corresponding to 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 of 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, 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 in 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 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, which has a value of 1 when k-1=iT, and 0 otherwise; the residual term R iT The partial derivative with respect to displacement is Will get and the adjoint vector multiplier ψ iT Substituting it into the sensitivity expression, we can get 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 cell sensitivity values ​​can avoid checkerboard phenomenon and grid dependency problems. Cell sensitivity filtering is performed according to the following relationship: Here, nel is based on unit e as the filter center, r min 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 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 the sensitivity information from the previous iterations: Where m represents the number of optimizations to establish the 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 material 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, and the sensitivity of the element immediately after the 20th added unit is obtained according to the sensitivity order, so as to determine

[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 process, the threshold value is taken in this embodiment. It 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 so, output the optimization results; if not, return to step S3.

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

[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 scope of protection of the present invention.

Claims

1. A BESO-based viscoelastic isogeometric topology optimization method under transient loads, characterized in that: The steps include: S1: Perform isogeometric meshing on the model to be designed according to the initially set grid size to obtain an isogeometric grid model, divide the isogeometric grid model, obtain the coordinates and numbers of the control points and Gaussian points corresponding to each unit; determine the convergence condition and evolution ratio parameter; 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 grid model; 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 integration space coordinate, ξ i0+1 and ξ i0 is the coordinate of the adjacent control points 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; S3: Determine the optimized target volume according to the evolution ratio parameter and the current design domain volume; S4: Based on 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. The specific solution steps for the viscoelastic and other geometric analysis in step S4 are as follows: First, the tangent stiffness matrix k is obtained by 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 relaxation time and relaxation modulus respectively; for linear viscoelastic materials, the total stress σ in the n+1th iteration step n+1 Defined as C ∞ Based on the long-term elastic modulus E ∞ The constitutive matrix obtained by solution, time step iT=1,2,...,n,n+1,...,Nt, Respectively represent the unit displacement vectors corresponding to time steps iT-1, iT, n+1 and n, t n+1 , t n-1 , t iT Indicates the time required to reach the n+1, n-1, and iT 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; 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; Assemble the above element 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 of 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; S5; According to the optimization problem definition, calculate the target value of the structure and use the adjoint method to calculate the sensitivity of the objective function with respect to the design variables; 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: Among them, 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 in 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, which has a value of 1 when k-1=iT, and 0 otherwise; 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 get the sensitivity of the objective function relative to the design variable; 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; S7: Update the design variables according to the set evolution ratio parameters; S8: Check whether the updated design variables meet the convergence conditions. If so, output the optimization results; if not, return to step S3.

2. The viscoelastic and geometric topology optimization method under transient load based on BESO according to claim 1, characterized in that: Step S1 is to use the H-refinement method to refine the design domain according to the preset initial grid size to generate an isogeometric grid model of the 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. The convergence condition, filter radius, and evolution ratio parameters are determined according to the requirements. The relationship between the convergence condition is: τ and N0 represent the allowed convergence error and positive integer respectively, E D Represents the target value under a certain iteration, i0 is the current iteration number, J∈N0.

3. The viscoelastic and geometric topology optimization method under transient load based on BESO according to claim 1 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.

4. The viscoelastic and geometric topology optimization method under transient load based on BESO according to claim 1 is characterized in that: Step S6 includes the following contents: 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; S62: Utilizing the threshold parameters of material addition and removal obtained through the iterative algorithm and implementing the updating of topological variables.

5. The viscoelastic and geometric topology optimization method under transient load based on BESO according to claim 4 is characterized in that: The content of step S61 is to filter the unit sensitivity values ​​according to the following relationship: Where nel is the filter center of unit e, r min For 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 jth unit in ei, respectively. ij is the NURBS basis function, α i is the sensitivity value corresponding to the i-th control point; the control point sensitivity values ​​are averaged so that the updated sensitivity value contains all the sensitivity information of the previous iteration: Where m represents the number of optimizations to establish the global control point sensitivity field.

6. The viscoelastic and geometric topology optimization method under transient load based on BESO according to claim 5, characterized in that: The specific steps of the iterative algorithm in step S62 are as follows: 1) Threshold parameter for adding material 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 ​​where the control point density is 0 The number of control points that switch from 0 to 1 will be equal to the evolution ratio er multiplied by the total number of control points in the current design, is the sensitivity of the element immediately following the last added control point; then determine Make the volume removed equal to (V i0+1 -V i0 + added control point volume), V i0 is the volume fraction of the i0th iteration; The updating of topological variables is performed according to the following scheme:

7. The BESO-based viscoelastic and other geometric topology optimization method under transient loads according to claim 6, characterized in that: The allowed convergence error τ is 0.001; the positive integer N0 is 5;

Images