Efficient fatigue topological optimization method for isogeometric analysis

Abstract

The invention provides an isogeometric analysis-oriented high-efficiency fatigue topological optimization method, which belongs to the technical field of geometric analysis and mainly comprises the following steps of: constructing a fatigue problem model, dividing the number of layers of a multi-layer grid method and the size of a design domain according to the requirements of the design domain, and generating an isogeometric node vector of the lowest layer by using h insertion; according to the isogeometric node vector, assembling each level of interpolation operator and continuation operator of the multi-grid method; decomposing and storing independent load components in the load spectrum, and solving a combined weight of each independent load component; an overall stress spectrum is obtained according to weight assembly; evaluating the fatigue constraint of the fatigue damage of the unit; approximately expressing the fatigue constraint topological optimization problem according to an augmented Lagrangian method; and performing updating iteration of the isogeometric analysis design variable until the convergence of the fatigue constraint topological optimization problem is judged according to the design variable variation. According to the scheme, the simplified calculation and counting problems of the multi-cycle general fatigue load spectrum are improved, and the solving cost is reduced.

Description

Technical Field

[0001] The invention relates to the technical field of geometric analysis, and in particular to an efficient fatigue topology optimization method oriented to isogeometric analysis. Background Art

[0002] With the rapid development of modern industrial technology, the demand for lightweight, long life and high reliability in engineering optimization design continues to increase, and structural fatigue has gradually become an important factor in engineering practice that cannot be ignored. Fatigue refers to the phenomenon that a structure gradually forms damage and fails after being subjected to enough cyclic disturbances. Even if the stress level is lower than the yield strength of the material, long-term cyclic disturbances can cause the structure to break or fail. Fatigue failure problems are widely present in the fields of machinery, aerospace, automobiles, etc., and have a decisive impact on the safety and service life of the structure.

[0003] Fatigue topology optimization is a technique that maximizes the fatigue resistance of a structure under cyclic disturbances by optimizing material distribution. Unlike static loads, variable disturbance analysis in fatigue problems requires multiple cyclic analyses. Due to the complexity of the physical damage accumulation process, the computational cost of fatigue analysis is usually much higher than that of static mechanical analysis, and this problem is more serious under irregular general cyclic disturbances. However, the traditional finite element analysis (FEA) method is limited by the model and calculation accuracy, and it is difficult to take into account both computational efficiency and accuracy. In recent years, isogeometric analysis (Isogeometric Analysis, IGA) as an emerging numerical analysis method has gradually been applied to the field of structural optimization. IGA directly uses CAD (Computer Aided Design, CAD) geometric models for analysis, avoiding the conversion steps between geometric models and analysis models, thereby improving structural continuity and geometric accuracy.

[0004] Although IGA has significant advantages in geometric accuracy and continuity, it still faces challenges in dealing with fatigue topology optimization problems. Fatigue analysis requires accurate prediction of the cumulative damage of the structure under cyclic loads, which usually requires a large amount of simulation calculations. For non-periodic loads and even general loads that are closer to engineering practice, the amount of calculation for fatigue analysis will increase further. Due to the high-order continuity of isogeometric analysis, the stiffness matrix involved in the optimization calculation is larger, which further increases the computational cost in fatigue analysis. Therefore, it is very necessary to provide an efficient fatigue topology optimization method for isogeometric analysis to simplify the calculation and counting problems of multi-cycle general fatigue load spectra and reduce the cost of matrix solution under large-scale fatigue problems. Summary of the invention

[0005] In view of this, the present invention proposes an efficient fatigue topology optimization method for isogeometric analysis, which effectively reduces the amount of calculation, improves calculation efficiency, and reduces calculation cost.

[0006] The present invention provides an efficient fatigue topology optimization method for isogeometric analysis, comprising the following steps:

[0007] Step 1: Construct fatigue problem model;

[0008] Step 2: Divide the number of multi-layer grid layers and the size of the design domain according to the design domain requirements, and use h insertion to generate the lowest level of isogeometric node vectors;

[0009] Step 3: Assemble the interpolation operators and extension operators of each level of the multigrid method according to the isogeometric node vectors;

[0010] Step 4: Decompose and store the independent load components in the load spectrum, and solve the combined weights of the independent load components;

[0011] Step 5: Perform isogeometric analysis to solve each independent load component;

[0012] Step 6: Calculate the stress spectrum of each unit of the independent load component, assemble it according to the weight to get the overall stress spectrum, and then perform the signed von-mises stress correction of the unit stress;

[0013] Step 7: The peak and valley values ​​of each cycle and their indexes are extracted from the modified stress spectrum of each unit using the rain flow counting method, and the average stress and stress amplitude of the unit in each cycle are calculated, and then the Morrow correction is performed to obtain the equivalent stress amplitude;

[0014] Step 8: Evaluate fatigue constraints of element fatigue damage based on the Palmgren-Miner linear cumulative damage model;

[0015] Step 9: Approximate the fatigue constrained topology optimization problem according to the augmented Lagrangian method;

[0016] Step 10: Repeat steps 3 to 9 to iterate the update of the isogeometric analysis design variables until the fatigue constrained topology optimization problem converges according to the change in the design variables;

[0017] Step 11: According to the preset number of layers of the multigrid method, the design variables on the coarse grid are mapped to the fine grid using the continuation operator P, the design domain is updated for solution, and steps 3 to 10 are repeated;

[0018] Step 12: Repeat steps 3 to 11 until the preset final level is reached and output the final solution result.

[0019] Based on the above technical solution, preferably, the specific content of step 1 is: i is the numbered i-th design variable, and the design variable ρ i After distance smoothing, the smoothed design variables are obtained Then the design variables Obtaining projection design variables through Heaviside projection The smoothing process expression is: W i，j = {r min -dist(i,j)},{j∈N i |dist(i,j)≤r min}, ρ j is related to the design variable ρ i The set of design variables within a certain Euclidean distance, N i is the total number of design variables in the design variable set, W i，j is the corresponding weight coefficient of the j-th design variable to the i-th design variable, r min is the filter radius, dist(i, j) is the design variable ρ i and ρ j The Euclidean distance between them; the expression of the Heaviside projection process is Where β and η are the steepness and threshold of the Heaviside projection respectively. The isogeometric analysis topology optimization model is established, and the corresponding mathematical expression is as follows: findρ=[ρ1,ρ2,…,ρ n ], subject to:g i (ρ)=D i -1<0, i=1,...,Ne, 0≤ρ e ≤1, e=1, ..., Ne, t=1,……NF,where ρ is the initial design variable vector, v e is the area or volume corresponding to the e-th design variable component, Ne is the total number of elements, NF is the total number of moments in the load spectrum, and D i is the damage variable of the ith unit, g i (ρ) is the fatigue constraint for the ith element, For projection-based design variables Assemble the obtained global stiffness matrix, U t and F t are the global displacement and load at the tth moment.

[0020] Preferably, the specific content of step 3 is: for the case where the topology optimization analysis domain is a single dimension, assuming that there are two sets of node vectors and The subscripts i and j represent the i-th and j-th numbers of the node vectors. When used as symbols, they are a general reference to a set of sequences. The superscript H represents a coarse grid, and the superscript h represents a fine grid. Depend on After h insertion, the NURBS basis function obtained by two sets of node vectors and Satisfies the following linear transformation relationship: Where k is the repeatability of the B-spline, n h is the number of nodes in the grid, r = 1, 2, ..., k, r is the extension matrix The number of layers, i, j are used as symbols to refer to the extension matrix With the intermediate interpolation matrix The rows and columns of , in specific calculations, represent the position of the matrix element in the i-th row and j-th column. The following solution relationship is satisfied for r:

[0021] Obtain the restriction matrix R for transforming from fine grid to coarse grid, and use the transformation rank of the restriction matrix as the interpolation matrix P = R for transforming from coarse grid to fine grid T ;

[0022] For the case where the topology optimization analysis domain is two-dimensional, for the coarse mesh node vector and The corresponding B-spline basis function is and The fine grid node vector obtained after conversion is and The corresponding B-spline basis function is and The restriction matrices for the conversion from fine grid to coarse grid in two dimensions are R 1，i，s and R 2，j，t , i, s, j, t represent the element positions in the two-dimensional matrix respectively, and the two-dimensional transformation relationship is obtained: R a，b =R 1，i，s R 2，j，t , a=i+(j-1)n, b=s+(t-1)n′, n is R 1，i，s The number of coarse grid control points of the corresponding dimension, n′ is R 2，j,t The number of coarse grid control points of the corresponding dimension; R 1，i，s and R 2，j，t The transpose of is used as the interpolation matrix for transforming from coarse grid to fine grid;

[0023] The data u corresponding to the coarse grid H The data u corresponding to the fine grid h The conversion expression is abbreviated as: H=Ru h ,u h =Pu H . .

[0024] Further preferably, the specific content of step 4 is: the form of decomposing the load spectrum is as follows: F t =a t *F a +b t *F b +…+n t *F n , where F t is the load component at the tth load spectrum moment, F a 、F b , ..., F n is the independent load component obtained by decomposition, a t 、b t ,...,n t It is the combined weight of each independent load component, and the combined weight is stored.

[0025] Further preferably, the specific content of step 6 is: for the independent load components solved in step 5, assemble them according to the stored combination weights, and the assembled expression is as follows:

[0026]

[0027] is the combined weight of unit e, is the independent load component that has been solved; then the stress spectrum assembled for each independent load component that has been solved is Perform a signed stress von-mises correction, specifically, a signed stress von-mises correction Δσ e (t) VM The expression is as follows: Δσ e (t) = σ e (t)-σ e (t0), where t0 is the loading start time; sign(·) is the sign function, which takes the value of 1 or -1; matrix W = [1 1 0] T ; V is the stress coefficient matrix, σ e (t) is the stress vector at element e, q is the stress penalty parameter, is the design variable ρ of unit e e The result after smoothing and Heaviside projection, B e and U e are the strain displacement matrix and displacement matrix of unit e respectively, D0 is the elastic matrix of solid material, E0 is Young's modulus of the solid material, and v is Poisson's ratio.

[0028] More preferably, the specific content of step 7 is: for the reorganized and corrected stress spectrum, the expressions of the average stress and stress amplitude of unit e in the i-th cycle are:

[0029]

[0030] and are the maximum stress and minimum stress of unit e in the i-th cycle respectively; Morrow correction is then performed to obtain the equivalent stress amplitude The expression is σ′ f is the fatigue strength coefficient of the material.

[0031] Further preferably, the specific content of step 8 is: the unit fatigue damage D must satisfy the following fatigue constraints: Where N f，i is the number of fatigue life cycles corresponding to the stress of the i-th reorganized and corrected stress spectrum; n i is the cycle number corresponding to the stress of the i-th reorganized and corrected stress spectrum; b is the fatigue strength index.

[0032] Further preferably, the specific content of step 9 is: according to the augmented Lagrangian method, the model of the fatigue constrained topology optimization problem is approximated as the following expression: findρ=[ρ1,ρ2,…,ρ n ], subject to:g i (ρ)=D i -1, 0≤ρ e ≤1, in, is the global stiffness matrix, which is abbreviated as K in the following, i.e. KU t =F t , where K is the abbreviated global stiffness matrix, λ e and μ e is the Lagrange multiplier corresponding to the fatigue constraint of the e-th unit, and the update criterion of the Lagrange multiplier is:

[0033] μ e,l+1 =max(μ e,l *1.1,10000),λ e,l+1 =λ e，l +μ e,l *g e (ρ) where λ e，l and μ e，l is the Lagrange multiplier at the lth iteration, μ e，l+1 and λe，l+1 is the Lagrange multiplier at the l+1th iteration.

[0034] Further preferably, the fatigue constraint is integrated into the objective function to solve the sensitivity of the objective function to the design variables. The mathematical expression of the sensitivity is as follows:

[0035]

[0036] According to the chain rule, The sensitivity expressions for the design variables are: After transformation, we get

[0037] For the single-corrected element stress σ extracted by the rainflow counting method e，i (t) The expression of the design variable is:

[0038] where σ e，i (t) is the overall stress state of the e element at time t (maximum or minimum stress moment) of the i-th rainflow cycle. For the two-dimensional case, σ e，i (t) contains three components, namely They are the stress component in the x direction, the stress component in the y direction, and the shear stress on the xy plane;

[0039]

[0040] Where L e A 0-1 sparse matrix is ​​selected so that U e =L e U, U is the global displacement.

[0041] Further preferably, the sensitivity of the global displacement U to the design variables is solved by the adjoint method, and the expression is: F is the global load, Update the expression for sensitivity to:

[0042] in and is the adjoint vector λ e The maximum and minimum values ​​of U max and U min are the maximum and minimum values ​​of the global displacement U.

[0043] The invention provides an efficient fatigue topology optimization method for isogeometric analysis, which has the following beneficial effects compared with the prior art:

[0044] (1) The present invention applies a multi-layer grid method based on isogeometry, utilizes the node refinement relationship of isogeometry, realizes isogeometry density information mapping, enables fine grid iteration to start directly from a more accurate intermediate optimization structure, thereby shortening the number of convergence iterations; it also adopts an independent linear component decomposition method to solve the node displacement, avoids the linear increase of the solution cost with the number of cycles, significantly reduces the redundant displacement calculation amount of the general load spectrum, and thereby improves the efficiency of structural topology optimization;

[0045] (2) The present invention proposes a multi-grid conjugate gradient method based on isogeometry, and applies the augmented Lagrangian method according to the actual solution process of the fatigue problem. It significantly reduces the amount of calculation while ensuring the accuracy of the optimization results, greatly improves the computational efficiency of inverting large-scale large condition number matrices, and makes it possible to solve large-scale fatigue problems considering general loads. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] 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.

[0047] Figure 1 A flowchart of the steps of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention;

[0048] Figure 2 A schematic diagram of one-dimensional and two-dimensional NURBS basis functions of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention;

[0049] Figure 3 A mapping schematic diagram of a multi-layer grid method of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention;

[0050] Figure 4 A schematic diagram of the stress-time history of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention.

[0051] Figure 5 A schematic diagram of rainflow reversal statistics of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention.

[0052] Figure 6 A histogram diagram of rain flow statistics for an efficient fatigue topology optimization method for isogeometric analysis according to the present invention.

[0053] Figure 7It is a schematic diagram of the SN fatigue strength curve of an efficient fatigue topology optimization method for isogeometric analysis of the present invention.

[0054] Figure 8 A schematic diagram of the Palmgren-Miner rule statistical cumulative damage of an efficient fatigue topology optimization method for isogeometric analysis of the present invention. DETAILED DESCRIPTION

[0055] 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.

[0056] like Figure 1 As shown, the present invention provides an efficient fatigue topology optimization method for isogeometric analysis, comprising the following steps:

[0057] Step 1: Construct a fatigue problem model.

[0058] Here, an isogeometric topology optimization model for fatigue optimization problem is constructed with the goal of minimizing material usage and the fatigue damage at each unit of the structure as a constraint.

[0059] Specifically, the isogeometric analysis model based on the variable density method is used to establish the optimal design variable ρ and the structural density field m v (ρ) and Young's modulus field m E (ρ), the corresponding relationship between the three is introduced below: ρ = [ρ 1，1 , ρ 1，2 ,…,ρ i，j ,…,ρ n0，m0 ], m E (ρ)=E mim +(E0-E min )m v (ρ) p0 , where i and j represent the control point numbers along the parameterized coordinates ξ0 and η0, ρ i，j is the design variable corresponding to the i-th and j-th control points in the corresponding direction, represents the NURBS basis function corresponding to the i-th and j-th control points, where r and q are the repetitions along the parameterized coordinates ξ0 and η0, respectively, and n0 and m0 are the number of control points along the parameterized coordinates ξ0 and η0, respectively. is the design variable after distance smoothing and heaviside projection processing; E0 is the Young's modulus of the solid material, Emin is the given minimum Young's modulus to avoid singular phenomena in the solution process; the superscript p0 is the penalty parameter of the variable density method, which is set to 3.

[0060] Considering that the principle of isogeometric multigrid assembly will be introduced later, the assembly method of NURBS basis functions is given below. For a non-decreasing real number sequence vector U = [u1, u2, ..., u n+p+1 ], where p is the degree and n is the number of control nodes, which is also the number of basis functions. The basis function B of the B-spline can be obtained by the Cox-de Boor recursion formula i，p The expression of (u) is as follows:

[0061] When p = 0,

[0062] When p>0,

[0063] In the above formula, the subscript i is the serial number of the corresponding vector sequence, and the maximum value is the number of nodes in the vector.

[0064] NURBS basis function is an extended form of B-spline basis function. The positive weight w is introduced into the univariate B-spline basis function. i , we can get a univariate NURBS basis function: Correspondingly, according to the tensor product property, the two-dimensional or multi-dimensional NURBS basis function N i，p (u) is constructed as: In the formula represents parameter coordinates, where i0 = [i1, i2, ..., i dp ] represents the index position in the tensor product, p = [p1, p2, ..., p dp ] represents the order of each parameter direction. On the right side of the formula, is the univariate NURBS basis function in the mth parameter direction. dp represents the dimension of the parameter space, is a dp univariate knot vector, n m is the number of control points, p m is the order along the parameter direction m. Figure 2 As shown, for the node vector ξ0 = η0 = [0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1], the values ​​of the two-dimensional NURB basis function form the corresponding surface.

[0065] ρ i is the numbered i-th design variable, and the design variable ρ i After distance smoothing, we get Then the smoothed design variables Obtained via Heaviside projection The smoothing expression is W i，j = {r min -dist(i, j)},{j∈N i |dist(i,j)≤r min}, ρ j is related to the design variable ρ i The set of design variables within a certain distance, N i is the total number of design variables in the design variable set, W i，j is the corresponding weight coefficient of the j-th design variable to the i-th design variable, r min is the filter radius, dist(i, j) is the design variable ρ i and ρ j The Euclidean distance between them; the expression of Heaviside projection is Where β and η are the steepness and threshold of the Heaviside projection, respectively. The larger the β, the more effective it is in removing grayscale elements. However, a larger β value will lead to poor convergence of topology optimization. In this embodiment, β is set to 5 and η is set to 0.5; an isogeometric analysis topology optimization model is established, and the corresponding mathematical expression is as follows: findρ=[ρ1,ρ2,…,ρ n ], subject to: g i (ρ)=D i -1<0, i=1,..., Ne, 0≤ρ e ≤1, e=1, ..., Ne, t=1,……,NF,where ρ is the design variable vector, v e is the area or volume corresponding to the e-th design variable component, Ne is the total number of elements, NF is the total number of moments in the load spectrum, and D i is the damage variable of the ith unit, g i (ρ) is the fatigue constraint of the ith element, is the global stiffness matrix, U t and F t are the global displacement and load at the tth moment.

[0066] Step 2: Divide the number of multi-layer grid layers and the size of the design domain according to the design domain requirements, and use h insertion to generate the lowest level of equal geometry node vectors.

[0067] Taking the two-dimensional design domain as an example, limited by the basic principle of h insertion, assuming that the number of layers of the multi-layer grid method is L and the preset design domain size is n×m, the design domain of each layer of grid satisfies the following expression:

[0068] n c =n / 2 c-1 ;mc =m / 2 c-1 ;n c ∈N,m c ∈N; c = 1, ..., L; where n c and m c is the grid design domain size of the cth level and is a positive integer. For the calculated design domain size at the lowest level, h nodes are inserted with equal geometry, and the node vector expressions in the ξ0 and η0 directions are obtained as follows: Here, the subscript i represents that a node corresponds to the i-th level, the subscript r is the repetition degree of ξ0, and the subscript q is the repetition degree of η0.

[0069] Step 3: According to the isogeometric node vectors, assemble the interpolation operators and extension operators of each level of the multigrid method.

[0070] The multigrid method itself is based on h insertion, and the mapping of node information is realized according to the relationship between control points. The following implementation method gives the assembly method of the interpolation operator and the extension operator:

[0071] For the case where the topology optimization analysis domain is single-dimensional, it is assumed that there are two sets of node vectors and The subscripts i and j represent the i-th and j-th numbers of the node vectors. When used as symbols, they are a general reference to a set of sequences. The superscript H represents a coarse grid, and the superscript h represents a fine grid. Depend on After h insertion, the NURBS basis function obtained by two sets of node vectors and Satisfies the following linear transformation relationship: Where k is the repeatability of the B-spline, n h is the number of nodes in the grid, r = 1, 2, ..., k, r is the extension matrix The number of layers, i, j are used as symbols to refer to the extension matrix With the intermediate interpolation matrix The rows and columns of , in specific calculations, represent the position of the matrix element in the i-th row and j-th column. The following solution relationship is satisfied for r: Obtain the restriction matrix R for transforming from fine grid to coarse grid, and use the transpose of the restriction matrix as the interpolation matrix P = R for transforming from coarse grid to fine grid T ;

[0072] For the case where the topology optimization analysis domain is two-dimensional, for the coarse mesh node vector and The corresponding B-spline basis function is and The fine grid node vector obtained after conversion is and The corresponding B-spline basis function is and The restriction matrices for the conversion from fine grid to coarse grid in two dimensions are R 1，i，s and R 2，j，t , i, s, j, t represent the element positions in the two-dimensional matrix respectively, and the two-dimensional transformation relationship is obtained: R a，b =R 1，i，s R 2，j，t , a=i+(j-1)n, b=s+(t-1)n′, n is R 1，i，s The number of coarse grid control points of the corresponding dimension, n′ is R 2，j，t The number of coarse grid control points of the corresponding dimension; R 1，i，s and R 2，j，t The transpose of is used as the interpolation matrix for transforming from coarse grid to fine grid;

[0073] The data u corresponding to the coarse grid H The data u corresponding to the fine grid h The conversion expression is abbreviated as: H =Ru h ,u h =Pu H . .

[0074] Step 4: Decompose and store the independent load components in the load spectrum, and solve the combined weights of the independent load components.

[0075] The decomposed load spectrum is in the form of: F t =a t *F a +b t *F b +…+n t *F n , where F t is the load component at the tth load spectrum moment, F a 、F b , ..., F n is the independent load component obtained by decomposition, a t 、b t ,...,n t It is the combined weight of each independent load component, and the combined weight is stored.

[0076] Step 5: Perform isogeometric analysis to solve each independent load component.

[0077] For the global displacement The solution is solved by assembling a multi-grid preprocessing matrix using interpolation operators and extension operators. The first iteration solves the problem by reusing the last solution result to achieve accelerated iterative convergence. Multi-grid is a multi-level iterative method for solving linear equations. This method divides the computational domain into grids of different sizes and solves the problem on the coarsest grid of the smallest size. The V-loop solution process of the multi-grid conjugate gradient method with reused initial values ​​is as follows:

[0078] Input: F 1 ,K 1 ,L,Upre,R,P

[0079] U 1 =Upre

[0080] for l=1:L-1

[0081] smooth(K l ,U l ,F l )

[0082] r l =F l -K l U l

[0083] F l+1 =P T r l

[0084] K l+1 =RK l P

[0085] U l+1 =0

[0086] end for

[0087] U l =K l \F l

[0088] for l=1:L-1

[0089] U l =U l +PU l+1

[0090] smooth(K l ,U l ,F l )

[0091] end for

[0092] Output: U 1

[0093] Where Upre is the last iterative solution stored, K l , U l 、F l is the stiffness matrix, displacement vector, and force vector of the lth layer, and the damped Jacobian method is used in the smoothing process.

[0094] Step 6: Calculate the stress spectrum of each unit of the independent load component, assemble it according to the weight to get the overall stress spectrum, and then perform the signed von-mises stress correction on the unit stress.

[0095] The specific content is: for the independent load components solved in step 5, assemble them according to the stored combination weights. The assembled expression is as follows: a e，t 、b e，t ,...,n e，t is the combined weight of unit e, is the independent load component that has been solved; then the stress spectrum assembled for each independent load component that has been solved is Perform a signed stress von-mises correction, specifically, a signed stress von-mises correction Δσ e (t) VM The expression is as follows: Δσ e (t) = σ e (t)-σ e (t0), where t0 is the loading start time; sign(·) is the sign function, which takes the value of 1 or -1; matrix W = [1 1 0] T ; V is the stress coefficient matrix, σ e (t) is the stress vector at element e, q is the stress penalty parameter, is the design variable ρ of unit e e The result after smoothing and Heaviside projection, B e and U e are the strain displacement matrix and displacement matrix of unit e respectively, D0 is the elastic matrix of solid material, E0 is Young's modulus of the solid material, and v is Poisson's ratio.

[0096] Step 7: The peak and valley values ​​of each cycle and their indexes are extracted from the modified stress spectrum of each unit using the rain flow counting method, and the average stress and stress amplitude of the unit in each cycle are calculated, and then the Morrow correction is performed to obtain the equivalent stress amplitude.

[0097] Figure 4 The stress time history shown is Figure 5This is the preliminary completed rain flow reversal extreme point diagram, the horizontal and vertical coordinates are the extracted rain flow reversal time and the corresponding reversal extreme stress, Figure 6 The rain flow counting results are given in the form of a histogram, and the relevant parameters of the rain flow reversal obtained statistically are given. Figure 7 The following is a schematic diagram of the applied S-N curve range. The specific content is: for the reorganized and corrected stress spectrum, the expressions of the average stress and stress amplitude of unit e in the i-th cycle are: and are the maximum stress and minimum stress of unit e in the i-th cycle respectively; Morrow correction is then performed to obtain the equivalent stress amplitude The expression is σ′ f is the fatigue strength coefficient of the material.

[0098] Step 8: Evaluate the fatigue constraint of the unit fatigue damage based on the Palmgren-Miner linear cumulative damage model. The statistical results are as follows Figure 8 As shown, the calculated cumulative unit damage is given.

[0099] The specific content is: The unit fatigue damage D must satisfy the following fatigue constraints: Where N f，i is the number of fatigue life cycles corresponding to the stress of the i-th reorganized and corrected stress spectrum; n i is the cycle number corresponding to the stress of the i-th reorganized and corrected stress spectrum; b is the fatigue strength index.

[0100] Step 9: Approximate the fatigue constrained topology optimization problem according to the augmented Lagrangian method.

[0101] The specific content is: According to the augmented Lagrangian method, the model of fatigue constrained topology optimization problem is approximated as follows: findρ=[ρ1,ρ2,…,ρ n ], subjectto: g i (ρ)=D i -1, 0≤ρ e ≤1, in, is the global stiffness matrix, which is abbreviated as K in the following, i.e. KU t =F t , where K is the abbreviated global stiffness matrix, λ e and μ e is the Lagrange multiplier corresponding to the fatigue constraint of the e-th unit. According to the KKT condition, the update criterion of the Lagrange multiplier is:

[0102] μe,l+1 =max(μ e,l *1.1,10000),λ e,l+1 =λ e，l +μ e,l *g e (ρ), where λ e，l and μ e，l is the Lagrange multiplier at the lth iteration, μ e，l+1 and λ e，l+1 is the Lagrange multiplier at the l+1th iteration.

[0103] In this implementation method, the constraint function is integrated into the objective function using the augmented Lagrangian method, that is, it is necessary to solve the sensitivity of the objective function to the design variables. The mathematical expression of the sensitivity is as follows:

[0104]

[0105] According to the chain rule, The sensitivity expressions for the design variables are:

[0106]

[0107]

[0108] After transformation, we get:

[0109]

[0110] Single corrected element stresses extracted by the rainflow counting method The sensitivity expression to the design variables is:

[0111] where σ e，i (t) is the overall stress state of the e-unit at the time t of the i-th rainflow cycle, i.e., the time of maximum or minimum stress. For the two-dimensional case, σ e，i (t) contains three components, namely They are the stress component in the x direction, the stress component in the y direction, and the shear stress on the xy plane;

[0112]

[0113] Where L e A 0-1 sparse matrix is ​​selected so that U e =L e U, U is the displacement. The sensitivity of displacement U to the design variable is unknown, so this implementation method uses the adjoint method to solve it, and the expression is as follows: F is the global load,

[0114]

[0115] Update the expression for sensitivity to:

[0116] in and is the adjoint vector λ e The maximum and minimum values ​​of U max and U min are the maximum and minimum values ​​of displacement U.

[0117] Step 10: Repeat steps 3 to 9 to iterate the update of the isogeometric analysis design variables until the fatigue constrained topology optimization problem converges based on the change in the design variables.

[0118] Based on the obtained sensitivity information, the unconstrained MMA algorithm is used to update the design variables, and the augmented Lagrangian operator is updated according to the KKT condition until the approximate model converges, thereby obtaining the convergence of the isogeometric fatigue optimization problem at the current level.

[0119] Step 11: Preset the number of layers according to the multigrid method, use the continuation operator P to map the design variables on the coarse grid to the fine grid, update the design domain for solution, and repeat steps 3 to 10.

[0120] The data mapping and design domain relationship between multi-layer grids are as follows:

[0121] u h =Pu H ;n c =n / 2 c-1 , m c =m / 2 c-1 , c=1,……,L。

[0122] Step 12: Repeat steps 3 to 11 until the preset final level is reached and output the final solution result.

[0123] The present invention provides an efficient fatigue topology optimization method for isogeometric analysis. The independent linear decomposition method is used to improve the simplified calculation and counting problems of multi-cycle general fatigue load spectra, and the multi-grid conjugate gradient method is used to avoid the huge cost of solving the adjoint equation under large-scale fatigue problems. In addition, the augmented Lagrangian method and the multi-layer grid method are used to solve the local constraint problem of large-scale fatigue problems.

[0124] 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. An efficient fatigue topology optimization method for isogeometric analysis, characterized in that: The steps include: Step 1: Construct fatigue problem model; Step 2: Divide the number of multi-layer grid layers and the size of the design domain according to the design domain requirements, and use h insertion to generate the lowest level of isogeometric node vectors; Step 3: Assemble the interpolation operators and extension operators of each level of the multigrid method according to the isogeometric node vectors; Step 4: Decompose and store the independent load components in the load spectrum, and solve the combined weights of the independent load components; Step 5: Perform isogeometric analysis to solve each independent load component; Step 6: Calculate the stress spectrum of each unit of the independent load component, assemble it according to the weight to get the overall stress spectrum, and then perform the signed von-mises stress correction of the unit stress; Step 7: The peak and valley values ​​of each cycle and their indexes are extracted from the modified stress spectrum of each unit using the rain flow counting method, and the average stress and stress amplitude of the unit in each cycle are calculated, and then the Morrow correction is performed to obtain the equivalent stress amplitude; Step 8: Evaluate fatigue constraints of element fatigue damage based on the Palmgren-Miner linear cumulative damage model; Step 9: Approximate the fatigue constrained topology optimization problem according to the augmented Lagrangian method; Step 10: Repeat steps 3 to 9 to iterate the update of the isogeometric analysis design variables until the fatigue constrained topology optimization problem converges according to the change in the design variables; Step 11: According to the preset number of layers of the multigrid method, the design variables on the coarse grid are mapped to the fine grid using the continuation operator P, the design domain is updated for solution, and steps 3 to 10 are repeated; Step 12: Repeat steps 3 to 11 until the preset final level is reached and output the final solution result.

2. The efficient fatigue topology optimization method for isogeometric analysis according to claim 1 is characterized in that: The specific content of step 1 is: i is the numbered i-th design variable, and the design variable ρ i After distance smoothing, the smoothed design variables are obtained Then smooth the design variables Obtaining projection design variables through Heaviside projection The smoothing process expression is: W i，j = {r min -dist(i, j)},{j∈N i |dist(i,j)≤r min }, ρ j is related to the design variable ρ i The set of design variables within a certain Euclidean distance, N i is the total number of design variables in the design variable set, W i，j is the corresponding weight coefficient of the j-th design variable to the i-th design variable, r min is the filter radius, dist(i, j) is the design variable ρ i and ρ j The Euclidean distance between them; the expression of the Heaviside projection process is Where β and η are the steepness and threshold of the Heaviside projection respectively. The isogeometric analysis topology optimization model is established, and the corresponding mathematical expression is as follows: findρ=[ρ1,ρ2,…,ρ n ], subject to: g i (ρ)=D i -1<0, i=1,..., Ne, 0≤ρ e ≤1, e=1, ..., Ne, t=1,……NF,where ρ is the initial design variable vector, v e is the area or volume corresponding to the e-th design variable component, Ne is the total number of elements, NF is the total number of moments in the load spectrum, and D i is the damage variable of the ith unit, g i (ρ) is the fatigue constraint for the ith element, For projection-based design variables Assemble the obtained global stiffness matrix, U t and F t are the global displacement and load at the tth moment.

3. The efficient fatigue topology optimization method for isogeometric analysis according to claim 2 is characterized in that: The specific content of step 3 is: For the case where the topology optimization analysis domain is a single dimension, assume that there are two sets of node vectors and The subscripts i and j represent the i-th and j-th numbers of the node vectors. When used as symbols, they are a general reference to a set of sequences. The superscript H represents a coarse grid, and the superscript h represents a fine grid. Depend on After h insertion, the NURBS basis function obtained by two sets of node vectors and Satisfies the following linear transformation relationship: Where k is the repeatability of the B-spline, n h is the number of nodes in the grid, r = 1, 2, ..., k, r is the extension matrix The number of layers, where i and j are used as symbols to refer to the extension matrix With the intermediate interpolation matrix The rows and columns of , in specific calculations, represent the position of the matrix element in the i-th row and j-th column. The following solution relationship is satisfied for r: Obtain the restriction matrix R for transforming from fine grid to coarse grid, and use the transpose of the restriction matrix as the interpolation matrix P = R for transforming from coarse grid to fine grid T ; For the case where the topology optimization analysis domain is two-dimensional, for the coarse mesh node vector and The corresponding B-spline basis function is and The fine grid node vector obtained after conversion is and The corresponding B-spline basis function is and The restriction matrices for the conversion from fine grid to coarse grid in two dimensions are R 1，i，s and R 2，j，t , where i, s, j, and t represent the element positions in the two-dimensional matrix, and the two-dimensional transformation relationship is obtained: R a，b =R 1，i，s R 2，j，t , a=i+(j-1)n, b=s+(t-1)n′, n is R 1，i，s The number of coarse grid control points of the corresponding dimension, n′ is R 2，j，t The number of coarse grid control points of the corresponding dimension; R 1，i，s and R 2，j，t The transpose of is used as the interpolation matrix for transforming from coarse grid to fine grid; The data u corresponding to the coarse grid H The data u corresponding to the fine grid h The conversion expression is abbreviated as: H =Ru h ,u h =Pu H .

4. The efficient fatigue topology optimization method for isogeometric analysis according to claim 3 is characterized in that: The specific content of step 4 is: The decomposition load spectrum is as follows: F t =a t *F a +b t *F b +…+n t *F n , where F t is the load component at the tth load spectrum moment, F a 、F b , ..., F n is the independent load component obtained by decomposition, a t , b t ,...,n t It is the combined weight of each independent load component, and the combined weight is stored.

5. The efficient fatigue topology optimization method for isogeometric analysis according to claim 4, characterized in that: The specific content of step 6 is: for the independent load components solved in step 5, assemble them according to the stored combination weights. The assembled expression is as follows: a e，t , b e，t , ..., n e，t is the combined weight of unit e, is the independent load component that has been solved; then the stress spectrum assembled for each independent load component that has been solved is Perform a signed stress von-mises correction, specifically, a signed stress von-mises correction Δσ e (t) VM The expression is as follows: Δσ e (t) = σ e (t)-σ e (t0), where t0 is the loading start time; sign(·) is the sign function, which takes the value of 1 or -1; matrix W = [1 1 0] T ; V is the stress coefficient matrix, σ e (t) is the stress vector at element e, q is the stress penalty parameter, is the design variable ρ of unit e e The result after smoothing and Heaviside projection, B e and U e are the strain displacement matrix and displacement matrix of unit e respectively, D0 is the elastic matrix of solid material, E0 is Young's modulus of the solid material, and v is Poisson's ratio.

6. The efficient fatigue topology optimization method for isogeometric analysis according to claim 5, characterized in that: The specific content of step 7 is: For the reorganized and corrected stress spectrum, the expressions of the average stress and stress amplitude of unit e in the i-th cycle are: and are the maximum stress and minimum stress of unit e in the i-th cycle respectively; Morrow correction is then performed to obtain the equivalent stress amplitude The expression is σ′ f is the fatigue strength coefficient of the material.

7. The efficient fatigue topology optimization method for isogeometric analysis according to claim 6, characterized in that: The specific content of step 8 is: the unit fatigue damage D must satisfy the following fatigue constraints: Where N f，i is the number of fatigue life cycles corresponding to the stress of the i-th reorganized and corrected stress spectrum; n i is the cycle number corresponding to the stress of the ith reorganized and corrected stress spectrum; b is the fatigue strength index.

8. The efficient fatigue topology optimization method for isogeometric analysis according to claim 7, characterized in that: The specific content of step 9 is: according to the augmented Lagrangian method, the model of fatigue constrained topology optimization problem is approximated as follows: findρ=[ρ1,ρ2,…,ρ n ], subjectto: g i (ρ)=D i -1, 0≤ρ e ≤1, in, is the global stiffness matrix, which is abbreviated as K in the following. e and μ e is the Lagrange multiplier corresponding to the fatigue constraint of the e-th unit, and the update criterion of the Lagrange multiplier is: μ e，l+1 =max(μ e，l *1.1,10000),λ e，l+1 +μ e，l *g e (ρ), where λ e，l and μ e，l is the Lagrange multiplier at the lth iteration, μ e，l+1 and λ e，l+1 is the Lagrange multiplier at the l+1th iteration.

9. The efficient fatigue topology optimization method for isogeometric analysis according to claim 8, characterized in that: Integrate the fatigue constraint into the objective function and solve the sensitivity of the objective function to the design variables. The mathematical expression of the sensitivity is as follows: According to the chain rule, The sensitivity expressions for the design variables are: After transformation, we get Single corrected element stresses extracted by the rainflow counting method The expressions for the design variables are: where σ e，i (t) is the overall stress state of the e element at time t (maximum or minimum stress moment) of the i-th rainflow cycle. For the two-dimensional case, σ e，i (t) contains three components, namely They are the stress component in the x direction, the stress component in the y direction, and the shear stress on the xy plane; Where L e A 0-1 sparse matrix is ​​selected so that U e =L e U, U is the global displacement.

10. The efficient fatigue topology optimization method for isogeometric analysis according to claim 9, characterized in that: The sensitivity of the global displacement U to the design variables is solved by the adjoint method, and the expression is: F is the global load, Update the expression for sensitivity to: in and is the adjoint vector λ e The maximum and minimum values ​​of U max and U min are the maximum and minimum values ​​of displacement U.

Images