An isogeometric robust topology optimization method for composite structures considering band constraints

Abstract

The application discloses a method for equal geometric robust topology optimization of a composite material structure considering band constraints, and comprises the following steps: firstly, constructing an equal geometric robust topology optimization model of the composite material structure considering material attribute correlation uncertainty and band constraints, calculating the mean value and standard deviation of the structural dynamic flexibility under the correlation uncertainty through Nataf transformation and Gauss-Hermite integral; then, solving the worst limit state approximation function corresponding to each limit state function and the uncertainty parameter point vector through a gradient search algorithm; finally, calculating the sensitivity of the objective function, the worst limit state approximation function and the volume constraint function, and obtaining the robust optimal topology structure satisfying the dynamic performance requirement by combining the moving asymptote method. The equal geometric robust topology optimization model of the composite material structure established by the application can truly reflect the correlation of material attribute uncertainty and the dynamic performance requirement, can efficiently obtain the optimization result, and has good engineering application value.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of structure optimization, and particularly relates to an equal geometric robust topology optimization method for a composite material structure considering frequency band constraints. BACKGROUND

[0002] In the fields of equipment manufacturing, automobiles, etc., the key structures of milling machines, engines and other equipment often bear a certain frequency of harmonic load during the working process. When the natural frequency of the key structures is close to the frequency of the harmonic load, resonance occurs, which may cause large deformation or fracture of the structures. Therefore, the topology optimization of such structures needs to consider the reliability of their dynamic performance.

[0003] Existing dynamic topology optimization researches are often directed at single-material structures, and less at composite material structures. Due to the complexity of the preparation process of composite materials, the material properties of the composite materials are uncertain, and there is a correlation between the properties of the same kind of raw materials. In addition, the current dynamic topology optimization researches on composite material structures mainly use structural dynamic flexibility, displacement and fundamental frequency as optimization performance indicators, and less use frequency band constraint methods. Therefore, there is an urgent need for a robust topology optimization method for composite material structures considering frequency band constraints and correlation uncertainty to meet the actual engineering needs. SUMMARY

[0004] To solve the problems of the prior art and achieve the purpose of avoiding the frequency of the harmonic load borne by the composite material structure and improving the stiffness, strength and vibration resistance of the structure, the present application adopts the following technical solutions:

[0005] The equal geometric robust topology optimization method for a composite material structure considering frequency band constraints comprises the following steps:

[0006] Step 1: An equal geometric analysis model of the composite material structure is constructed to obtain the density of control points.

[0007] Step 2: The related material uncertainty description of the composite material structure is obtained; a frequency band constraint is constructed according to the vibration resistance requirement of the composite material structure; the density of the control points is taken as a design variable, and the equal geometric robust topology optimization model for the composite material structure considering the frequency band constraint is constructed by combining the uncertainty description, the frequency band constraint, the dynamic stress constraint and the utilization rate constraint of the optimization space; and the minimum dynamic flexibility mean value and the dynamic flexibility standard deviation are obtained according to the displacement and the dynamic stiffness of the control points under the constraint conditions.

[0008] Step 3: The dynamic stiffness of the composite material structure under the action of the harmonic load is constructed.

[0009] Step 4: The global failure coefficient of the composite material structure considering the correlation uncertainty is calculated through the dynamic stress of the composite material structure at the control points.

[0010] Step 5: define a weighted objective function to calculate the mean and standard deviation of the structural dynamic flexibility under the influence of relevant uncertainties; at the same time, based on the frequency band constraint and the global failure coefficient, solve the uncertainty parameter point corresponding to each worst limit state function by a gradient search algorithm to calculate the corresponding material properties and worst limit state approximation function;

[0011] Step 6: update the control point density according to the objective function, the worst limit state approximation function, and the density sensitivity of the optimization space utilization constraint with respect to the control point, to obtain the optimal design variable.

[0012] Further, the composite structure isogeometric robust topology optimization model considering frequency band constraint in step 2 is constructed as follows:

[0013]

[0014] s.t.

[0015] g V (ρ)=V(ρ) / V0≤[V]

[0016]

[0017] P(g ω (ρ,X)≤exp(-0.5))≥P t

[0018] K d (ρ,X)U=F

[0019]

[0020] ρ={ρ i,j},10 -9 ≤ρ i,j ≤1 (i=1,2,…,m;j=1,2,...,n)

[0021] X=(V f ,E m ,ν m ,E f ,ν f ) T

[0022] wherein, represents the mean of the dynamic flexibility, represents the dynamic flexibility standard deviation, N e represents the number of units, u e (ρ,X) represents the displacement vector of unit e, e represents the unit corresponding to the control point, and the following is abbreviated as represents the dynamic stiffness matrix of unit e, and the following is abbreviated as g V (ρ) represents the optimization space utilization constraint function, V(ρ) represents the current optimization space size, V0 and [V] represent the design domain size and the upper limit of the space utilization constraint, respectively, and P(·) represents the probability calculation function. and [ω] TH [ ] represents the global failure coefficient of the structure constructed by the Tsai-Hill criterion and the P-mean condensation function, and the set upper limit value of the global failure coefficient, respectively. t G represents the given constraint reliability requirement. ω (ρ,X) represents the band-constrained performance function, exp(·) represents the exponential function, and K d (ρ,X) represents the global dynamic stiffness matrix, U represents the global displacement amplitude vector, F represents the global harmonic load amplitude vector, K(ρ,X) represents the global static stiffness matrix, and φ k Let ω represent the k-th structural mode. k (ρ,X) represents the k-th natural frequency, M(ρ) represents the overall mass matrix, and N ω ρ represents the pre-selected natural frequency order, and ρ represents the material density corresponding to the control point. i,j Indicates control point P i,j The corresponding material density, X represents a probability vector that follows a normal distribution, V f E f ν f E m ν m These represent the volume fraction of the reinforcing phase, Young's modulus, and Poisson's ratio of the composite material, respectively, and the Young's modulus and Poisson's ratio of the matrix material.

[0023] Furthermore, in step 2, the volume fraction V of the reinforcing phase in the composite material structure is considered. f Young's modulus E f Compared with Poisson's ratio ν f Young's modulus E of matrix material m Compared with the matrix Poisson ratio ν m The uncertainty and correlation between properties of the same material are described as a probability vector X = (V) following a normal distribution. f E m ,ν m E f ,ν f ) T ,and in and These represent the mean values ​​of Young's modulus and Poisson's ratio of the matrix material, respectively. and These represent the standard deviations of Young's modulus and Poisson's ratio of the matrix material, respectively. denotes the correlation coefficient of Young's modulus and Poisson's ratio of the matrix material, denotes the mean value of Young's modulus and Poisson's ratio of the reinforcing phase, denotes the standard deviation of Young's modulus and Poisson's ratio of the reinforcing phase, denotes the correlation coefficient of Young's modulus and Poisson's ratio of the reinforcing phase. denotes the correlation coefficient of Young's modulus and Poisson's ratio of the reinforcing phase. denotes the correlation coefficient of Young's modulus and Poisson's ratio of the reinforcing phase.

[0024] Further, in step 2, according to the vibration resistance requirement of the composite structure, a frequency band constraint function based on the smooth Heaviside function and the P-mean condensation function is constructed, specifically:

[0025]

[0026] wherein g ω (ρ,X) denotes the frequency band constraint performance function, ρ denotes the material density corresponding to the control point, ω k (ρ,X) denotes the k(k=1, 2, …, N ω )th order natural frequency, which is abbreviated as ω N ω denotes the pre-selected natural frequency order, p b denotes the coefficient of the P-mean condensation function corresponding to the frequency band constraint, denotes the smooth Heaviside function of the kth order natural frequency with respect to the constraint frequency band, specifically:

[0027]

[0028] wherein ω c denotes the midpoint of the working frequency band determined according to the simple harmonic excitation, δ denotes the coefficient of the frequency band constraint performance function, ω w denotes half of the width of the limited frequency band region, specifically:

[0029] ω w =(ω upp -ω low ) / 2

[0030] wherein ω low denotes the lower limit value of the frequency band constraint; ω upp denotes the upper limit value of the frequency band constraint.

[0031] Further, the step 3 includes the following steps:

[0032] Step 3.1: calculating the element static stiffness matrix k e (ρ,X) and the element mass matrix m e (ρ) of the composite structure based on the plate model;

[0033] Step 3.2: Calculate the overall stiffness matrix K(ρ, X) of the structure, the overall mass matrix M(ρ) of the structure and the overall damping matrix C(ρ, X) of the structure by the unit stiffness matrix and the unit mass matrix based on Rayleigh damping method;

[0034] Step 3.3: Calculate the overall dynamic stiffness matrix K(ρ, X) of the composite structure based on the overall stiffness matrix K(ρ, X), the overall mass matrix M(ρ) of the structure and the overall damping matrix C(ρ, X). d (ρ, X).

[0035] Further, the isogeometric analysis model in the step 1 is a non-uniform rational B-spline (NURBS) mesh model, the mesh model is entirely Gauss subdivided, and the density of the i-th row and j-th column element Gauss point in the e-th element is obtained from the corresponding control point of the element e.

[0036] Further, the step 4 comprises the following steps:

[0037] Step 4.1: Calculate the dynamic stress of the composite structure at the element Gauss point

[0038] Step 4.2: Calculate the global failure coefficient of the composite structure based on Tsai-Hill failure criterion and P-mean condensation function. Specifically,

[0039]

[0040] Wherein, ρ represents the material density corresponding to the control point, X represents the probability vector satisfying the normal distribution, represents the failure coefficient of the element Gauss point considering the related uncertainty, N e , N G respectively represent the number of elements of the NURBS mesh model, the row or column number of the internal Gauss point of the element, and specifically,

[0041]

[0042] Wherein, represents the failure coefficient of the element Gauss point considering the related uncertainty, and the following is abbreviated as represents the strength coefficient matrix of the material, and the following is abbreviated as p t represents the P-mean condensation function coefficient of the stress constraint function.

[0043] ​​Further, in the step 5, the mean value and standard deviation of the structural dynamic flexibility under the influence of the related uncertainty are calculated based on the Nataf transformation, and specifically comprising the following steps:

[0044] Step 5.1.1: Introducing weight coefficients w X Define the weighted objective function

[0045]

[0046] Step 5.1.2: Calculate the structural dynamic flexibility c d The approximate value of the origin moment of (p, X) is specifically as follows:

[0047] Step 5.1.2.1: Based on the Nataf transformation, convert the related uncertainty vector X into an independent standard normal random vector

[0048] Z = (Z1, Z2, Z3, Z4, Z5) T ;

[0049] Step 5.1.2.2: Based on the single variable dimension reduction method, use the independent standard normal random vector Z to rewrite the origin moment of the structural dynamic flexibility

[0050] Specifically, it is as follows:

[0051]

[0052] Wherein, m represents the order of the origin moment of the target performance to be solved, represents the marginal probability density distribution function of the independent standard normal random variable Z t ; c N (p, Z) represents the composite material structure dynamic flexibility represented by the independent standard normal random vector Z;

[0053] Step 5.1.2.3: Based on Gauss-Hermite integral, calculate the approximate value of the origin moment when m = 1, 2, and specifically as follows:

[0054]

[0055] Wherein, H represents the number of selected Gauss-Hermite integral points; w h represents the weight corresponding to the hth Gauss-Hermite integral point; Wherein is the corresponding value of the hth Gauss-Hermite integral point; represents the nominal value of the structural dynamic flexibility obtained by taking the mean value of each variable in the vector Z.​​​

[0056] Step 5.1.2.4: Based on the Nataf inverse transformation, the vector obtained by Gauss-Hermite integration point valuation Determine the corresponding material uncertainty vector Use denotes the origin matrix approximation value when m = 1, 2;

[0057] Step 5.1.3: Calculate the mean and standard deviation of the structural dynamic flexibility under the influence of the relevant uncertainty X, specifically:

[0058]

[0059] where, denotes the mean of the dynamic flexibility, denotes the standard deviation of the dynamic flexibility.

[0060] Further, in the step 5, the uncertainty parameter point corresponding to each worst limit state function is solved by a gradient search algorithm, specifically including the following steps:

[0061] Step 5.2.1: Based on the frequency band constraint performance function g ω (ρ, X) and the stress constraint performance function Construct the limit state function, denoted as:

[0062]

[0063] where, ρ represents the material density corresponding to the control point, X represents the probability vector satisfying the normal distribution, exp(·) represents the exponential function, denotes the limit state function corresponding to the stress constraint performance function, denotes the limit state function corresponding to the frequency band constraint performance function;

[0064] Step 5.2.2: Rewrite the limit state function into the form represented by an independent standard normal vector Z and

[0065] Step 5.2.3: Establish a reliability model respectively to obtain the uncertainty parameter point of the limit state function on the reliability index surface, specifically:

[0066]

[0067] s.t. ||Z|| = β t

[0068]

[0069] s.t. ||Z|| = βt

[0070] where β t represents the reliability index corresponding to the constraint performance function, and P t is the corresponding reliability requirement; t = Φ -1 (P t ) ; Φ -1 (·) represents the inverse function of the standard normal cumulative distribution function;

[0071] Step 5.2.4: Calculate the stress limit state function and the frequency band limit state function about the sensitivity of each standard normal variable Z t (t = 1, 2,..., 5); by the gradient search algorithm, the uncertainty parameter point vector of each worst limit state function on the reliability index surface is calculated and and the material property vector corresponding to the above uncertainty parameter points is obtained by Nataf inverse transformation and and the worst limit state function and

[0072] Step 5.2.5: Based on each worst limit state function and establish the limit state approximation function in normalized form Specifically:

[0073]

[0074] Based on Nataf transformation, the independent standard normal random vector Z is used to replace the related uncertainty vector X to establish the reliability analysis model, which is specifically:

[0075]

[0076] s.t. ||Z|| = β t

[0077] Step 5.2.6: Solve the function about the sensitivity of variable Z t , based on the above sensitivity solving results combined with the gradient search algorithm, the uncertainty parameter point approximation value vector Z ap on the reliability index surface corresponding to Z when taking the minimum value is calculated ap , and the material property vector X and the worst limit state approximation function and

[0078] Further, in step 6, the target function is calculated Each worst limit state approximation function With And the optimization space utilization rate constraint g V (ρ) Sensitivity of each unit control point density ρ e,a,b According to the control point density sensitivity calculation result, the control point density is updated combined with the moving asymptote method, and the convergence condition is checked, that is, the average of the relative change degree of the target function In the continuous number of iteration steps is not higher than the threshold value, if the convergence condition is not met, return to step 3 iteration.

[0079] The advantages and beneficial effects of the present application are:

[0080] Compared with the equal geometry robust topology optimization method without frequency band constraint, the present application considers the influence of harmonic load, and the natural frequency of the obtained optimization structure can avoid the frequency of the received harmonic load, and has good stiffness, strength and vibration resistance; The present application describes the material uncertainty as a random variable satisfying the binary normal distribution, uses the correlation coefficient to describe the correlation between the uncertainties, and proposes a target performance statistical value solving method based on Nataf transformation, converts the correlated variables into independent standard normal variables, and then uses the single variable dimension reduction method combined with Guass-Hermite transformation to solve the target performance statistical moment. Compared with the method without considering correlation, the target performance statistical value can be more accurately evaluated; The present application proposes a multi-constraint performance reliability evaluation method considering the worst constraint performance approximation, uses the worst limit state function obtained by each probability constraint to construct a normalized reliability analysis model, and obtains the worst limit state approximation function, so that each probability constraint is in the same optimization background, and compared with the method without approximation processing, a more optimal structure target performance can be obtained. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 The method flowchart of the embodiment of the present application is shown in the figure.

[0082] Figure 2a The schematic diagram of the engine connecting rod structure design domain in the embodiment of the present application is shown in the figure.

[0083] Figure 2b The schematic diagram of the initial design of the engine connecting rod structure in the embodiment of the present application is shown in the figure.

[0084] Figure 3 The optimal topology diagram of the engine connecting rod structure in the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0085] The specific embodiments of the present application are described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely intended to illustrate and explain the present application, and are not intended to limit the present application.

[0086] As shown in Figure 1 , the equal geometric robust topology optimization method of composite structure considering frequency band constraint, considering the influence of simple harmonic load and related uncertainty on the composite structure under actual working condition, the sample sufficient material property uncertainty is described as variable satisfying binary normal distribution, the correlation coefficient is used to describe the correlation between different material properties of the same material, and the dynamic flexibility mean value and standard deviation are used as the objective function, the equal geometric robust topology optimization model of composite structure under dynamic stress constraint, frequency band constraint and material consumption constraint is established, which fully reflects the requirements of high stiffness, high strength, high vibration resistance and lightweight design of composite structure. At the same time, a reliability evaluation method considering worst constraint performance approximation is proposed, which is used to obtain the unified worst constraint performance approximation point. On this basis, the sensitivity of the objective function, the limit state approximation function and the volume constraint function is derived and solved, and the optimal topology design scheme is obtained by using the moving asymptote method. The optimization method includes the following steps:

[0087] Step 1: Construct the Non-Uniform Rational B-Splines (NURBS) equal geometric analysis model of the composite structure, and perform Gaussian subdivision on the NURBS grid model as a whole. The density of the i-th element in the ii-th row and jj-th column element Gaussian point is obtained by the control point corresponding to the element e

[0088] As shown in Figure 2a 、 Figure 2b , the design information in the embodiment of the present application is the actual application data in the topology optimization design of the connecting rod in a certain complex component manufactured by using TM800s / M21 carbon fiber reinforced composite material. The left end of the connecting rod is fixed, and the right end is subjected to three horizontal left simple harmonic loads F1, F2 and F3 with a frequency of 80Hz. The load amplitude f1=f2=f3=300N; Figure 2b The upper base length of the design domain is 27, the lower base length is 36, and the height is 130mm.

[0089] On the basis of constructing the NURBS curve, two direction node vectors u=(u0, u1,..., u 101+2+1 ) and v=(v0, v1,..., v 21+2+1 ) are introduced, 101×21 control points P i,j ,i=1,2,...,101; j=1,2,...,21 and 101×21 weights w​i,j i = 1,2,...,101; j = 1,2,...,21; the overall NURBS mesh model is subjected to Gaussian subdivision, and the ith row and jj column unit Gaussian points in the e unit are obtained from the control points corresponding to the unit e Density

[0090] Step 2: Construct an isogeometric robust topology optimization model of the composite structure considering the correlation of uncertainties and band constraints, specifically comprising the following steps:

[0091] Step 2.1: Consider the volume fraction V of the reinforcing phase of the composite structure f , Young's modulus E f and Poisson's ratio v f , the Young's modulus E of the matrix material m and the Poisson's ratio v of the matrix m uncertainties and the correlation between the properties of the same material are described as a probability vector X = (V f , E m , v m , E f , v f ) that satisfies the normal distribution T , and wherein and represent the mean values of the Young's modulus and Poisson's ratio of the matrix material, and represent the standard deviations of the Young's modulus and Poisson's ratio of the matrix material, represents the correlation coefficient of the Young's modulus and Poisson's ratio of the matrix material; and represent the mean values of the Young's modulus and Poisson's ratio of the reinforcing phase, and represent the standard deviations of the Young's modulus and Poisson's ratio of the reinforcing phase, represents the correlation coefficient of the Young's modulus and Poisson's ratio of the reinforcing phase.

[0092] In the embodiment of the present application, the uncertainties of the reinforcing fiber volume fraction V f , Young's modulus E f and Poisson's ratio v f , the Young's modulus E of the matrix material m and the Poisson's ratio v of the matrix m and the correlation between the Poisson's ratio and the Young's modulus of the same material are described as V f ~ N(0.6, 0.01 2 ), (E m , v m ) ~ N(3500, 0.35, 1502 0.01 2 0.8) and (E f , v f ) ~ N(240000, 0.30, 4800 2 0.01 2 0.8).

[0093] Step 2.2: According to the vibration resistance requirements of the composite structure, a frequency band constraint function based on the smooth Heaviside function and the P-mean condensation function is constructed, which is specifically:

[0094]

[0095] where g ω (p, X) represents the frequency band constraint performance function, p represents the material density corresponding to the control point, ω k (p, X) represents the k(k = 1, 2, …, N ω ) order natural frequency, which is abbreviated as ω N ω represents the pre-selected natural frequency order; p b represents the coefficient of the P-mean function corresponding to the frequency band constraint; represents the smooth Heaviside function of the k order natural frequency with respect to the constraint frequency band, which is specifically:

[0096]

[0097] where ω c represents the midpoint of the working frequency band determined according to the simple harmonic excitation, δ represents the coefficient of the frequency band constraint performance function, ω w represents half of the width of the limited frequency band region, which is specifically:

[0098] ω w = (ω upp - ω low ) / 2 (3)

[0099] where ω low represents the lower limit value of the frequency band constraint; ω upp represents the upper limit value of the frequency band constraint.

[0100] Step 2.3: Taking the control point density as the design variable, a composite structure isogeometric robust topology optimization model considering frequency band constraint is constructed:

[0101]

[0102] where represents the mean value of dynamic flexibility, represents the standard deviation of dynamic flexibility, N eNumber of elements, u e (ρ, X) represents the displacement vector of element e, which is abbreviated as (ρ, X) represents the dynamic stiffness matrix of element e, which is abbreviated as g V (ρ) represents the optimization space utilization constraint function, V(ρ) represents the current optimization space size, V0 and [V] represent the design domain size and the upper limit value of the space utilization constraint respectively, P(·) represents the probability calculation function, and [ω TH ] respectively represent the structure global failure coefficient composed of Tsai-Hill criterion and P-mean condensation function, and the set upper limit value of global failure coefficient, P t represents the given constraint reliability requirement, K d (ρ, X) represents the overall dynamic stiffness matrix, U represents the overall displacement amplitude vector, F represents the overall harmonic load amplitude vector, K(ρ, X) represents the overall static stiffness matrix, φ k represents the kth order structure modal, M(ρ) represents the overall mass matrix, N ω represents the pre-selected natural frequency order, ρ i,j is the control point P i,j corresponding to the material density.

[0103] In the embodiment of the application, and represent the mean and standard deviation of the dynamic flexibility performance of the connecting rod in a certain complex component; K d (ρ, X) represents the overall dynamic stiffness matrix of the connecting rod; U represents the overall displacement amplitude vector of the connecting rod; F represents the harmonic load amplitude vector acting on the connecting rod; K(ρ, X) represents the overall static stiffness matrix of the connecting rod; φ k represents the kth order modal of the connecting rod; M(ρ) represents the overall mass matrix of the connecting rod.

[0104] Step 3: Construct the dynamic stiffness matrix of the composite structure (connecting rod) under the action of the harmonic load, specifically:

[0105] Step 3.1: Calculate the element static stiffness matrix k e (ρ, X) and the element mass matrix m e (ρ) of the composite structure (connecting rod) based on the plate model;

[0106] Step 3.2: Based on the Rayleigh damping method, calculate the overall stiffness matrix K(ρ, X), the overall mass matrix M(ρ) and the overall damping matrix C(ρ, X) of the structure through the element stiffness matrix and the element mass matrix;

[0107] Step 3.3: Calculate the overall dynamic stiffness matrix K of the composite structure (link) based on K(p, X), M(p) and C(p, X) d (p, X).

[0108] Step 4: Calculate the global failure coefficient of the composite structure (link) considering the related uncertainty, specifically:

[0109] Step 4.1: Calculate the element Gauss point The dynamic stress of the composite structure (link) at the element Gauss point

[0110] Step 4.2: Calculate the global failure coefficient of the composite structure (link) based on the Tsai-Hill failure criterion and the P-mean function Specifically:

[0111]

[0112] wherein, represents the failure coefficient at the element Gauss point considering the related uncertainty, N e = 2000, N G = 3 respectively represent the number of elements of the NURBS mesh model, the row or column number of the element Gauss point, specifically:

[0113]

[0114] wherein, represents the failure coefficient at the element Gauss point considering the related uncertainty, hereinafter abbreviated as represents the strength coefficient matrix of the material, hereinafter abbreviated as p t = 12 represents the P-mean condensation function coefficient of the stress constraint function.

[0115] Step 5: Calculate the mean and standard deviation of the structure dynamic flexibility under the influence of the related uncertainty based on the Nataf transformation, and solve the uncertainty parameter point vector corresponding to each worst limit state function through the gradient search algorithm;

[0116] wherein, the mean and standard deviation of the structure dynamic flexibility under the influence of the related uncertainty are calculated based on the Nataf transformation, specifically:

[0117] Step 5.1.1: Introduce the weight coefficient w X Define the weighted objective function

[0118]

[0119] Step 5.1.2: Calculate the structure dynamic flexibility c at m = 1, 2 d The approximate value of the origin moment of (p, X) is Specifically, it is

[0120] Step 5.1.2.1: Based on the Nataf transformation, convert the related uncertainty vector X into an independent standard normal random vector Z = (Z1, Z2, Z3, Z4, Z5) T ;

[0121] Step 5.1.2.2: Based on the single variable dimension reduction method, rewrite the origin moment of the structure dynamic flexibility using the independent standard normal random vector Z is Specifically, it is

[0122]

[0123] Where m represents the order of the target performance origin moment to be solved, represents the marginal probability density distribution function of the independent standard normal random variable Z t ; c N (p, Z) represents the dynamic flexibility of the composite structure (connecting rod) represented by the independent standard normal random vector Z;

[0124] Step 5.1.2.3: Calculate the approximate value of the origin moment at m = 1, 2 based on Gauss-Hermite integration, specifically:

[0125]

[0126] Where H represents the number of selected Gauss-Hermite integration points; w h represents the weight corresponding to the hth Gauss-Hermite integration point; Where is the corresponding value of the hth Gauss-Hermite integration point; represents the nominal value of the structure dynamic flexibility when each variable in the vector Z takes the mean value;

[0127] Step 5.1.2.4: Based on the Nataf inverse transformation, the vector determined by the Gauss-Hermite integration point value corresponds to the material uncertainty vector Using represents the approximate value of the origin moment at m = 1, 2;

[0128] Step 5.1.3: Calculate the mean and standard deviation of the structure dynamic flexibility under the influence of the related uncertainty X, specifically:

[0129]

[0130] The uncertainty parameter point vector corresponding to each worst limit state function is solved by a gradient search algorithm, specifically:

[0131] Step 5.2.1: Based on the frequency band constraint performance function g ω (ρ, X) and the stress constraint performance function The limit state function is constructed and expressed as:

[0132]

[0133]

[0134] Step 5.2.2: Rewrite the limit state function into a form represented by an independent standard normal vector Z and

[0135] Step 5.2.3: Establish a reliability model respectively to obtain the uncertainty parameter points of the limit state function on the reliability index surface, specifically:

[0136]

[0137] where β t represents the reliability index corresponding to the constraint performance function, and its corresponding reliability requirement P t is calculated, β t = Φ -1 (P t ); Φ -1 (·) represents the inverse function of the standard normal cumulative distribution function;

[0138] Step 5.2.4: Calculate the stress limit state function and the frequency band limit state function The sensitivity of each standard normal variable Z t (t = 1, 2,..., 5) is calculated. Through a gradient search algorithm, the uncertainty parameter point (Minimum Performance Target Point, MPTP) vector of each worst limit state function on the reliability index surface is calculated and and the material property vector corresponding to the above uncertainty parameter points is obtained using Nataf inverse transformation and and the worst limit state functions and

[0139] Step 5.2.5: Based on each worst limit state function and The limit state approximation function in normalized form is established, specifically as follows:

[0140]

[0141] Based on the Nataf transformation, the independent standard normal random vector Z is used to replace the correlated uncertainty vector X, and a reliability analysis model is established, specifically as follows:

[0142]

[0143] Step 5.2.6: solving the function The sensitivity of the variable Z t is calculated based on the above sensitivity solving result combined with the gradient search algorithm, and the approximation value vector Z of the uncertainty parameter point on the reliability index surface corresponding to Z when the minimum value is taken ap The corresponding material attribute vector X ap and the worst limit state approximation function are calculated through the Nataf inverse transformation

[0144] Step 6: calculating the objective function The worst limit state approximation function and and the volume constraint function g V (ρ) are calculated based on the sensitivity of each unit control point density ρ e,a,b ; according to the sensitivity calculation results of the above objective function, the worst limit state approximation function and the volume constraint function with respect to the control point density, the control point density is updated combined with the moving asymptote method; the convergence condition (the average of the relative change degree of the objective function in continuous number of iteration steps is not higher than the threshold value) is checked, and if the convergence condition is not met, the step 3 iteration is returned.

[0145] In the embodiment of the application, the convergence condition is that the average of the relative change degree of the objective function in continuous 3 iteration steps is lower than 0.05; after the 30th iteration, the convergence condition is met, the average of the structure dynamic flexibility after optimization and the standard deviation global failure coefficient meet the design requirements and working requirements of the connecting rod in the complex component, thereby verifying the effectiveness of the method of the application, and the corresponding engine connecting rod structure optimal topology is shown as Figure 3 .

[0146] The above examples are only used to illustrate the technical solutions of the present application, and are not intended to limit the present application; although the present application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that the technical solutions recorded in the foregoing examples can be modified, or some or all of the technical features can be replaced by equivalents; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. An isogeometric robust topology optimization method for composite structures considering frequency band constraints, characterized in that The method comprises the following steps: Step 1: constructing an isogeometric analysis model of the composite structure to obtain the density of control points; Step 2: obtaining a description of related material uncertainties of the composite structure; constructing a frequency band constraint according to the vibration resistance requirement of the composite structure; taking the density of the control points as design variables, combining the uncertainty description, the frequency band constraint, a dynamic stress constraint and an optimization space utilization constraint, constructing an isogeometric robust topology optimization model of the composite structure considering the frequency band constraint, and obtaining the minimum mean and standard deviation of dynamic flexibility under the constraint condition according to the displacement and dynamic stiffness of the control points; Step 3: constructing the dynamic stiffness of the composite structure under the action of a simple harmonic load; Step 4: calculating the global failure coefficient of the composite structure considering the related uncertainties through the dynamic stress of the composite structure at the control points; Step 5: defining a weighted objective function to calculate the mean and standard deviation of the dynamic flexibility of the structure under the influence of the related uncertainties; simultaneously, based on the frequency band constraint and the global failure coefficient, solving the uncertainty parameter points corresponding to each worst limit state function by a gradient search algorithm to calculate the corresponding material properties and worst limit state approximation function; Step 6: updating the control point density according to the sensitivity of the objective function, the worst limit state approximation function and the optimization space utilization constraint with respect to the density of the control points to obtain the optimal design variables.

2. The method of isogemetric robust topology optimization of composite structures considering band constraints according to claim 1, characterized in that: The isogeometric robust topology optimization model of the composite structure considering the frequency band constraint in step 2 is as follows: s.t. g V (p) = V(p) / V0≤ [V] P(g ω (ρ,X)≤exp(-0.5))≥P t K d (ρ,X)U = F p = {p i,j}, 10 -9 ≤ p i,j ≤ 1 i = 1, 2,..., m; j = 1, 2,..., n X = (V f ,E m ,ν m ,E f ,ν f ) T where, denotes the mean of dynamic compliance, denotes the standard deviation of dynamic compliance, N e denotes the number of elements, u e denotes the displacement vector of element e, e denotes the control point corresponding to element, denotes the dynamic stiffness matrix of element e, g V denotes the optimization space utilization constraint function, V(ρ) denotes the current optimization space size, V0 and [V] denote the design domain size and the upper limit of space utilization constraint, respectively, P(·) denotes the probability calculation function, and [ω TH ] denote the structure global failure coefficient and the upper limit of the set global failure coefficient, respectively, P t denotes the given constraint reliability requirement, g ω denotes the frequency band constraint performance function, exp(·) denotes the exponential function, K d denotes the overall dynamic stiffness matrix, U denotes the overall displacement amplitude vector, F denotes the overall harmonic load amplitude vector, K(ρ, X) denotes the overall static stiffness matrix, φ k denotes the kth order structure mode, ω k denotes the kth order natural frequency, M(ρ) denotes the overall mass matrix, N ω denotes the pre-selected natural frequency order, ρ denotes the material density corresponding to the control point, ρ i,j denotes the control point P i,j corresponding material density, X denotes the probability vector satisfying the normal distribution, V f , E f , ν f , E m , ν m denote the volume fraction of the reinforcing phase, Young's modulus and Poisson's ratio of the composite structure, and the Young's modulus and Poisson's ratio of the matrix material, respectively.

3. The method of isogemetric robust topology optimization of composite structures considering band constraints according to claim 2, characterized in that: In step 2, the volume fraction V of the reinforcing phase of the composite structure is considered f , Young's modulus E f , and Poisson's ratio v f , Young's modulus E m of the matrix material, and Poisson's ratio v m of the matrix material, the uncertainties and correlations between the properties of the same material are described as a probability vector X = (V f , E m , v m , E f , v f ) that satisfies a normal distribution T ; and wherein and E and v represent the mean value of the Young's modulus and the Poisson's ratio of the matrix material, respectively, and and v represent the standard deviation of the Young's modulus and the Poisson's ratio of the matrix material, respectively, r represents the correlation coefficient of the Young's modulus and the Poisson's ratio of the matrix material; and E and v represent the mean value of the Young's modulus and the Poisson's ratio of the reinforcing phase, respectively, and and v represent the standard deviation of the Young's modulus and the Poisson's ratio of the reinforcing phase, respectively, r represents the correlation coefficient of the Young's modulus and the Poisson's ratio of the reinforcing phase.

4. The method of isogemetric robust topology optimization of composite structures considering band constraints according to claim 2, characterized in that: In step 2, a frequency band constraint function based on a smoothing function and a condensation function is constructed according to the vibration resistance requirement of the composite structure, and specifically is as follows: where g ω (ρ,X) represents the frequency band constraint performance function, ρ represents the material density corresponding to the control point, ω k (ρ,X) represents the k-th = 1,2,…,N ω natural frequency, hereinafter abbreviated as N ω represents the pre-selected natural frequency order, p b represents the coefficient of the condensation function corresponding to the frequency band constraint, represents the smoothing function of the k-th natural frequency with respect to the constraint frequency band, specifically: where ω c represents the midpoint of the operating frequency band determined according to the simple harmonic excitation, δ represents a coefficient of the frequency band constraint performance function, ω w represents half of the width of the limited frequency band region, specifically: ω w = (ω upp - ω low ) / 2 where ω low represents a lower limit value of the frequency band constraint; and ω upp represents an upper limit value of the frequency band constraint.

5. The method of isogemetric robust topology optimization of composite structures considering band constraints of claim 1, wherein: Step 3 comprises the following steps: Step 3.1: calculating the element static stiffness matrix and the element mass matrix of the composite structure based on a plate model; Step 3.2: calculating the overall stiffness matrix, the overall mass matrix and the overall damping matrix of the structure based on the Rayleigh damping method through the element stiffness matrix and the element mass matrix; Step 3.3: calculating the overall dynamic stiffness matrix of the composite structure based on the overall stiffness matrix, the overall mass matrix and the overall damping matrix.

6. The method of isogemetric robust topology optimization of composite structures considering band constraints of claim 1, wherein: The isogeometric analysis model in the step 1 is a non-uniform rational B-spline grid model, and the whole grid model is subjected to Gauss subdivision, so that the corresponding control points of the element e are used to obtain the element Gauss points in the i-th row and j-th column in the e-th element density 7. The method of isogemetric robust topology optimization of composite structures considering band constraints according to claim 6, characterized in that: Step 4 comprises the following steps: Step 4.1: Compute unit Gaussian points Dynamic stress at composite structures Step 4.2: Calculate the global failure coefficient of the composite structure based on the failure criterion and the cohesion function Specifically: wherein p represents the material density corresponding to the control point, X represents the probability vector satisfying the normal distribution, represents the failure coefficient at the unit Gaussian point considering the related uncertainty, N e , N G respectively represent the number of units, the row or column number of the Gaussian points inside the unit, and specifically: wherein, represents a failure coefficient at a unit Gaussian point considering the correlation uncertainty; represents a strength coefficient matrix of the material; p t represents a condensation function coefficient of the stress constraint function.

8. The method of isogemetric robust topology optimization of composite structures considering band constraints of claim 1, wherein: In step 5, the mean and standard deviation of the dynamic flexibility of the structure under the influence of the related uncertainties are calculated based on Nataf transformation, and specifically comprise the following steps: Step 5.1.1 : Introducing the weight coefficient w X Defining the weighted objective function Step 5.1.2: Calculate the structure dynamic flexibility c for m = 1, 2 d The approximation of (p, X) origin matrix is as follows: , and the approximation of is as follows: Step 5.1.2.1 : Transform the correlated uncertainty vector X into an independent standard normal random vector Z = (Z1, Z2, Z3, Z4, Z5) based on the Nataf transformation T ; Step 5.1.2.2: Rewrite the origin point of the structural dynamic flexibility based on the univariate dimension reduction method using independent standard normal random vector Z For Specifically: where m denotes the order of the sought target performance origin moment, denotes the marginal probability density distribution function of the independent standard normal random variable Z t c N (ρ, Z) denotes the dynamic compliance of the composite structure represented by the independent standard normal random vector Z Step 5.1.2.3: calculating the origin matrix approximation value when m = 1, 2 based on integration, and specifically is as follows: wherein H represents the number of selected integration points; w h represents the weight corresponding to the hth integration point wherein is the corresponding value for the hth integration point; denotes the nominal value of the structural dynamic flexibility obtained by taking the average of each variable in the vector Z. Step 5.1.2.4: Vector obtained by evaluating the integral points based on the Nataf inverse transform determining a corresponding material uncertainty vector using denotes the zeroth order approximation for m = 1,2; Step 5.1.3: calculating the mean and standard deviation of the dynamic flexibility of the structure under the influence of the related uncertainty X, and specifically is as follows: wherein, represents the mean of the dynamic compliance, represents the standard deviation of the dynamic compliance.

9. The method of isogemetric robust topology optimization of composite structures considering band constraints of claim 1, wherein: In step 5, the uncertainty parameter points corresponding to each worst limit state function are solved by a gradient search algorithm, and specifically comprise the following steps: Step 5.2.1 : Band-based constraint performance function g ω ( p, X ) and stress constraint performance function Constructing the limit state function, denoted as: wherein p represents the material density corresponding to the control point, X represents a probability vector satisfying a normal distribution, exp(·) represents an exponential function, represents a limit state function corresponding to the stress constraint performance function, represents a limit state function corresponding to the frequency band constraint performance function; Step 5.2.2: Rewrite the limit state function in terms of the independent standard normal vectors Z and Step 5.2.3: respectively establishing a reliability model to obtain the uncertainty parameter points of the limit state function on the reliability index surface, and specifically is as follows: s.t. ||Z|| = β t s.t. ||Z|| = β t Where, β t This represents the reliability index corresponding to the constraint performance function, and is derived from its corresponding reliability requirement P. t Calculations show that β t =Φ -1 (P t );Φ -1 (·) denotes the inverse function of the standard normal cumulative distribution function; Step 5.2.4: Calculate the stress limit state function With the band limit state function About each standard normal variable Z t Sensitivity of each worst limit state function with respect to each standard normal variable Zt, t = 1, 2,..., 5; the uncertainty parameter point vector of each worst limit state function on the reliability index surface is calculated by the gradient search algorithm And And the material attribute vector corresponding to the above uncertainty parameter point is obtained by using the Nataf inverse transformation And And the worst limit state function And Step 5.2.5: Based on each worst limit state function and Establishing the normalized form of the limit state approximation function Specifically: Based on Nataf transformation, the independent standard normal random vector Z is used to replace the related uncertainty vector X to establish a reliability analysis model, and specifically is as follows: s.t. ||Z|| = β t Step 5.2.6: solving the function Regarding the sensitivity of the variable Z t Based on the above sensitivity solving results combined with the gradient search algorithm, the calculation of The minimum value of Z corresponds to the uncertainty parameter point approximation value vector Z on the reliability index surface ap , through the Nataf inverse transformation, the corresponding material attribute vector X ap And the worst limit state approximation function And 10. The method of isogemetric robust topology optimization of composite structures considering band constraints of claim 1, wherein: In step 6, the objective function, each worst limit state approximation function and the sensitivity of the utilization rate of the optimization space with respect to the density of each unit control point are calculated; according to the calculation result of the sensitivity of the control point density, the control point density is updated in combination with the moving asymptote method, and a convergence condition is checked, i.e. the average of the relative change degree of the objective function in consecutive iterations is not higher than a threshold value, and if the convergence condition is not met, the iteration returns to step 3.

Images