Isogeometric robust topological optimization method of composite material structure considering frequency band constraint

Abstract

The invention discloses an isogeometric robust topological optimization method of a composite material structure considering frequency band constraints, which comprises the following steps of: firstly, constructing an isogeometric robust topological optimization model of the composite material structure considering material attribute related uncertainty and frequency band constraints; the mean value and the standard deviation of the dynamic flexibility of the structure under the related uncertainty are calculated through Nataf transformation and Gauss-Hermite integration; then solving a worst limit state approximation function and an uncertainty parameter point vector corresponding to each limit state function through a gradient search algorithm; and finally, calculating the sensitivities of the target function, the worst limit state approximation function and the volume constraint function, and obtaining a robust and optimal topological structure meeting the dynamic performance requirement in combination with a moving asymptote method. The isogeometric robust topological optimization model of the composite material structure established by the invention can truly reflect the correlation and dynamic performance requirements of the material attribute uncertainty, can efficiently obtain the optimization result, and has good engineering application value.

Description

Technical Field

[0001] The present invention belongs to the field of structural optimization, and particularly relates to an isogeometric robust topology optimization method for composite material structures considering frequency band constraints. Background Art

[0002] In the fields of equipment manufacturing, automobiles, etc., during the working process of equipment such as milling machines and engines, their key structures often bear harmonic loads of a certain frequency. When their natural frequencies approach the frequencies of the harmonic loads they receive, resonance will occur, which may cause large deformations or fractures in the structure. Therefore, the research on topology optimization of such structures needs to consider the reliability of their dynamic performance.

[0003] Existing research on dynamic topology optimization often focuses on single-material structures and rarely conducts research on composite material structures. Due to the complexity of the composite material preparation process, there are uncertainties in their material properties, and there are correlations between the properties of the same type of raw materials. In addition, the current research on dynamic topology optimization of composite material structures mainly uses structural dynamic compliance, displacement, and fundamental frequency as optimization performance indicators, and rarely uses 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 related uncertainties to meet the actual engineering needs. Summary of the Invention

[0004] In order to solve the deficiencies of the existing technology and achieve the purpose of avoiding the frequencies of the harmonic loads received by the composite material structure and improving the stiffness, strength, and vibration resistance of the structure, the present invention adopts the following technical solutions:

[0005] An isogeometric robust topology optimization method for composite material structures considering frequency band constraints, comprising the following steps:

[0006] Step 1: Construct an isogeometric analysis model of the composite material structure to obtain the density of the control points;

[0007] Step 2: Obtain the relevant material uncertainty descriptions of the composite material structure; construct frequency band constraints according to the vibration resistance requirements of the composite material structure; use the density of the control points as the design variables, and combine the uncertainty descriptions, frequency band constraints, dynamic stress constraints, and optimization space utilization constraints to construct an isogeometric robust topology optimization model of the composite material structure considering frequency band constraints. Under the constraint conditions, obtain the mean value and standard deviation of the minimum dynamic compliance according to the displacement and dynamic stiffness of the control points;

[0008] Step 3: Construct the dynamic stiffness of the composite material structure under harmonic loads;

[0009] Step 4: Calculate the global failure coefficient of the composite material structure considering relevant uncertainties 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 compliance under the influence of relevant uncertainties; meanwhile, based on the frequency band constraint and the global failure coefficient, solve the uncertainty parameter points corresponding to each worst limit state function through a gradient search algorithm to calculate the corresponding material properties and the worst limit state approximation function;

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

[0012] Further, the geometrically robust topology optimization model of the composite material structure considering the frequency band constraint is constructed in Step 2 as follows:

[0013]

[0014] s.t.

[0015] g V (ρ) = V(ρ) / V 0 ≤ [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] Among them, represents the mean of the dynamic compliance, represents the standard deviation of the dynamic compliance, N e represents the number of elements, u e (ρ,X) represents the displacement vector of element e, and the element corresponding to e is the control point, hereinafter abbreviated as Denote the dynamic stiffness matrix of element e, hereinafter abbreviated as g V (ρ) represents the optimization space utilization constraint function, V(ρ) represents the current optimized space size, V 0 and [V] respectively represent the design domain size and the upper limit value of the space utilization constraint, P(·) represents the probability calculation function, and [ω TH respectively represent the global structural failure coefficient composed of the Tsai-Hill criterion and the P-mean aggregation function and the set upper limit value of the global failure coefficient, P t represents the given constraint reliability requirement, g ω (ρ,X) represents the frequency band constraint performance function, exp(·) represents the exponential function, 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, φ k represents the k-th order structural mode, ω k (ρ,X) represents the k-th order natural frequency, M(ρ) represents the global mass matrix, N ω represents the pre-selected number of natural frequency orders, ρ represents the material density corresponding to the control point, ρ i,j represents the control point P i,j corresponding material density, X represents the probability vector that satisfies the normal distribution, V f , E f , ν f , E m , ν m respectively represent the volume fraction of the reinforcement phase, Young's modulus and Poisson's ratio of the composite material structure, and Young's modulus and matrix Poisson's ratio of the matrix material.

[0023] Furthermore, in step 2, considering the uncertainties of the volume fraction V f , Young's modulus E f and Poisson's ratio ν f of the reinforcement phase of the composite material structure, Young's modulus E m and matrix Poisson's ratio ν m of the matrix material, and the correlation between the same material properties, describe them as a probability vector X = (V f , E m , ν m , E f , ν f ) T , and where and respectively represent the mean values of Young's modulus and Poisson's ratio of the matrix material, and respectively 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 respectively represent the means of the Young's modulus and Poisson's ratio of the reinforcement phase, and respectively represent the standard deviations of the Young's modulus and Poisson's ratio of the reinforcement phase, represents the correlation coefficient of the Young's modulus and Poisson's ratio of the reinforcement phase.

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

[0025]

[0026] wherein, g ω (ρ,X) represents the frequency band constraint performance function, ρ represents the material density corresponding to the control point, ω k (ρ,X) represents the k-th (k = 1, 2,..., N ω ) natural frequency, hereinafter abbreviated as N ω represents the pre-selected order number of the natural frequency, p b represents the coefficient of the P-mean aggregation function corresponding to the frequency band constraint, represents the smoothed Heaviside function of the k-th natural frequency with respect to the constraint frequency band, specifically:

[0027]

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

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

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

[0031] Furthermore, the said step 3 includes the following steps:

[0032] Step 3.1: Calculate the element static stiffness matrix k e (ρ,X) and the element mass matrix m e(ρ);

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

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

[0035] Furthermore, the isogeometric analysis model in Step 1 is a non-uniform rational B-spline (NURBS) grid model. The overall grid model is subjected to Gaussian subdivision, and the density of the Gaussian point at the ii-th row and jj-th column in the e-th element is obtained from the control points corresponding to the element e. density

[0036] Furthermore, Step 4 includes the following steps:

[0037] Step 4.1: Calculate the dynamic stress at the Gaussian point of the composite structure.

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

[0039]

[0040] where ρ represents the material density corresponding to the control points, X represents the probability vector that follows a normal distribution, represents the failure coefficient at the Gaussian point of the element considering the relevant uncertainty, N e , N G respectively represent the number of elements in the NURBS grid model and the number of rows or columns of the Gaussian points inside the element. Specifically:

[0041]

[0042] where represents the failure coefficient at the Gaussian point of the element considering the relevant uncertainty, hereinafter abbreviated as represents the strength coefficient matrix of the material, hereinafter abbreviated as p t represents the P-mean aggregation function coefficient of the stress constraint function.

[0043] Further, in step 5, the mean and standard deviation of the structural dynamic flexibility under the influence of relevant uncertainties are calculated based on the Nataf transformation, specifically including the following steps:

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

[0045]

[0046] Step 5.1.2: Calculate the structural dynamic flexibility c when m = 1, 2 d The approximate value of the origin moment of (ρ, X), specifically: is as follows:

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

[0048] Z = (Z 1 , Z 2 , Z 3 , Z 4 , Z 5 ) T ;

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

[0050] as Specifically:

[0051]

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

[0053] Step 5.1.2.3: Calculate the approximate value of the origin moment when m = 1, 2 based on the Gauss-Hermite integral, specifically:

[0054]

[0055] where H represents the number of Gauss-Hermite integration points selected; w h represents the weight corresponding to the hth Gauss-Hermite integration point; where is the corresponding value of the h-th Gauss-Hermite integration point; represents the nominal value of the structural dynamic flexibility obtained when the variables in vector Z are all taken as the mean value;

[0056] Step 5.1.2.4: Based on the Nataf inverse transformation, the vector obtained by taking values at the Gauss-Hermite integration points determine the corresponding material uncertainty vector Use represents the approximate value of the origin moment when m = 1, 2;

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

[0058]

[0059] Among them, represents the mean value of the dynamic flexibility, represents the standard deviation of the dynamic flexibility.

[0060] Furthermore, in the above Step 5, the uncertainty parameter points corresponding to each worst limit state function are solved by the 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, expressed as:

[0062]

[0063] Among them, ρ represents the material density corresponding to the control point, X represents the probability vector that satisfies the normal distribution, exp(·) represents the exponential function, represents the limit state function corresponding to the stress constraint performance function, represents 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 the independent standard normal vector Z and

[0065] Step 5.2.3: Establish reliability models respectively to obtain the uncertainty parameter points of the limit state function on the reliability index surface, specifically as follows:

[0066]

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

[0068]

[0069] such that ||Z|| = β t

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

[0071] Step 5.2.4: Calculate the sensitivity of the stress limit state function and the frequency band limit state function with respect to each standard normal variable Z t (t = 1, 2,..., 5); through the gradient search algorithm, calculate the uncertainty parameter point vectors and of each worst limit state function on the reliability index surface, and use the Nataf inverse transformation to obtain the material property vectors and corresponding to the above uncertainty parameter points, as well as the worst limit state functions and

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

[0073]

[0074] Based on the Nataf transformation, use the independent standard normal random vector Z to replace the relevant uncertainty vector X, and establish a reliability analysis model, specifically:

[0075]

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

[0077] Step 5.2.6: Solve the sensitivity of the function with respect to the variable Z t Based on the above sensitivity solution results and combined with the gradient search algorithm, calculate the approximation value vector Z of the uncertainty parameter point on the reliability index surface corresponding to Z when taking the minimum value ap, calculate the corresponding material property vector X through the Nataf inverse transformation ap and the worst limit state approximation function and

[0078] Furthermore, in the step 6, calculate the objective function for each worst limit state approximation function and and the optimization space utilization constraint g V (ρ) with respect to the density ρ of each element control point e,a,b sensitivity; according to the calculation result of the control point density sensitivity, combine the moving asymptote method to update the control point density, check the convergence condition, that is, the average value of the relative change degree of the objective function in several consecutive iteration steps is not higher than the threshold. If the convergence condition is not satisfied, return to step 3 for iteration.

[0079] The advantages and beneficial effects of the present invention are as follows:

[0080] Compared with the isogeometric robust topology optimization method without frequency band constraints, the present invention considers the influence of harmonic loads, and the natural frequency of the optimized structure can avoid the frequency of the harmonic loads received, having good stiffness, strength and vibration resistance; the present invention describes material uncertainty as a random variable satisfying the bivariate normal distribution, uses the correlation coefficient to describe the correlation between uncertainties, and proposes a method for solving the statistical value of the target performance based on the Nataf transformation. After converting the correlated variables into independent standard normal variables, the statistical moments of the target performance are solved by combining the single-variable dimensionality reduction method and the Guass-Hermite transformation. Compared with the method that does not consider the correlation, it can more accurately evaluate the statistical value of the target performance; the present invention proposes a multi-constraint performance reliability evaluation method considering the approximation of the worst constraint performance. Using the worst limit state function obtained from each probability constraint, a normalized reliability analysis model is constructed to obtain the worst limit state approximation function, so that each probability constraint is in the same optimization background. Compared with the method without approximation processing, it can obtain better structural target performance. Description of the Drawings

[0081] Figure 1 is the flowchart of the method of the embodiment of the present invention.

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

[0083] Figure 2b is a schematic diagram of the initial design of the engine connecting rod structure in the embodiment of the present invention.

[0084] Figure 3 is the optimal topology diagram of the engine connecting rod structure in the embodiment of the present invention. Specific implementation manners

[0085] The following is a detailed description of the specific implementation manners of the present invention with reference to the accompanying drawings. It should be understood that the specific implementation manners described herein are only for explaining and understanding the present invention, and are not used to limit the present invention.

[0086] As Figure 1 shown, for the isogeometric robust topology optimization method of a composite material structure considering frequency band constraints, considering that the composite material structure is affected by harmonic loads and related uncertainties under actual working conditions, the material property uncertainties with sufficient samples are described as variables satisfying the bivariate normal distribution, and the correlation between different material properties of the same material is described by the correlation coefficient. At the same time, with the dynamic compliance mean and standard deviation as the objective functions, an isogeometric robust topology optimization model of the composite material structure under dynamic stress constraints, frequency band constraints and material usage constraints is established, which fully reflects the design requirements for high stiffness, high strength, high vibration resistance and light weight of the composite material structure. At the same time, a reliability evaluation method considering the approximation of the worst constraint performance is proposed to obtain a unified approximation point of the worst constraint performance. On this basis, the sensitivities of the objective function, the limit state approximation function and the volume constraint function are derived and solved, and the optimal topology design scheme is obtained by using the moving asymptote method. The optimization method specifically includes the following steps:

[0087] Step 1: Construct a non-uniform rational B-splines (NURBS) isogeometric analysis model of the composite material structure, perform Gaussian subdivision on the overall NURBS grid model, and obtain the density of the \(i\)-th row and \(j\)-th column of the \(e\)-th element in the \(e\)-th element from the corresponding control points of the element \(e\) density

[0088] As Figure 2a 、 Figure 2b shown, the design information in the embodiment of the present invention is the actual application data in the topology optimization design of a connecting rod in a certain complex component made of TM800s / M21 carbon fiber reinforced composite material. The left end of the connecting rod is fixed, and the right end is subjected to three harmonic loads \(F\) 1 、\(F\) 2 and \(F\) 3 in the horizontal left direction with a frequency of 80 Hz, and the load amplitudes \(f\) 1 \(=f\) 2 \(=f\) 3 \(= 300\ N\); Figure 2b The upper base length of the design domain in [the figure] is 27, the lower base length is 36, and the height is 130 mm.

[0089] On the basis of constructing the NURBS curve, two-direction knot vectors \(u=(u0 , u 1 ,..., u 101+2+1 ), and v = (v 0 , v 1 ,..., v 21+2+1 ), introduce 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; perform Gaussian subdivision on the overall NURBS grid model, and the Gaussian points of the ii-th row and jj-th column of the e-th element can be obtained from the control points corresponding to the element e Density

[0090] Step 2: Construct an isogeometric robust topology optimization model for composite structures considering relevant uncertainties and frequency band constraints, which specifically includes the following steps:

[0091] Step 2.1: Considering the uncertainties of the volume fraction V f , Young's modulus E f and Poisson's ratio ν f of the reinforcement phase of the composite structure, the Young's modulus E m and matrix Poisson's ratio ν m of the matrix material, and the correlation between the same material properties, describe them as a probability vector X = (V f , E m , ν m , E f , ν f ) T , and where and respectively represent the means of the Young's modulus and Poisson's ratio of the matrix material, and respectively represent the standard deviations of the Young's modulus and Poisson's ratio of the matrix material, represents the correlation coefficient between the Young's modulus and Poisson's ratio of the matrix material; and respectively represent the means of the Young's modulus and Poisson's ratio of the reinforcement phase, and respectively represent the standard deviations of the Young's modulus and Poisson's ratio of the reinforcement phase, represents the correlation coefficient between the Young's modulus and Poisson's ratio of the reinforcement phase.

[0092] In the embodiment of the present invention, considering the volume fraction V f , Young's modulus E f and Poisson's ratio ν f, the Young's modulus E of the matrix material m and the Poisson's ratio ν m of the uncertainty and the correlation between the Poisson's ratio and Young's modulus of the same material, which is described as V f ~N(0.6, 0.01 2 ), (E m , ν m )~N(3500, 0.35, 150 2 , 0.01 2 , 0.8) and (E f , ν f )~N(240000, 0.30, 4800 2 , 0.01 2 ), 0.8).

[0093] Step 2.2: According to the anti-vibration performance requirements of the composite material structure, construct a frequency band constraint function based on the smooth Heaviside function and the P-mean aggregation function, specifically:

[0094]

[0095] 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 (k = 1, 2,..., N ω ) natural frequency, hereinafter abbreviated as N ω represents the order of the natural frequency pre-selected; p b represents the coefficient of the P-mean function corresponding to the frequency band constraint; represents the smooth Heaviside function of the k-th natural frequency with respect to the constraint frequency band, specifically:

[0096]

[0097] where ω c represents the midpoint of the working frequency band determined according to the harmonic excitation, δ represents the coefficient of the frequency band constraint performance function, ω w represents half of the width of the restricted frequency band region, 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, construct a geometrically robust topology optimization model of the composite structure considering the frequency band constraint:

[0101]

[0102] Among them, represents the mean value of dynamic compliance, represents the standard deviation of dynamic compliance, N e represents the number of elements, u e (ρ, X) represents the displacement vector of element e, hereinafter abbreviated as represents the dynamic stiffness matrix of element e, hereinafter abbreviated as g V (ρ) represents the constraint function of the optimization space utilization rate, V(ρ) represents the current optimization space size, V 0 and [V] respectively represent the design domain size and the upper limit value of the space utilization rate constraint, P(·) represents the probability calculation function, and [ω TH respectively represent the global structural failure coefficient composed of the Tsai-Hill criterion and the P-mean aggregation function and the set upper limit value of the 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 k-th order structural mode, M(ρ) represents the overall mass matrix, N ω represents the preselected number of natural frequency orders, ρ i,j is the material density corresponding to the control point P i,j corresponding to.

[0103] In the embodiment of the present invention, and represent the mean value and standard deviation of the dynamic compliance 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 overall harmonic load amplitude vector applied to the connecting rod; K(ρ, X) represents the overall static stiffness matrix of the connecting rod; φ k represents the k-th order mode 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 harmonic load, specifically:

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

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

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

[0108] Step 4: Calculate the global failure coefficient of the composite structure (connecting rod) considering relevant uncertainties, specifically:

[0109] Step 4.1: Calculate the dynamic stress of the composite structure (connecting rod) at the element Gaussian point at

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

[0111]

[0112] Among them, represents the failure coefficient at the element Gaussian point considering relevant uncertainties, N e = 2000, N G = 3 respectively represent the number of elements of the NURBS mesh model and the number of rows or columns of the element Gaussian points, specifically:

[0113]

[0114] Among them, represents the failure coefficient at the element Gaussian point considering relevant uncertainties, 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 structural dynamic flexibility under the influence of relevant uncertainties 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] Among them, calculate the mean and standard deviation of the structural dynamic flexibility under the influence of relevant uncertainties 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 structural dynamic flexibility c when m = 1, 2 d The origin moment of (ρ, X) The approximate value of, specifically:

[0120] Step 5.1.2.1: Based on the Nataf transformation, transform the relevant uncertainty vector X into an independent standard normal random vector Z = (Z 1 , Z 2 , Z 3 , Z 4 , Z 5 ) T ;

[0121] Step 5.1.2.2: Based on the univariate dimension reduction method, rewrite the origin moment of the structural dynamic flexibility using the independent standard normal random vector Z as Specifically:

[0122]

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

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

[0125]

[0126] where H represents the number of Gauss-Hermite integration points selected; w h represents the weight corresponding to the h-th Gauss-Hermite integration point; where is the corresponding value of the h-th Gauss-Hermite integration point; represents the nominal value of the structural dynamic flexibility obtained by taking the mean of each variable in the vector Z;

[0127] Step 5.1.2.4: The vector obtained by taking values at Gauss-Hermite integration points based on the Nataf inverse transformation determine the corresponding material uncertainty vector Use to represent the approximate values of the origin moments when m = 1, 2;

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

[0129]

[0130] Solve the vector of uncertainty parameter points corresponding to each worst limit state function through the gradient search algorithm, specifically:

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

[0132]

[0133]

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

[0135] Step 5.2.3: Establish reliability models 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, which is calculated from its corresponding reliability requirement P t , β 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 band limit state function with respect to each standard normal variable Z tThe sensitivity of (t = 1, 2,..., 5); through the gradient search algorithm, calculate the vector of uncertainty parameter points (Minimum Performance Target Points, MPTP) of each worst limit state function on the reliability index surface and And use the Nataf inverse transformation to obtain the material property vectors corresponding to the above uncertainty parameter points and And the worst limit state function and

[0139] Step 5.2.5: Based on each worst limit state function and Establish a normalized limit state approximation function, specifically:

[0140]

[0141] Based on the Nataf transformation, use the independent standard normal random vector Z to replace the correlated uncertainty vector X to establish a reliability analysis model, specifically:

[0142]

[0143] Step 5.2.6: Solve the function Regarding the variable Z t Of the sensitivity, based on the above sensitivity solution results combined with the gradient search algorithm, calculate The approximation value vector Z of the uncertainty parameter point on the reliability index surface corresponding to Z when taking the minimum value ap , through the Nataf inverse transformation, calculate the corresponding material property vector X ap And the worst limit state approximation function And

[0144] Step 6: Calculate the objective function Each worst limit state approximation function And And the volume constraint function g V (ρ) Regarding the density ρ of each element control point e,a,b Of the sensitivity; according to the above calculation results of the sensitivity of the objective function, the worst limit state approximation function and the volume constraint function with respect to the control point density, update the control point density by combining the moving asymptote method; check the convergence condition (the average relative change degree of the objective function Within several consecutive iteration steps is not higher than the threshold), if the convergence condition is not satisfied, return to step 3 for iteration.

[0145] In the embodiment of the present invention, the convergence condition is checked that the average relative change degree of the objective function within three consecutive iteration steps is lower than 0.05; the convergence condition is satisfied after the 30th iteration, and the average value of the dynamic flexibility of the optimized structure and the standard deviation of the 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 present invention. The optimal topology of the engine connecting rod structure is shown in Figure 3 as follows.

[0146] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than limiting them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. An isogeometric robust topology optimization method for composite structures considering frequency band constraints, characterized by The steps include: Step 1: Construct an isogeometric analysis model of the composite structure to obtain the density of control points; Step 2: Obtain the uncertainty description of the relevant materials of the composite structure; construct the frequency band constraint according to the vibration resistance requirement of the composite structure; take the density of the control points as the design variable, combine the uncertainty description, frequency band constraint, dynamic stress constraint and optimization space utilization constraint, and construct an isogeometric robust topology optimization model of the composite structure considering the frequency band constraint. Under the constraint conditions, according to the displacement and dynamic stiffness of the control points, obtain the minimized dynamic flexibility mean and dynamic flexibility standard deviation; Step 3: Construct the dynamic stiffness of the composite structure under simple harmonic loads; Step 4: Calculate the global failure coefficient of the composite structure taking into account the relevant uncertainties by controlling the dynamic stress of the composite structure at the control point; 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 points corresponding to each worst limit state function through a gradient search algorithm to calculate the corresponding material properties and the worst limit state approximation function; 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 variables.

2. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: In step 2, a geometrically robust topology optimization model of composite material structure considering frequency band constraints is constructed as follows: in, represents the mean value of dynamic flexibility, represents the standard deviation of dynamic compliance, N e Indicates the number of units, u e (ρ,X) represents the displacement vector of unit e, and the unit corresponding to the control point represented by e, represents the dynamic stiffness matrix of element e, g V (ρ) represents the optimization space utilization constraint function, V(ρ) represents the current optimization space size, V0 and [V] represent the upper limit of the design domain size and space utilization constraint, respectively, P(·) represents the probability calculation function, and [ω TH ] represent the global failure coefficient of the structure and the upper limit of the global failure coefficient, respectively. t represents the given constraint reliability requirement, g ω (ρ,X) represents the band-constrained performance function, exp(·) represents the exponential function, and 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 k-th structural mode, ω k (ρ,X) represents the kth order natural frequency, M(ρ) represents the overall mass matrix, N ω represents the pre-selected natural frequency order, ρ represents the material density corresponding to the control point, ρ i,j represents the control point P i,j The corresponding material density, X represents the probability vector that satisfies the normal distribution, V f 、E f , ν f 、E m , ν m They represent the volume fraction of the reinforcement phase, Young's modulus and Poisson's ratio of the composite material structure, as well as the Young's modulus and Poisson's ratio of the matrix material, respectively.

3. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 2, characterized in that: In step 2, the volume fraction V of the composite material structure reinforcement phase is considered. f , Young's modulus E f and Poisson's ratio ν f , Young's modulus of matrix material E m Poisson's ratio of the matrix m The uncertainty and correlation between the properties of the same material are described as a probability vector X=(V f ,E m ,ν m ,E f ,ν f ) T , and in and represent the mean values ​​of Young’s modulus and Poisson’s ratio of the matrix material, and represent the standard deviation of Young’s modulus and Poisson’s ratio of the matrix material, Represents the correlation coefficient between Young's modulus and Poisson's ratio of the matrix material; and are the means of Young’s modulus and Poisson’s ratio of the reinforcement phase, and represent the standard deviation of Young’s modulus and Poisson’s ratio of the reinforcement phase, It represents the correlation coefficient between Young's modulus and Poisson's ratio of the reinforcement phase.

4. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 2, characterized in that: In step 2, according to the vibration resistance requirements of the composite material structure, a frequency band constraint function based on a smoothing function and a cohesive function is constructed, specifically: Among them, g ω (ρ,X) represents the band constraint performance function, ρ represents the material density corresponding to the control point, ω k (ρ,X) represents the kth (k=1,2,…,N ω ) order natural frequency, hereinafter abbreviated as N ω represents the pre-selected natural frequency order, p b represents the coefficient of the aggregation function corresponding to the frequency band constraint, represents the smooth function of the k-th order natural frequency with respect to the constrained frequency band, specifically: Among them, ω c represents the midpoint of the operating frequency band determined according to simple harmonic excitation, δ represents the coefficient of the frequency band constraint performance function, ω w Indicates half of the width of the restricted frequency band area, specifically: oh w =(ω upp -oh low ) / 2 Among them, ω low Represents the lower limit of the frequency band constraint; ω upp Indicates the upper limit value of the band constraint.

5. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: The step 3 comprises the following steps: Step 3.1: Calculate the unit static stiffness matrix and unit mass matrix of the composite structure based on the plate model; Step 3.2: Based on the Rayleigh damping method, the structural overall stiffness matrix, the structural overall mass matrix and the overall damping matrix are calculated through the unit stiffness matrix and the unit mass matrix; Step 3.3: Calculate the overall dynamic stiffness matrix of the composite structure based on the overall stiffness matrix, the overall mass matrix of the structure and the overall damping matrix.

6. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: The isogeometric analysis model in step 1 is a non-uniform rational B-spline grid model. The grid model is Gaussian subdivided as a whole, and the Gaussian points of the units in the iith row and jjth column in the eth unit are obtained from the control points corresponding to the unit e. Density 7. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 6, characterized in that: The step 4 comprises the following steps: Step 4.1: Calculate the unit Gaussian points Dynamic stress in composite structures Step 4.2: Calculate the global failure coefficient of the composite structure based on the failure criterion and the aggregation function Specifically: Among them, ρ represents the material density corresponding to the control point, X represents the probability vector that satisfies the normal distribution, represents the failure coefficient at the Gaussian point of the unit considering the relevant uncertainty, N e 、N G They represent the number of cells in the NURBS mesh model and the number of rows or columns of Gaussian points inside the cells, respectively. Specifically: in, represents the failure coefficient at the Gaussian point of the unit considering the relevant uncertainty; Represents the strength coefficient matrix of the material; p t Represents the aggregation function coefficient of the stress constraint function.

8. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: In step 5, the mean and standard deviation of the dynamic flexibility of the structure under the influence of relevant uncertainties are calculated based on the Nataf transformation, which specifically includes the following steps: Step 5.1.1: Introduce weight coefficient w X Define the weighted objective function Step 5.1.2: Calculate the dynamic flexibility c of the structure when m = 1,2 d (ρ,X) origin moment The approximate value of is: Step 5.1.2.1: Based on the Nataf transformation, transform the related uncertainty vector X into an independent standard normal random vector Z = (Z1, Z2, Z3, Z4, Z5) T ; Step 5.1.2.2: Based on the univariate dimensionality reduction method, use the independent standard normal random vector Z to rewrite the origin moment of the structural dynamic flexibility for Specifically: Where m represents the order of the target performance origin moment. represents independent standard normal random variables Z t The marginal probability density distribution function of N (ρ, Z) represents the dynamic compliance of the composite structure represented by an independent standard normal random vector Z; Step 5.1.2.3: Calculate the approximate origin moment based on the integral when m = 1, 2, specifically: Where H represents the number of selected integral points; w h represents the weight corresponding to the hth integral point; in is the corresponding value of the hth integration point; It represents the nominal value of the structural dynamic flexibility obtained by taking the mean value of each variable in the vector Z; Step 5.1.2.4: Based on the Nataf inverse transform, the vector obtained by taking the value of the integral point Determine the corresponding material uncertainty vector use It represents the approximate value of the origin moment when m=1,2; Step 5.1.3: Calculate the mean and standard deviation of the structural dynamic flexibility under the influence of the relevant uncertainty X, specifically: in, represents the mean value of dynamic flexibility, Represents the standard deviation of dynamic compliance.

9. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: In step 5, the uncertainty parameter points corresponding to each worst limit state function are solved by a gradient search algorithm, which specifically includes the following steps: Step 5.2.1: Based on the band-constrained performance function g ω (ρ,X) and stress constraint performance function Construct the limit state function, expressed as: Where ρ represents the material density corresponding to the control point, X represents the probability vector that satisfies the normal distribution, and exp(·) represents the exponential function. represents the limit state function corresponding to the stress constraint performance function, represents the limit state function corresponding to the band constraint performance function; Step 5.2.2: Rewrite the limit state function in the form of an independent standard normal vector Z and Step 5.2.3: Establish reliability models respectively to obtain the uncertainty parameter points of the limit state function on the reliability index surface, specifically: Among them, β t Represents the reliability index corresponding to the constraint performance function, and its corresponding reliability requirement P t Calculated, β t =Φ -1 (P t );Φ -1 (·) represents the inverse function of the standard normal cumulative distribution function; Step 5.2.4: Calculate the stress limit state function Band limit state function About each standard normal variable Z t (t=1,2,...,5) sensitivity; by using the gradient search algorithm, the uncertainty parameter point vector of each worst limit state function on the reliability index surface is calculated and The material property vector corresponding to the above uncertainty parameter point is obtained by using the Nataf inverse transform and And the worst limit state function and Step 5.2.5: Based on the worst limit state function and Establish 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, which is as follows: Step 5.2.6: Solve the function About variable Z t The sensitivity of is calculated based on the above sensitivity solution combined with the gradient search algorithm. The approximate value vector Z of the uncertainty parameter point on the reliability index surface corresponding to Z when taking the minimum value ap , through the Nataf inverse transform, calculate the corresponding material property vector X ap And the worst limit state approximation function 10. The isogeometric robust topology optimization method for composite material structures considering frequency band constraints according to claim 1, characterized in that: In step 6, the sensitivity of the objective function, the worst limit state approximation functions and the optimization space utilization constraint to the density of the control points of each unit is calculated; based on the calculation result of the control point density sensitivity, the control point density is updated in combination with the moving asymptote method, and the convergence condition is checked, that is, the average relative change degree of the objective function in several consecutive iteration steps is not higher than the threshold. If the convergence condition is not met, return to step 3 to iterate.

Images