Efficient fatigue topology optimization method for isogeometric analysis

Abstract

The application provides an efficient fatigue topology optimization method for isogeometric analysis, and belongs to the technical field of geometric analysis, and mainly comprises the following steps: constructing a fatigue problem model, dividing multiple layers of grid methods according to the design domain requirement, and generating the lowest level isogeometric node vector by using h insertion; assembling each level interpolation operator and continuation operator of the multi-grid method according to the isogeometric node vector; decomposing and storing independent load components in a load spectrum, and solving the combination weight of each independent load component; assembling the overall stress spectrum according to the weight; evaluating the fatigue constraint of the unit fatigue damage; expressing the fatigue constraint topology optimization problem by using the augmented Lagrange method; performing the update iteration of the isogeometric analysis design variable until the fatigue constraint topology optimization problem converges according to the design variable change amount. The scheme improves the simplified calculation and counting problem of the multi-cycle general fatigue load spectrum, and reduces the solving cost.

Description

Technical Field

[0001] The present invention relates to the field of geometric analysis technology, 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 engineering optimization design continues to increase. Structural fatigue has gradually become a key factor in engineering practice. Fatigue refers to the phenomenon in which a structure gradually suffers damage and failure after being subjected to sufficient cyclic disturbances. Even if the stress level is below the yield strength of the material, long-term cyclic disturbances can cause the structure to fracture or fail. Fatigue failure is a widespread problem in the fields of machinery, aerospace, automotive, etc., and has a decisive impact on the safety and service life of structures.

[0003] Fatigue topology optimization is a technique that maximizes the fatigue resistance of a structure under cyclic perturbations by optimizing material distribution. Unlike static loads, variable perturbation 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 perturbations. 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 (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 the geometric model and the analysis model, 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 or even general loads that are closer to engineering practice, the computational complexity of 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 matrix solution cost 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 element of the independent load component, assemble it according to the weight to obtain the overall stress spectrum, and then perform the signed von-mises stress correction of the element stress;

[0013] Step 7: Use the rain flow counting method to extract the peak and valley values ​​and their indexes of each cycle from the modified stress spectrum of each element, and calculate the average stress and stress amplitude of the element in each cycle. Then perform Morrow correction to obtain the equivalent stress amplitude.

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

[0015] Step 9: Approximate the fatigue constrained topology optimization problem using 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 based on the change in the design variables.

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

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

[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, D i is the damage variable of the i-th unit, g i (ρ) is the fatigue constraint for the i-th 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 there are two sets of node vectors and The subscripts i and j represent the i-th and j-th numbers of the node vector, and are used as symbols to refer to a set of sequence symbols. 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 conversion 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 as symbols are used to refer to the extension matrix and the intermediate interpolation matrix The rows and columns of represent the position of the matrix element in the i-th row and j-th column in the specific calculation. The following solution relationship is satisfied for r:

[0021] Obtain the restriction matrix R for the transformation from fine grid to coarse grid, and use the transformation rank of the restriction matrix as the interpolation matrix P = R for the transformation 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 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: decomposing the load spectrum 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 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 signed stress von-mises correction, specifically, 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 the Young's modulus of the solid material, and v is the 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; the equivalent stress amplitude is obtained by Morrow correction 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, that is, 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. 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+1)th iteration.

[0034] Further preferably, the fatigue constraint is integrated into the objective function, and the sensitivity of the objective function to the design variable is solved, and the mathematical expression of the sensitivity is as follows:

[0035]

[0036] According to the chain rule, The expression of the sensitivity of the design variable is as follows: After transformation, the following is obtained

[0037] The single corrected unit stress σ e，i (t) is expressed as follows:

[0038] Where σ e，i (t) is the overall stress state of the e unit at the t moment (stress maximum moment or minimum moment) of the i rainflow cycle, for the two-dimensional case, σ e，i (t) contains three components, namely The stress component in the x direction, the stress component in the y direction and the shear stress in the xy plane respectively.

[0039]

[0040] Where L e is a 0-1 sparse selection matrix, so that U e = L e U, and U is the global displacement.

[0041] Further preferably, the sensitivity of the global displacement U to the design variable is solved by using the adjoint method, and the expression is as follows: F is the global load, The expression of the sensitivity is updated as follows: , Where and are the maximum value and the minimum value of the adjoint vector λ e ; U max and U min are the maximum value and the minimum value of the global displacement U.

[0042] The high-efficiency fatigue topology optimization method for isogeometric analysis provided by the application has the following beneficial effects compared with the prior art:

[0043] (1) The present invention applies a multi-layer grid method based on isogeometry, which utilizes the node refinement relationship of isogeometry to realize isogeometric density information mapping, so that the fine grid iteration starts directly from a more accurate intermediate optimization structure, thereby shortening the number of convergence iterations; it also adopts the method of independent linear component decomposition to solve the node displacement, avoiding the linear increase of the solution cost with the number of cycles, significantly reducing the redundant displacement calculation amount of the general load spectrum, and thus improving the efficiency of structural topology optimization;

[0044] (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 fatigue problems. While ensuring the accuracy of the optimization results, the amount of calculation is significantly reduced, and the computational efficiency of large-scale large condition number matrix inversion is greatly improved, making it possible to solve large-scale fatigue problems considering general loads. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0046] Figure 1 This is a flowchart of the steps of an efficient fatigue topology optimization method for isogeometric analysis according to the present invention;

[0047] Figure 2 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;

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

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

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

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

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

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

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

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

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

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

[0058] Specifically, the isogeometric analysis model based on the variable density method is established to optimize the design variables ρ 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.

[0059] Considering that the principle of isogeometric multigrid assembly will be introduced later, the assembly method of NURBS basis functions is given below. 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:

[0060] When p = 0,

[0061] When p>0,

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

[0063] NURBS basis function is an extension 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 follows: 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.

[0064] ρ i is the numbered i-th design variable, and the design variable ρ i After distance smoothing, we get Then the smoothed design variables Obtained by 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 the removal of 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, D i is the damage variable of the i-th 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.

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

[0066] Taking the two-dimensional design domain as an example, limited by the basic principle of h insertion, assuming that the number of layers in 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:

[0067] 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 mesh design domain size at level c and is a positive integer. For the calculated design domain size at the lowest level, perform isogeometric h-node insertion to obtain the node vector expressions in the ξ0 and η0 directions 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.

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

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

[0070] 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 vector, and are used as symbols to refer to a set of sequence symbols. 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 conversion 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 as symbols are used to refer to the extension matrix and the intermediate interpolation matrix The rows and columns of represent the position of the matrix element in the i-th row and j-th column in the specific calculation. The following solution relationship is satisfied for r: Obtain the restriction matrix R for the transformation from fine grid to coarse grid, and use the transpose of the restriction matrix as the interpolation matrix P = R for the transformation from coarse grid to fine grid T ;

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

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

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

[0074] 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 is the combined weight of each independent load component, and the combined weight is stored.

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

[0076] For the global displacement To solve the problem, interpolation and extension operators are used to assemble a multigrid preprocessing matrix. The first iteration reuses the previous solution to accelerate iterative convergence. Multigrid is a multi-level iterative method for solving linear equations. This method divides the computational domain into grids of different sizes and solves on the smallest, coarsest grid. The V-loop solution process of the multigrid conjugate gradient method with reused initial values ​​is as follows:

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

[0078] U 1 =Upre

[0079] for l=1:L-1

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

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

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

[0083] K l+1 =RK l P

[0084] U l+1 =0

[0085] end for

[0086] U l =K l \F l

[0087] for l=1:L-1

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

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

[0090] end for

[0091] Output: U 1

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

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

[0094] 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 signed stress von-mises correction, specifically, 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 the Young's modulus of the solid material, and v is the Poisson's ratio.

[0095] Step 7: Use the rain flow counting method to extract the peak and valley values ​​of each cycle and their indexes from the modified stress spectrum of each unit, and calculate the average stress and stress amplitude of the unit in each cycle, and then perform Morrow correction to obtain the equivalent stress amplitude.

[0096] 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 moment and the corresponding reversal extreme stress, Figure 6 The rain flow counting results are given in the form of a histogram, which shows the relevant parameters of the rain flow reversal obtained by statistics. Figure 7 The following is an illustration of the applied S-N curve range. Specifically, for the reorganized and corrected stress spectrum, the expressions for 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; the equivalent stress amplitude is obtained by Morrow correction The expression is σ′ f is the fatigue strength coefficient of the material.

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

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

[0099] Step 9: Approximate the fatigue constrained topology optimization problem using the augmented Lagrangian method.

[0100] The specific content is: According to the augmented Lagrangian method, the model of 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, that is, 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:

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

[0102] In the embodiment, the constraint function is integrated into the objective function by using the augmented Lagrange method, that is, the sensitivity of the objective function to the design variable needs to be solved, and the mathematical expression of the sensitivity is as follows:

[0103]

[0104] According to the chain rule, The sensitivity expression of the design variable is as follows:

[0105]

[0106]

[0107] After transformation, the following is obtained:

[0108]

[0109] For the single correction unit stress extracted by the rainflow counting method The sensitivity expression of the design variable is as follows:

[0110] Where σ e，i (t) is the overall stress state of the e unit at t moment of the i rainflow cycle, that is, the maximum stress moment or the minimum value moment, for the two-dimensional case, σ e，i (t) contains three components, that is, respectively, the stress component in the x direction, the stress component in the y direction, and the shear stress on the xy plane;

[0111]

[0112] Where L e is a 0-1 sparse selection matrix, so that U e = L e U, U is displacement. The sensitivity of the displacement U to the design variable is unknown, so the adjoint method is used to solve it, and the expression is as follows: F is the global load,

[0113]

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

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

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

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

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

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

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

[0121] This paper provides an efficient fatigue topology optimization method for isogeometric analysis. It employs an independent linear decomposition method to improve the simplified calculation and counting of multi-cycle general fatigue load spectra. It also employs a multigrid conjugate gradient method to avoid the significant cost of solving the adjoint equations for large-scale fatigue problems. Furthermore, it employs an augmented Lagrangian method and a multi-layer grid approach to solve local constraints in large-scale fatigue problems.

[0122] 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. 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 element of the independent load component, assemble it according to the weight to obtain the overall stress spectrum, and then perform the signed von-mises stress correction of the element stress; Step 7: Use the rain flow counting method to extract the peak and valley values ​​of each cycle and their indexes from the modified stress spectrum of each element, and calculate the average stress and stress amplitude of the element in each cycle, and then perform Morrow correction to obtain the equivalent stress amplitude; Step 8: Evaluate the fatigue constraints of the element fatigue damage based on the Palmgren-Miner linear cumulative damage model; Step 9: Approximate the fatigue constrained topology optimization problem using 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 based on the change in the design variables. 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; Step 12: Repeat steps 3 to 11 until the preset final level is reached and output the final solution.

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, D i is the damage variable of the i-th unit, g i (ρ) is the fatigue constraint for the i-th element, For projection-based design variables The global stiffness matrix obtained by assembling, 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 single-dimensional, 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 vector, and are used as symbols to refer to a set of sequence symbols. 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 conversion 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 and the intermediate interpolation matrix The rows and columns of represent the position of the matrix element in the i-th row and j-th column in the specific calculation. The following solution relationship is satisfied for r: Obtain the restriction matrix R for the transformation from fine grid to coarse grid, and use the transpose of the restriction matrix as the interpolation matrix P = R for the transformation 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 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: decompose the load spectrum 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 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 is 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 signed stress von-mises correction, specifically, 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 the Young's modulus of the solid material, and v is the 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; the equivalent stress amplitude is obtained by Morrow correction 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 i-th 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 the following expression: findρ=[ρ1,ρ2,…,ρ n ], subject to: g i (ρ)=D i -1, 0≤ρ e ≤1, in, is the global stiffness matrix, subsequently abbreviated as K, λ e and μ e is the Lagrange multiplier corresponding to the fatigue constraint of the e-th unit. 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 are the stress component in the x direction, the stress component in the y direction, and the shear stress on the xy plane; Among them L e For 0-1 sparse selection matrix, 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 using 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