Isogeometric topological optimization method based on dynamic response driving and application thereof

Abstract

The invention belongs to the field of material structure optimization, and particularly relates to an isogeometric topological optimization method based on dynamic response driving, which comprises the following steps of: (1) establishing a thin shell mid-plane geometric representation model based on a T spline and Bessel extraction technology; (2) establishing a frequency response-oriented isogeometric shell unit based on the Kirchhoff theory; (3) solving the aerodynamic response of the surface of the structure based on a high-order panel method of a T spline; (4) constructing a topology description model of the thin shell based on a Bessel unit; (5) establishing a structure dynamic flexibility minimum topological optimization model oriented to the dynamic response of the pneumatic thin shell; (6) dynamic flexibility minimization sensitivity analysis and derivation; and (7) updating optimization model design variables based on an optimal criterion method. According to the method provided by the invention, a pneumatic coupling thin shell dynamic isogeometric analysis framework is constructed, and the problems of poor geometric adaptability, difficulty in complex thin shell dynamic design, lack of pneumatic environment related to optimization and the like in an existing method are effectively solved.

Description

Technical Field

[0001] This invention belongs to the field of material structure optimization, and more specifically, relates to a dynamic response-driven isogeometric topology optimization method and its application. Background Technology

[0002] Thin-shell structures, due to their excellent lightweight and high-strength properties, are widely used in aerospace, automotive, shipbuilding, and civil engineering. These structures, with thicknesses much smaller than other characteristic dimensions, are typically used to resist bending, torsion, and multi-directional combined loads, exhibiting complex coupled behavior under static and dynamic environments. Especially in scenarios such as aircraft skins, submarine hulls, and large machinery enclosures, thin-shell structures are frequently subjected to periodic, harmonic, or random excitations during service. Therefore, their dynamic performance directly determines the structure's safety, fatigue life, and vibration and noise characteristics. Thus, designing for the dynamic performance of thin-shell structures has become a key aspect of improving the reliability of high-performance structures. Topology optimization, as an advanced structural design method, can seek the optimal material distribution configuration within the design domain to improve the structure's performance under specific load and constraint conditions. Traditional topology optimization methods mostly rely on finite element modeling. However, when dealing with thin-shell structures with complex curvatures and arbitrary shapes, the finite element method has limitations in geometric modeling accuracy and element division, often leading to discontinuous optimization results or geometric distortions, affecting the feasibility of engineering applications. The development of isogeometric analysis (IGA) has brought new opportunities to structural design. This method unifies computational geometry and analytical models, allowing high-order smooth surfaces to be directly used for structural analysis and optimization, particularly suitable for the accurate description and modeling of complex thin-shell structures. However, traditional isogeometric modeling relies on NURBS, whose global tensor product properties limit the flexible construction of complex design domains. T-splines, as a parametric modeling tool with local tensor product properties, can more flexibly represent complex thin-shell structures with arbitrary geometries, overcoming the shortcomings of NURBS in local refinement and free modeling.

[0003] Meanwhile, aerodynamic loads, as the main external forces that thin-shell structures must face in the aerospace environment, place higher demands on structural optimization design. Most current structural topology optimization methods fail to effectively consider the complex aerodynamic loads they bear in real service environments, resulting in unstable performance of optimization results in real engineering environments. The main reasons why most current thin-shell structure topology optimizations fail to effectively incorporate aerodynamic loads include the following: (1) Difficulty in coupling analysis methods: Aerodynamic loads are usually calculated using fluid dynamics methods. Especially when considering three-dimensional aerodynamic fields, problems such as velocity potential and pressure distribution need to be solved. This has a complex coupling relationship with structural mechanics solutions, making the construction and solution of optimization models extremely complex; (2) High computational cost: Aerodynamic-structural coupling analysis usually requires a large number of iterative solutions, and the computational cost is much higher than that of pure structural optimization problems. This poses a challenge to the existing topology optimization framework, especially when facing high-resolution models; (3) Limited modeling expression: Traditional NURBS or finite element modeling methods are difficult to flexibly represent arbitrarily complex geometries, limiting the ability to freely design thin-shell structures under complex aerodynamic conditions. In addition, existing models suffer from accuracy loss when representing multi-scale and local details, and cannot accurately convey aerodynamic information. Based on the above problems, there is an urgent need for a unified method that integrates aerodynamic analysis, complex geometric expression and efficient structural dynamic performance optimization, so as to truly realize the dynamic optimization design of thin-shell structures oriented towards aerodynamic environment. This also shows that the development of topology optimization methods for thin-shell structures that consider aerodynamic loads and excitation frequencies has important theoretical significance and engineering application value. Summary of the Invention

[0004] To address the shortcomings and improvement needs of existing technologies, this invention provides a dynamic response-driven isogeometric topology optimization method and its application. Its purpose is to overcome technical bottlenecks such as insufficient aerodynamic load coupling capability, low geometric modeling continuity, and difficulty in solving dynamic responses. It breaks through the limitations of traditional methods in modeling and optimizing complex thin-shell structures, achieving accurate simulation and optimization of structural aerodynamic behavior. This further uncovers the dynamic performance advantages of thin-shell structure configurations, effectively improving structural dynamic stiffness and stability, and demonstrating broad application prospects in aerospace, complex equipment, and functional structures.

[0005] To achieve the above objectives, according to one aspect of the present invention, a dynamic response-driven isogeometric topology optimization method is provided, comprising: S1. Model the geometric model of a complex thin-shell structure in CAD software based on T-splines, and export the model; S2. Construct a unified geometric representation model for T-spline geometric models based on Bézier extraction technology; S3. Construct isogeometric thin-shell elements based on Kirchhoff theory to address structural frequency response, and establish a dynamic analysis model for the thin-shell structure. S4. The aerodynamic response of the structural surface is solved using a high-order panel method based on T-splines. The collocation method is used to satisfy the Neumann condition required for aerodynamics, and the source density distribution of the structural surface is solved based on the T-spline mixing function. Furthermore, the velocity potential of the structural surface is described by the Bessel element basis function, thereby obtaining the pressure distribution of the thin shell surface. S5. Apply the aerodynamic pressure distribution as a load directly to the Bessel control points, construct a local density distribution function based on the Bessel elements, and assemble it into a global density distribution function. S6. Solve the element stiffness and mass matrices of the uniformly geometric shell elements using the Gaussian integral method and assemble them into global stiffness and mass matrices. Solve the displacement field of the thin shell under aerodynamic load based on the reduced-order model. S7. Based on the aerodynamic response field, establish an iso-geometric topology optimization model, calculate the objective function and derive the complete analytical sensitivity, iteratively update the design variables based on the optimal criterion method, and finally obtain the optimized configuration after iteration.

[0006] Furthermore, the implementation of step S1 includes: In Rhinoceros 3D software, the geometric model of any thin shell structure can be converted into a T-spline-based surface geometry model by combining it with the Autodesk T-splines plugin; the integration of this plugin with Rhinoceros 3D makes the geometric modeling process visual.

[0007] Furthermore, step S2 is implemented in the following ways: The mid-surface geometry of an arbitrary thin-shell structure can be represented using a T-spline mixture function as follows: In the formula, A represents the vertex number in the T mesh. and It is the local node vector that defines the mixing function, which can be obtained through the Cox-de Boor recursive formula when p = 0. When p≥1 In the formula, It is the i-th node; Based on the Bessel extraction technique, a unified geometric representation model can be constructed for any thin shell structure. Therefore, the T-spline mixture function can be further expressed as: In the formula, For the parameter coordinates in the Bessel parent unit domain, select the parameter coordinates in the interval [-1,1] for Gaussian integration, so as to facilitate subsequent analysis and optimization; Let be a set containing T-spline mixture functions, whose support points are located in Bessel element e, where a represents the relevant local index of the control point, and n is the number of control points; Extract the Bessel operator for element e, assuming the polynomial order is the same in all directions; the extraction matrix can be obtained through the following steps: Establish a one-dimensional extraction matrix of order p and q in each direction. ; Combining in the form of tensor products: ; For special T-spline cases, if the local junction vectors are irregular, then after expanding the T-junction, the effective support region is taken, and the corresponding structure is constructed. ; To obtain the Bessel expression for the T-spline blending function, a cell mapping from the parent cell domain to the physical domain can be defined, described as follows: In the formula, Describes a rational T-spline mixture function of a Bessel expression; This represents a control point matrix containing unit control points; Indicates the unit weight; The T-spline mixture function can be further expressed as: (The diagonal matrix form represents the unit weights.) The control points of the global T-spline can be mapped to the control points of the Bézier element by the transpose of the extraction operator; therefore, the Bézier control points and weights can be represented by the element extraction operator, and the T-spline curve corresponding to the control points and weights is: In the formula, Unit weight The T-spline element is in diagonal matrix form; this means that each T-spline element has a corresponding equivalent Bessel element, and the further expression for the T-spline mixing function is: In the formula, The formulas for calculating the first and second derivatives of the T-spline mixture function with respect to the local parameter coordinates are as follows: as well as The derivatives related to physical coordinates can be calculated using the chain rule.

[0008] Furthermore, step S3 is implemented in the following ways: Based on Kirchhoff's thin shell theory, the geometric shell element is constructed, and its deformation is represented by the mid-face of the shell. Therefore, along the thickness direction on the mid-face... The displacement vector of any point can be expressed as: In the formula, x and These are the position vectors of material points in deformed and undeformed states, respectively. and The unit normal vectors at points in the mid-plane of the thin shell before and after deformation are represented by the following basis vectors: The upper and lower scales of the second-order tensor represent the inverse covariance components, respectively; assuming thickness... Then, higher-order terms can be ignored. To simplify, it is: Similarly, this also applies to deformable configurations. and These are the first and second basic forms, respectively: The expressions for the components of the Green-Lagrange strain tensor can then be further expressed as: In the formula, Represents mid-surface strain, which describes mid-surface stretching or compression. Representing the change in curvature, used to describe bending, is: The virtual work principle for thin-shell structures is derived using the Galerkin weak form, and the weak form is obtained through variational methods as follows: In the formula, f represents body force, and ρ represents material density. and These are the fourth-order tensors for membrane stiffness and bending stiffness, respectively, which can be derived from the fourth-order tensor of the material. get: In isogeometric analysis, the element displacement field and virtual displacement field Interpolation using the same set of basis functions can be described as follows: In the formula, The total number of basis functions, , To control the displacement of the control point; according to Kirchhoff's theory, the membrane strain and curvature change can be further expressed as: In the formula, Indicates the displacement of the control point; The membrane strain matrix is... The bending strain matrix is: Considering the Rayleigh damping integral over the entire parameter domain and assembled, the discretized dynamic control equations are obtained as follows: In the formula, M is the global mass matrix, C is the Rayleigh damping term, K is the global stiffness matrix, and F is the global external load vector, which is: In the formula, and Here, A and D are the damping coefficients, and A and D are the membrane stiffness and bending stiffness matrices, respectively: Q can be calculated as The global mass and stiffness matrices can be assembled from the element mass and stiffness matrices, which are calculated using Gaussian point numerical integration, as follows: In the formula, The total number of units, Let J be the total number of Gaussian points within the element, and J be the mapping from the parameter domain to the physical domain. Jacobian matrix, These are the corresponding Gaussian point weights; in the proposed method, Based on SIMP, dynamic equivalent material interpolation of thin-shell structures is implemented. Simultaneously, equivalent elastic parameters and linear density interpolation are set to avoid matrix singularities and excessive mass penalty. The results are as follows: In the formula, To minimize the elastic Young's modulus and prevent singularities, Young's modulus of solid materials The penalty index is generally taken as... To suppress intermediate density; For the density of solid materials, The minimum density value is generally taken as... To avoid zero mass in the empty region; based on the above equation, the element mass and stiffness matrix can be further expressed as: The high-order continuity based on T-spline geometry achieves smoothness and accuracy of the structural dynamic response field, improving numerical stability and sensitivity reliability. Therefore, the proposed framework can achieve stable and accurate coupling between aerodynamic loads and structural response in an efficient computational manner.

[0009] Furthermore, in step S4, the characteristic is that, S4.1 Consider a non-viscous, incompressible, irrotational potential flow field with velocity potential. Satisfying the Laplace equation, it is: In the formula, Let the external fluid domain be represented. Applying Green's second identity, the velocity potential is transformed into a boundary integral form, which is: In the formula, For field points (including shell surface points). For the source point location, For free space Green's function, It is a Euclidean norm. Geometric coefficients; Introducing surface element source intensity function Its velocity potential at any point P on the boundary is: In the formula, It is the boundary surface The source intensity depends on the location of the source point; combining the Neumann boundary conditions, we obtain the non-singular expression for the normal fluid velocity as follows: In the formula, The outer surface of the shell, For the far field, The source point normal is given; the right-hand side of the above equation has no singular kernel, which facilitates higher-order numerical integration; this equation can be solved using the panel method to obtain the source intensity. And further calculate its velocity potential energy.

[0010] The proposed method describes the geometry of thin-shell structures based on T-splines. Therefore, the generalized expression for any point P on the spline surface is defined as: Furthermore, the unit normal vector of point P on the T-spline surface can be obtained. : In the formula, and Let denote the derivative of the surface with respect to the directions ξ and η; for ease of subsequent derivation, the expression for the symbol is defined as... The calculation formula is as follows: For the source points on the body surface and If the source point is on the boundary surface, then the normal velocity can be converted into an expression based on the T-spline description, as follows: The body surface, represented by T-splines, is discretized into several panels, with source points defined on each panel. The source intensities can then be derived using boundary integral equations. This method calculates the integrals on the panels based on the Gaussian orthogonal method, as follows: In the formula, , , For the i-th collocation point, and These are the collocation point and the source point normal, respectively. Let g be the Jacobian matrix of the g-th point in the e-th unit. The weights are Gaussian; the assembly is performed over all collocation points i, thus obtaining the weights for the source strength. A system of linear equations with coefficients, where the unknowns are all Gaussian points. and matching points The source strength is obtained by solving for the source strength. The velocity potential can then be further solved. ; Introducing an equipotential auxiliary function, the homogeneous boundary integral equation it satisfies is: In the formula, It is obtained through independent, singular iterative equations. This represents the distance from the source point to the selected origin, which only appears as a scale reference to eliminate odd kernels; the panel method based on T-splines yields the velocity potential reconstruction formula as follows: In the formula, .

[0011] The velocity expression for any point P in the flow field is: This equation solves for the velocity distribution of the fluid on the surface of the thin-shell structure, including the velocity components in all directions; after obtaining the velocity potential energy, its velocity can be further calculated. S4.2 However, when point P is located on the surface of the body, the integral exhibits singularity, therefore the above method is not suitable for solving the velocity field in this case; the higher-order continuity of the Bernstein polynomial allows for direct calculation of the derivative, thereby improving the accuracy of the velocity potential gradient calculation; the velocity potential at a point on the body surface is constructed as follows: In the formula, The Bernstein bivariate basis functions can be calculated using the tensor product of the Bernstein univariate basis functions, and the calculation formula is as follows: In the formula, p defines the order of a polynomial; This represents the velocity potential energy value at the control point; the velocity U of the flow field can be directly obtained as: Therefore, the pressure on the surface of the thin-shell structure can be calculated as follows: In the formula, For free-flow static pressure, For fluid density, For the free flow velocity, according to Bernoulli's equation, the pressure coefficient... It can be calculated as: Projecting the pressure along the normal to the shell surface into surface force density And map this onto the control points to obtain the expression for the aerodynamic load: This solution process ensures the continuity of the surface velocity potential and effectively solves the problems of numerical instability and integral singularity in the panel method calculation process. For a thin shell with uniform geometry, the dynamic equation can be expressed as: In the formula, This is the aerodynamic load vector. For other mechanical loads, for low-speed, lift-dominated problems with relatively small geometric deflections, a unidirectional weakly coupled approximation can be used, where the aerodynamic forces are evaluated with the initial configuration and treated as constants within the iterations: In the formula, the global stiffness matrix K is updated with the design variables, while It does not update with displacement feedback.

[0012] Furthermore, in step S5, the characteristic is that, S5.1. Applying geometric analysis methods such as linear elastic materials to solve the structural response of thin-shell structures, the equilibrium equation of the structural response can be expressed as: In the formula, Let M be the excitation angular frequency, M, C, and K be the global mass, damping, and stiffness matrices, respectively, d be the reset displacement, and F be the global external load vector; let Let the dynamic stiffness matrix be denoted as , then the above equation can be expressed as: In frequency domain dynamics optimization, when the design variables are iteratively updated, the solution is directly repeated on the fully free system. The computational cost is extremely high; the MQSRV model order reduction strategy is adopted to reduce the computational cost while maintaining accuracy. Calculate the cost; will give Target frequency band Dividing the interval into N sub-intervals, the center frequency of each interval is defined as follows: In the formula, each The representative frequency of this interval is used to generate a local quasi-static basis; this uniform division method can flexibly handle multi-peak frequency domain response characteristics, ensuring high-precision approximation in each neighborhood; at the nth center frequency Define the frequency shift stiffness operator The basis derivation is as follows: The recursive process of the basis involves linear equations To avoid high-frequency local mode distortion; at the same time, to maintain the orthogonality of M, the basis is modified and normalized, thus the MQSRV basis... The calculation is as follows: This yields a set of M orthonormal bases, namely: By concatenating the basis vectors obtained from each frequency band, we can obtain the global reduced-order basis of MQSRV, which is: In the formula, Let the total reduced dimension be ; let the approximate solution Combined with the above formula and multiplied on the left The reduced dynamic control equations of MSQRV can be obtained, which significantly reduce the dimensionality of the original system: And it is equivalent to the full-order system in the sense of energy inner product; where, and Let the reduced-order mass, damping, stiffness matrix, and reduced-order external load vector be respectively, which can be derived as follows: Reduced-order dynamic stiffness matrix It can be deduced as: The reduced-order response can then be solved as follows: according to The structural response can then be solved. S5.2 Construct a global density function based on Bézier elements to meet the requirements of topological description of thin shell structures. Specifically, given a vector containing the initial Bézier control density... The smoothing mechanism improves the smoothness of the Bezier control density by using the Shepard function on the Bezier cells, as follows: In the formula, The Shepard function represents the Bessel control point; the smooth, continuous local density distribution function. Density can be controlled by smoothing. and Bessel unit basis functions Represented as: Topological boundary The isomorphic profile of a local DDF can be used. This can be represented as: The common control points between adjacent Bessel elements ensure that local DDFs can be connected into a topological description model to describe the overall structure; therefore, local DDFs can be assembled into a global DDF through the natural connections between adjacent elements to describe the entire thin-shell structure design domain, and its expression is: In the formula, It is the number of Bézier elements in the design domain, symbol […]. Represents the local DDF of the i-th Bessel element; the structural topology of the entire design domain. It can be composed of the local structural topology of each Bessel unit, as follows: Furthermore, step S6 is implemented in the following ways: Substitute the design variables into the S3 analysis model and calculate the reduced-order dynamic stiffness matrix using the MSQRV reduced-order model. Then solve for the reduced-order response. And then according to Calculate the displacement field.

[0013] Furthermore, step S7 is characterized in that, S7.1 The dynamic flexibility minimization topology optimization model is expressed as: In the formula, This refers to the initial node density of the Bézier control points; the design variable values ​​must be between their minimum values. The minimum value between 1 and 0 is to avoid singularities; This represents the total number of design variables. The objective function for calculating global dynamic compliance is given by: Calculate the size of a complex number; The range of excitation frequencies; Represent the global density distribution function; consider the excitation frequency as... Under the action of a simple harmonic load, and These are the complex forms of load and displacement, respectively; i is the imaginary unit; G is the volume constraint. Represents the volume fraction of a solid material. This represents the maximum value of material consumption. Represents a virtual displacement field, belonging to the kinematically acceptable space. ; Indicates Dirichlet boundary The specified displacement vector at the location; k, c, and m are the strain energy, damping energy, and mass bilinear energy, respectively: For a semi-linear load type l, the calculation formula is: In the formula, Represents body forces in thin-shell structures. This is the boundary traction force, belonging to the Neumann boundary. ; S7.2. The sensitivity analysis of the optimization model requires calculating the first derivative of the objective function relative to the density distribution function. Based on the MSQRV order reduction strategy, the objective function can be approximated by the frequency band numerical integral form: objective function If a variable is composed of a real part and an imaginary part, then its absolute value is defined as: In the formula, and Let represent the real and imaginary parts of the functional, respectively. The first derivative of the objective function with respect to the design variables and the dynamic compliance. related: The sensitivity analysis of the objective function is transformed into the derivation of dynamic compliance. The first derivative relative to the design variables; dynamic compliance is used to clearly derive the calculation results. It can be represented as: Dynamic flexibility The first derivative of the relative global density distribution function can be expressed as: In the formula, It is the derivative of the displacement field with respect to the global DDF. It is the derivative of the virtual displacement field with respect to the DDF; Combining the optimization model with the dynamic response analysis of the thin-shell structure, we can conclude that: The first derivative of the bilinear expression for strain energy k and mass m with respect to the global density distribution function is: because Then the following equations hold: Considering that the problem of minimizing dynamic compliance is self-adjoint, the virtual displacement field in the above equation can be eliminated: Then the flexibility The first derivative of the relative global density distribution function is: Further calculation of the first derivatives of the membrane stiffness A and bending stiffness D matrices with respect to the global DDF yields the dynamic compliance. The expression for the derivative of the global density distribution function is: Based on the construction principle of the global density distribution function, the first derivative of the density distribution function with respect to the design variables can be obtained as follows: Therefore, dynamic flexibility The final expression for sensitivity is: The sensitivity of the volume constraint function is: S7.3 Update the control density using the optimal criterion method. That is, design variables, repeat the process until the iteration termination condition is met, and obtain the enhanced structure topology optimization configuration based on the objective function and sensitivity calculated in the last iteration S7.2.

[0014] The application of a dynamic response-driven isogeometric topology optimization method is applied to the design of swept wing skin.

[0015] A computer-readable storage medium includes a stored computer program that performs the dynamic response-driven isogeometric topology optimization method as described above.

[0016] In summary, the above-described technical solutions conceived in this invention can achieve the following beneficial effects: (1) The dynamic response-driven isogeometric topology optimization method and its application provided by this invention, by introducing aerodynamic analysis based on the high-order panel method, can accurately and quickly calculate the velocity potential and aerodynamic pressure field of the thin shell structure surface, and use the aerodynamic field as the main driving factor for topology optimization, realizing the deep coupling between aerodynamic load and structural topology, thus solving the problem of neglecting aerodynamic influence in traditional topology optimization. The high-order panel method used only needs to process the structural surface mesh, avoiding the mesh generation and large-scale solution of the flow field domain, greatly reducing the computational cost of aerodynamic analysis, and is suitable for multi-round iterative processes in topology optimization, effectively improving the overall optimization efficiency.

[0017] (2) The dynamic response-driven isogeometric topology optimization method and its application provided by this invention, based on the T-spline modeling framework and combined with the Bezier element extraction technology, realizes the accurate representation of thin shell structures with arbitrary complex shapes. It has good local refinement ability and high-order continuity, significantly improves the geometric expression accuracy and flexibility in the optimization process, and breaks through the limitations of traditional NURBS and finite element methods in complex geometric expression. This invention does not rely on regular meshes or predefined structural shapes, but is based on the continuous representation of controllable parameters and the construction of driving functions, allowing for exploration of material layout with higher degrees of freedom, and can discover high-performance configurations that are difficult to obtain by traditional design methods, fully releasing the structural potential.

[0018] (3) The dynamic response-driven isogeometric topology optimization method and its application provided by the present invention introduce the MSQRV reduced-order model to realize a high-fidelity description of the dynamic response of multiple frequency bands without having to resolve the complete frequency domain equation in each optimization iteration. Compared with the traditional model update strategy, MSQRV can significantly reduce the system degrees of freedom while maintaining the accuracy of dynamic features, and improve the solution stability and convergence efficiency.

[0019] (4) The dynamic response-driven isogeometric topology optimization method and its application provided by this invention are applicable to the optimization design requirements of typical aerospace equipment structures such as aerospace vehicle skin, aircraft wing structure, satellite bulkhead, and lightweight ship hull. It can also be extended to the optimization of any thin-shell component that needs to take into account both aerodynamic performance and structural dynamic strength. It has broad engineering application value and promotion prospects. Attached Figure Description

[0020] Figure 1 A flowchart illustrating the dynamic response-driven isogeometric topology optimization method and its application provided in this embodiment of the invention; Figure 2 The following is a schematic diagram of a swept wing provided for an embodiment of the present invention, wherein (a) is a schematic diagram of the complete initial design domain of the swept wing, and (b) is a schematic diagram of the aerodynamic pressure distribution of three airfoil sections therein; Figure 3 A schematic diagram of the optimized configuration and iteration history of the swept wing optimization iteration process provided in an embodiment of the present invention; Figure 4 The following is a schematic diagram of the optimized swept wing provided in the embodiment of the present invention, wherein (a) is a schematic diagram of the structural topology, (b) is a schematic diagram of the optimized configuration reconstruction model, (c1) is a schematic diagram of the optimized configuration of the upper surface of the swept wing, and (c2) is a schematic diagram of the optimized configuration of the lower surface of the swept wing. Figure 5 The above is a schematic diagram comparing the frequency response functions of the swept wing before and after optimization, provided in an embodiment of the present invention. (a) is a schematic diagram comparing the frequency response curves without damping, and (b) is a schematic diagram comparing the frequency response curves with Rayleigh damping. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0022] The dynamic response-driven isogeometric topology optimization method provided by this invention has the following process: Figure 1 As shown, it includes the following steps: (1) Modeling the geometric model of complex thin shell structure in CAD software (Rhinoceros 3D and T-splines plugin) based on T-splines.

[0023] (2) Based on the geometric model in step (1), the Bessel extraction technique constructs a unified geometric representation model for the T-spline geometric model. Then, the T-spline mixture function can be expressed as: In the formula, For the parameter coordinates in the Bessel parent unit domain, select the parameter coordinates in the interval [-1,1] for Gaussian integration, so as to facilitate subsequent analysis and optimization. Let be a set containing T-spline mixture functions, whose support points are located in Bessel element e, where a represents the relevant local index of the control point, and n is the number of control points. Extract the Bessel operator for element e. Let It is a vector of the T-spline mixture function of element e, and the T-spline mixture function of element e can be expressed as: In the formula, The formulas for calculating the first and second derivatives of the T-spline mixture function with respect to the local parametric coordinates are as follows: as well as The derivatives related to physical coordinates can be calculated using the chain rule.

[0024] (3) Based on the geometric representation model in step (2), and combined with Kirchhoff's thin shell theory to describe and construct geometric shell elements, the thickness direction along the mid-surface is... The displacement vector of any point can be expressed as In the formula, x and These are the position vectors of material points in deformed (actual) and undeformed (reference) states, respectively. and This represents the unit normal vector at a point in the midplane of the thin shell before and after deformation. The expressions for the components of the Green-Lagrange strain tensor can be further expressed as: In the formula, Represents mid-surface strain, which describes mid-surface stretching or compression. Representing the change in curvature, used to describe bending, is: The virtual work principle for thin-shell structures is derived using the Galerkin weak form, and the weak form is obtained through variational methods as follows: In the formula, f represents body force, and ρ represents material density. and These are the fourth-order tensors for membrane stiffness and bending stiffness, respectively, which can be derived from the fourth-order tensor of the material. express.

[0025] In isogeometric analysis, the element Hexu The displacement field is described by interpolation using the same set of basis functions: In the formula, The total number of basis functions, , To control the displacement of the control point. According to Kirchhoff's theory, the membrane strain and curvature change can be further expressed as: In the formula, This indicates the displacement of the control point. The membrane strain matrix is... Let be the bending strain matrix. Considering the Rayleigh damping, integrating over the entire parameter domain and assembling, we obtain the discretized dynamic control equations as follows: In the formula, M is the global mass matrix, C is the Rayleigh damping term, K is the global stiffness matrix, and F is the global external load vector, which is: In the formula, and Here, A and D are the damping coefficients, and A and D are the membrane stiffness and bending stiffness matrices, respectively: Q can be calculated as The global mass and stiffness matrices can be assembled from the element mass and stiffness matrices, which are calculated using Gaussian point numerical integration, as follows: In the formula, The total number of units, Let J be the total number of Gaussian points within the element, and J be the mapping from the parameter domain to the physical domain. Jacobian matrix, The corresponding Gaussian point weights are used. In the proposed method, equivalent material interpolation for the dynamics of thin-shell structures is achieved based on SIMP, while equivalent elastic parameters and linear density interpolation are set to avoid matrix singularities and excessive mass penalty. In the formula, To minimize the elastic Young's modulus and prevent singularities, Young's modulus of solid materials The penalty index is generally taken as... To suppress intermediate density. For the density of solid materials, The minimum density value is generally taken as... To avoid zero mass in empty regions, the element mass and stiffness matrix can be further expressed as: (4) Consider an inviscid, incompressible, irrotational potential flow field whose velocity potential satisfies the Laplace equation. Applying Green's second identity, we can transform it into a boundary integral form: In the formula, For field points (including shell surface points). For the source point location, For free space Green's function, It is a Euclidean norm. These are geometric coefficients.

[0026] Introducing surface element source intensity function Its velocity potential at any point P on the boundary is In the formula, It is the boundary surface The source intensity is related to the location of the source point. Combining Neumann boundary conditions, the method proposed in this patent applies Gauss's formula to avoid directly calculating the singular integral when the source point coincides with the field point, thus obtaining the non-singular expression for the normal fluid velocity: In the formula, The outer surface of the shell, For the far field, The source point normal is given. The right-hand side of the above equation lacks a singular kernel, facilitating higher-order numerical integration. This equation can be solved using the panel method to obtain the source intensity. And further calculate its velocity potential energy.

[0027] The surface represented by T-splines is discretized into several panels, with source points defined on each panel. The source intensities can then be derived using boundary integral equations. This method calculates the integrals on the panels based on the Gaussian orthogonal method, as follows: In the formula, , , For the i-th collocation point, and These are the collocation point and the source point normal, respectively. Let g be the Jacobian matrix of the g-th point in the e-th unit. The weights are Gaussian. Assembled over all collocational points i, we obtain the weights for the source intensity. A system of linear equations with coefficients, where the unknowns are all Gaussian points. and matching points The source intensity is obtained by solving for the source intensity (since the Gaussian point is chosen as both the source point and the collocation point). The velocity potential can then be further solved. .

[0028] However, the passive function in the above velocity potential energy calculation is still singular. To address this issue, an equipotential auxiliary function is introduced (constructing the shell surface as an equipotential surface), which satisfies the homogeneous boundary integral equation as follows: In the formula, It is obtained through independent, singular iterative equations. This represents the distance from the source point to the selected origin, and it only appears as a scale reference to eliminate odd kernels. The panel method based on T-splines yields the velocity potential reconstruction formula as follows: In the formula, .

[0029] The velocity expression for any point P in the flow field is: Therefore, the pressure on the surface of the thin-shell structure can be calculated as follows: In the formula, For free-flow static pressure, For fluid density, For the free flow velocity, according to Bernoulli's equation, the pressure coefficient... It can be calculated as: Projecting the pressure along the normal to the shell surface into surface force density And map this onto the control points to obtain the expression for the aerodynamic load: (5) Applying geometric analysis methods such as linear elastic materials to solve the structural response of thin-shell structures, the equilibrium equation of the structural response can be expressed as: In the formula, Let be the excitation angular frequency, M, C, and K be the global mass, damping, and stiffness matrices, respectively, d be the reset displacement, and F be the global external load vector. Let Let the dynamic stiffness matrix be denoted as , then the above equation can be expressed as: Using the MQSRV model reduction strategy reduces the order while maintaining accuracy. Calculate the cost. This will be given... Target frequency band Dividing the interval into N sub-intervals, the center frequency of each interval is defined as follows: In the formula, each The representative frequency of this interval is used to generate a local quasi-static basis. This uniform division method can flexibly handle multi-peak frequency domain response characteristics, ensuring high-precision approximation within each neighborhood. At the nth center frequency... Define the frequency shift stiffness operator The basis derivation is as follows: To maintain the orthogonality of M, the basis is modified and normalized. The calculation is as follows: This yields a set of M orthonormal bases, which are... By concatenating the basis vectors obtained from each frequency band, we can obtain the global reduced-order basis of MQSRV, which is: In the formula, Let be the total reduced dimension. Let the approximate solution... Combined with the above formula and multiplied on the left The reduced dynamic control equations of MSQRV can be obtained as follows: And it is equivalent to the full-order system in the sense of energy inner product. Among them, and These represent the reduced-order mass, damping, stiffness matrix, and reduced-order external load vector, respectively. Reduced-order dynamic stiffness matrix. It can be deduced as The reduced-order response can then be solved as follows: according to The structural response can then be calculated.

[0030] A global density function is constructed based on Bézier elements to meet the requirements of topological description of thin-shell structures. Specifically, given a vector containing the initial Bézier control density... The smoothing mechanism improves the smoothness of the Bezier control density by using the Shepard function on the Bezier cells, as follows: In the formula, Represents the Shepard function at the Bessel control points. A smooth, continuous local density distribution function. Density can be controlled by smoothing. and Bessel unit basis functions Represented as: Topological boundary The isomorphic profile of a local DDF can be used. This can be represented as: The common control points between adjacent Bessel elements ensure that local DDFs can be connected into a topological description model to describe the overall structure. Therefore, local DDFs can be assembled into a global DDF through natural connections between adjacent elements to describe the entire thin-shell structure design domain, expressed as: In the formula, It is the number of Bézier elements in the design domain, symbol […]. Let represent the local DDF of the i-th Bessel element. The structural topology of the entire design domain. It can be composed of the local structural topology of each Bessel unit, as follows: (6) Substitute the design variables into the S3 analysis model and calculate the reduced-order dynamic stiffness matrix using the MSQRV reduced-order model. Then solve for the reduced-order response. And then according to Calculate the displacement field.

[0031] (7) Based on the model established in the above steps, the dynamic flexibility minimization topology optimization model is expressed as: In the formula, This refers to the initial node density of the Bézier control points; the design variable values ​​must be between their minimum values. The minimum value between 1 and 0 is to avoid singularities. This represents the total number of design variables. The objective function for calculating global dynamic compliance is given by: Calculate the size of a complex number. This refers to the excitation frequency range. This represents the global density distribution function. Consider an excitation frequency of... Under the action of a simple harmonic load, and These are the complex numbers representing the load and displacement, respectively. i is the imaginary unit. G is the volume constraint. Represents the volume fraction of a solid material. This represents the maximum value of material consumption. Represents a virtual displacement field, belonging to the kinematically acceptable space. . Indicates the Dirichlet boundary The specified displacement vector at the location. k, c, and m are the strain energy, damping energy, and mass bilinear energy, respectively: For a semi-linear load type l, the calculation formula is: In the formula, Represents body forces in thin-shell structures. This is the boundary traction force, belonging to the Neumann boundary. .

[0032] Based on the construction principle of the global density distribution function, the first derivative of the density distribution function with respect to the design variables can be obtained as follows: Therefore, dynamic flexibility The final expression for sensitivity is: The sensitivity of the volume constraint function is: Update the control density using the optimal criterion method. The updated design variables are substituted into the updated topology description model representing the global design domain. The updated model is used to update the displacement field, thereby updating the objective function and sensitivity calculations to obtain new control variables. It is then determined whether the convergence condition is met (the difference in node density between two consecutive iterations is less than a certain value). If not, the above iterative process continues. If it is met, the topology optimization configuration of the reinforced structure is obtained.

[0033] The following is combined with Figure 2 The following specific embodiment will be shown to illustrate the above-described results of the present invention: (1) Model the swept wing skin in Rhino software and convert the model into a T-spline geometry model, such as Figure 2 As shown in (a); (2) Import the model file. In order to clearly show the configuration of the reinforced structure, the design parameters are selected as shown in Table 1. Table 1. Optimization Design Parameters for Minimizing Dynamic Compliance of Swept Wings The model file is read, and a unified geometric representation model for the T-spline geometry is constructed using Bessel extraction techniques. Furthermore, geometric shell elements are constructed using Kirchhoff theory. The T-spline mixture function is then: In the formula, For the parameter coordinates in the Bessel parent unit domain, select the parameter coordinates in the interval [-1,1] for Gaussian integration, so as to facilitate subsequent analysis and optimization. Let be a set containing T-spline mixture functions, whose support points are located in Bessel element e, where a represents the relevant local index of the control point, and n is the number of control points. Extract the Bessel operator for element e. Let It is a vector of the T-spline mixture function of element e, and the T-spline mixture function of element e can be expressed as: In the formula, Combining Kirchhoff's theory with the construction of geometric shell elements, the deformation of a thin shell structure can be described by its mid-surface. Therefore, the deformation along the thickness direction on the mid-surface... The displacement vector of any point can be expressed as: In the formula, x and These are the position vectors of material points in deformed (actual) and undeformed (reference) states, respectively. and This represents the unit normal vector at a point in the midplane of the thin shell before and after deformation. The expressions for the components of the Green-Lagrange strain tensor can be further expressed as: In the formula, Represents mid-surface strain, which describes mid-surface stretching or compression. Representing the change in curvature, used to describe bending, is: The virtual work principle for thin-shell structures is derived using the Galerkin weak form, and the weak form is obtained through variational methods as follows: In the formula, f is the body force, and ρ is the material density. and These are the fourth-order tensors for membrane and bending stiffness, respectively.

[0034] (3) Discretize the surface of the T-spline geometry into several panels, and define source points on each panel. Then, solve the source intensity by solving the boundary integral equation. Calculate the integral on the panel based on the Gaussian orthogonal method; its normal velocity is: In the formula, , , For the i-th collocation point, and These are the collocation point and the source point normal, respectively. Let g be the Jacobian matrix of the g-th point in the e-th unit. The weights are Gaussian. Assembled over all collocational points i, we obtain the weights for the source intensity. A system of linear equations with coefficients, where the unknowns are all Gaussian points. and matching points The source intensity is obtained by solving for the source intensity (since the Gaussian point is chosen as both the source point and the collocation point). The velocity potential can then be further solved. .

[0035] Using nonsingular equations based on T-splines, the velocity potential at the body surface can be calculated more accurately and efficiently: In the formula, .

[0036] The velocity at any point P in the flow field can be calculated from the gradient of the velocity potential, as follows: However, when point P is located on the body surface, the integral exhibits singularity. Therefore, a T-spline mixture function based on Bessel extraction is used to describe the velocity potential at the body surface to address the singular integral in the boundary integral calculation. The higher-order continuity of the Bernstein polynomial allows for direct calculation of the derivative, thereby improving the accuracy of the velocity potential gradient calculation. The velocity potential at a point on the body surface is constructed as follows: In the formula, The Bernstein bivariate basis functions can be calculated using the tensor product of the Bernstein univariate basis functions, and the calculation formula is as follows: In the formula, p defines the order of the polynomial. This represents the velocity potential energy at the control point. The velocity U of the flow field can be directly obtained as: Therefore, the pressure on the surface of the thin-shell structure can be calculated as follows: In the formula, For free-flow static pressure, For fluid density, For the free flow velocity, according to Bernoulli's equation, the pressure coefficient... It can be calculated as: Projecting the pressure along the normal to the shell surface into surface force density And map this onto the control points to obtain the expression for the aerodynamic load: Aerodynamic loads of three airfoil sections of a swept wing, such as Figure 2 As shown in (b).

[0037] The aerodynamic load can be directly applied to the Bessel control points on the mid-surface of the thin shell. The structural response of the thin shell structure can be solved using geometric analysis methods such as linear elastic materials combined with the MSQRV reduced-order model. The MSQRV reduced-order dynamic control equations are: And it is equivalent to the full-order system in the sense of energy inner product. Among them, and These represent the reduced-order mass, damping, stiffness matrix, and reduced-order external load vector, respectively. Reduced-order dynamic stiffness matrix. It can be deduced as The reduced-order response can then be solved as follows: according to The structural response can then be calculated.

[0038] (4) Construct a global density function based on Bézier elements to meet the requirements of the topological description of the thin shell structure. Specifically, given a vector containing the initial Bézier control density The smoothing mechanism improves the smoothness of the Bezier control density by using the Shepard function on the Bezier cells, as follows: In the formula, Represents the Shepard function at the Bessel control points. A smooth, continuous local density distribution function. Density can be controlled by smoothing. and Bessel unit basis functions Represented as: Topological boundary The isomorphic profile of a local DDF can be used. This can be represented as: The common control points between adjacent Bessel elements ensure that local DDFs can be connected into a topological description model to describe the overall structure. Therefore, local DDFs can be assembled into a global DDF through natural connections between adjacent elements to describe the entire thin-shell structure design domain, expressed as: In the formula, It is the number of Bézier elements in the design domain, symbol […]. Let represent the local DDF of the i-th Bessel element. The structural topology of the entire design domain. It can be composed of the local structural topology of each Bessel unit, as follows: (5) In isogeometric analysis, the element Hexu The displacement field is described by interpolation using the same set of basis functions: In the formula, The total number of basis functions, , To control the displacement of the control point. According to Kirchhoff's theory, the membrane strain and curvature change can be further expressed as: In the formula, This indicates the displacement of the control point. The membrane strain matrix is... Let be the bending strain matrix. Considering the Rayleigh damping, integrating over the entire parameter domain and assembling, we obtain the discretized dynamic control equations as follows: In the formula, M is the global mass matrix, C is the Rayleigh damping term, K is the global stiffness matrix, and F is the global external load vector, which is: In the formula, and Here, A and D are the damping coefficients, and A and D are the membrane stiffness and bending stiffness matrices, respectively: Q can be calculated as: The global mass and stiffness matrices can be assembled from the element mass and stiffness matrices, which are calculated using Gaussian point numerical integration, as follows: In the formula, The total number of units, Let J be the total number of Gaussian points within the element, and J be the mapping from the parameter domain to the physical domain. Jacobian matrix, The corresponding Gaussian point weights are used. In the proposed method, equivalent material interpolation for the dynamics of thin-shell structures is achieved based on SIMP, while equivalent elastic parameters and linear density interpolation are set to avoid matrix singularities and excessive mass penalty. In the formula, To minimize the elastic Young's modulus and prevent singularities, Young's modulus of solid materials The penalty index is generally taken as... To suppress intermediate density. For the density of solid materials, The minimum density value is generally taken as... To avoid zero mass in empty regions, the element mass and stiffness matrix can be further expressed as: Subsequently, the displacement field of the thin-shell structure was calculated based on the reduced-order system.

[0039] (6) The topology optimization model that minimizes dynamic flexibility is expressed as: In the formula, This refers to the initial node density of the Bézier control points; the design variable values ​​must be between their minimum values. The minimum value between 1 and 0 is to avoid singularities. This represents the total number of design variables. The objective function for calculating global dynamic compliance is given by: Calculate the size of a complex number. This refers to the excitation frequency range. This represents the global density distribution function. Consider an excitation frequency of... Under the action of a simple harmonic load, and These are the complex numbers representing the load and displacement, respectively. i is the imaginary unit. G is the volume constraint. Represents the volume fraction of a solid material. This represents the maximum value of material consumption. Represents a virtual displacement field, belonging to the kinematically acceptable space. . Indicates the Dirichlet boundary The specified displacement vector at the location. k, c, and m are the strain energy, damping energy, and mass bilinear energy, respectively: For a semi-linear load type l, the calculation formula is: In the formula, Represents body forces in thin-shell structures. This is the boundary traction force, belonging to the Neumann boundary. .

[0040] Based on the construction principle of the global density distribution function, the first derivative of the density distribution function with respect to the design variables can be obtained as follows: Therefore, dynamic flexibility The final expression for sensitivity is: The sensitivity of the volume constraint function is: (8) Update the control density (design variables) using the OC method. Substitute the updated design variables into the updated topology description model representing the complex structural design domain. Update the displacement field using the updated model, thereby updating the objective function and sensitivity calculations to obtain new control variables. Determine whether the convergence condition is met. If not, continue the above iterative process. If it is met, the optimized topology configuration is obtained. The convergence criterion is that the difference in node density between two consecutive iterations is less than 1%, or the maximum iteration step size of 400 is reached.

[0041] The final structural topology of the swept wing is as follows Figure 4 As shown in (a), a continuous material distribution structure matching the aerodynamic load distribution is formed. The optimized configuration forms multiple truss structures interwoven along the spanwise and chordal directions on the upper and lower surfaces, with local holes distributed in a trapezoidal pattern. This topological feature effectively adjusts the bending-torsional coupling stiffness of the airfoil, causing the local strain energy density to redistribute along the mainstream aerodynamic direction, significantly weakening the resonance response of low-order modes. In addition, the convergence history of the objective function and constraint function at 0 angle of attack is as follows: Figure 3 As shown, the red and blue curves represent the evolution of the objective function and volume fraction, respectively. The figure also illustrates the optimized configurations under different iteration steps; by the 60th iteration, the characteristic structure of the topology is clearly visible. With changes in structural morphology and a gradual decrease in volume, the objective function converges smoothly and rapidly, and the volume fraction reaches the specified maximum consumption value. The final optimized value for the dynamic flexibility of the structure is 903.90. The iteration curves all converge smoothly and rapidly, demonstrating the effectiveness and efficiency of the proposed method for engineering thin-shell structures in aerodynamic environments.

[0042] Figure 5A comparison of the frequency response function curves before and after optimization is presented, clearly reflecting the dynamic differences. After optimization, the first-order natural frequency is significantly increased from 42Hz to 103Hz, and the strong single-peak resonance in the target frequency band of the original design is suppressed. The optimization results show a clear anti-resonance region near 100Hz, indicating that the vibration energy of the structure is redistributed among multiple modes, and the overall frequency response curve tends to be flatter. A comparison of the FRF curves before and after optimization considering Rayleigh damping shows that the damping effect further smooths the high-frequency response, the peak response of the optimized structure decreases by approximately 40%, and a continuous attenuation trend is observed in the high-frequency range. The swept-wing skin case fully validates the application potential of the frequency response-driven isogeometric topology optimization method in complex aerodynamic thin shells. The optimized skin structure, while maintaining a lightweight design with smooth boundaries, achieves low-frequency response suppression and high-frequency response attenuation, providing a new design approach for vibration suppression of aerospace thin shell structures under aerodynamic loads.

[0043] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A dynamic response-driven isogeometric topology optimization method, characterized in that, include: S1. Model the geometric model of a complex thin-shell structure in CAD software based on T-splines, and export the model; S2. Construct a unified geometric representation model for T-spline geometric models based on Bézier extraction technology; S3. Construct isogeometric thin-shell elements based on Kirchhoff theory to address structural frequency response, and establish an analytical model for the thin-shell structure. S4. The aerodynamic response of the structure surface is solved by the high-order panel method based on T-splines. The collocation method is used to satisfy the Neumann condition required for aerodynamics. The source density distribution of the structure surface is solved based on the T-spline mixing function. The velocity potential of the structure surface is described by the Bessel element basis function, and the pressure distribution of the thin shell surface is obtained. S5. Apply the aerodynamic pressure distribution as a load directly to the Bessel control points, construct a local density distribution function based on the Bessel elements, and assemble it into a global density distribution function. S6. Solve the element stiffness and mass matrices of the uniformly geometric shell elements using the Gaussian integral method and assemble them into global stiffness and mass matrices. Solve the displacement field of the thin shell under aerodynamic load based on the reduced-order model. S7. Based on the aerodynamic response field, establish an iso-geometric topology optimization model, calculate the objective function and derive the complete analytical sensitivity, iteratively update the design variables based on the optimal criterion method, and finally obtain the optimized configuration after iteration.

2. The isogeometric topology optimization method based on dynamic response-driven method according to claim 1, characterized in that, The implementation of step S1 includes: In Rhinoceros 3D software, the geometric model of any thin shell structure can be converted into a T-spline-based surface geometry model by combining it with the Autodesk T-splines plugin. The integration of this plugin with Rhinoceros 3D makes the geometric modeling process more visual.

3. The isogeometric topology optimization method based on dynamic response driving according to claim 2, characterized in that, The implementation of step S2 includes: The mid-surface geometry of an arbitrary thin-shell structure can be represented using a T-spline mixture function as follows: In the formula, A represents the vertex number in the T mesh. and It is the local node vector that defines the mixing function, which can be obtained through the Cox-de Boor recursive formula; When p = 0, When p≥1 In the formula, It is the i-th node; Based on the Bessel extraction technique, a unified geometric representation model can be constructed for any thin shell structure. Therefore, the T-spline mixture function can be further expressed as: In the formula, Given the parameter coordinates in the Bessel parent unit domain, select the parameter coordinates in the interval [-1, 1] for Gaussian integration; Let be a set containing T-spline mixture functions, whose support points are located in Bessel element e, where a represents the relevant local index of the control point, and n is the number of control points; To extract the Bessel operator for element e, assuming the polynomial order is the same in all directions, the extraction matrix can be obtained through the following steps: Establish a one-dimensional extraction matrix of order p and q in each direction. ; Combining in the form of tensor products: ; For special T-spline cases, if the local junction vectors are irregular, then after expanding the T-junction, the effective support region should be taken, and the corresponding structure should be constructed. ; To obtain the Bessel expression for the T-spline blending function, a cell mapping from the parent cell domain to the physical domain can be defined, described as follows: In the formula, Describes a rational T-spline mixture function of a Bessel expression; This represents a control point matrix containing unit control points; Indicates the unit weight; Representing the unit weights in diagonal matrix form, the T-spline mixture function can be further expressed as: The control points of the global T-spline can be mapped to the control points of the Bézier element by the transpose of the extraction operator. Therefore, the Bézier control points and weights can be represented by the element extraction operator, and the T-spline curve corresponding to the control points and weights is: In the formula, Unit weight The T-spline element is in diagonal matrix form. This means that each T-spline element has a corresponding equivalent Bessel element, and the further expression for the T-spline mixing function is: In the formula, The formulas for calculating the first and second derivatives of the T-spline mixture function with respect to the local parameter coordinates are as follows: as well as The derivatives related to physical coordinates can be calculated using the chain rule.

4. The isogeometric topology optimization method based on dynamic response driving according to claim 3, characterized in that, The implementation of step S3 includes: Based on Kirchhoff's thin shell theory, the geometric shell element is constructed, and its deformation is represented by the mid-face of the shell. Therefore, along the thickness direction on the mid-face... The displacement vector of any point can be expressed as: In the formula, x and These are the position vectors of material points in deformed and undeformed states, respectively. and The unit normal vectors at points in the mid-plane of the thin shell before and after deformation are represented by the following basis vectors: Where the upper and lower scales of the second-order tensor represent the inverse covariance components, assuming thickness... Then, higher-order terms can be ignored. To simplify, it can be represented as: Similarly, this also applies to deformable configurations; and These are the first and second basic forms, respectively, represented as: The expressions for the components of the Green-Lagrange strain tensor can then be further expressed as: In the formula, Represents mid-surface strain, which describes mid-surface stretching or compression. Representing the change in curvature, used to describe bending, is: The virtual work principle for thin-shell structures is derived using the Galerkin weak form, and the weak form is obtained through variational methods as follows: In the formula, f represents body force, and ρ represents material density. and These are the fourth-order tensors for membrane stiffness and bending stiffness, respectively, which can be derived from the fourth-order tensor of the material. get: In isogeometric analysis, the element displacement field and virtual displacement field Interpolation described using the same set of basis functions In the formula, The total number of basis functions, , To control the displacement of the control point; according to Kirchhoff's theory, the membrane strain and curvature change can be further expressed as: In the formula, Indicates the displacement of the control point. The membrane strain matrix is... The bending strain matrix is: Considering the Rayleigh damping integral over the entire parameter domain and assembled, the discretized dynamic control equations are obtained as follows: In the formula, M is the global mass matrix, C is the Rayleigh damping term, K is the global stiffness matrix, and F is the global external load vector, which is: In the formula, and This is the damping coefficient, and A and D are the membrane stiffness and bending stiffness matrices, respectively. Q can be calculated as: The global mass and stiffness matrices can be assembled from the element mass and stiffness matrices, which are calculated using Gaussian point numerical integration, as follows: In the formula, The total number of units, Let J be the total number of Gaussian points within the element, and J be the mapping from the parameter domain to the physical domain. Jacobian matrix, For the corresponding Gaussian point weights, the proposed method uses SIMP to achieve dynamic equivalent material interpolation of thin-shell structures. Simultaneously, it sets equivalent elastic parameters and linear density interpolation to avoid matrix singularities and excessive mass penalty, as follows: In the formula, To minimize the elastic Young's modulus and prevent singularities, Young's modulus of solid materials The penalty index is generally taken as... To suppress intermediate density, For the density of solid materials, The minimum density value is generally taken as... To avoid zero mass in the empty region; based on the above equation, the element mass and stiffness matrix can be further expressed as: The high-order continuity based on T-spline geometry enables the smoothing and accuracy of the structural dynamic response field.

5. The isogeometric topology optimization method based on dynamic response driving according to claim 4, characterized in that, Step S4 consists of two parts: S4.1 Consider a non-viscous, incompressible, irrotational potential flow field with velocity potential. Satisfying the Laplace equation, it is: In the formula, For the external fluid domain; applying Green's second identity, the velocity potential is transformed into a boundary integral form, as follows: In the formula, For field points (including shell surface points). The source point location, For free space Green's function, It is a Euclidean norm. Geometric coefficients; Introducing surface element source intensity function Its velocity potential at any point P on the boundary is: In the formula, It is the boundary surface The source intensity is related to the location of the source point. Combining the Neumann boundary conditions, the method proposed in this patent applies the Gaussian formula to avoid directly calculating the singular integral when the source point and field point coincide, thus obtaining the non-singular expression for the normal fluid velocity: In the formula, The outer surface of the shell, For the far field, The source point normal is given. The right-hand side of the above equation has no singular kernel, facilitating higher-order numerical integration; this equation can be solved using the panel method to obtain the source intensity. And further calculate its velocity potential energy; The proposed method describes the geometry of thin-shell structures based on T-splines. Therefore, the generalized expression for any point P on the spline surface is defined as: Furthermore, the unit normal vector of point P on the T-spline surface can be obtained. : In the formula, and Let denote the derivative of the surface with respect to the directions ξ and η; for ease of subsequent derivation, the expression for the symbol is defined as... The calculation formula is as follows: For the source points on the body surface and If the source point is on the boundary surface, then the normal velocity can be converted into an expression based on the T-spline description, as follows: The body surface, represented by T-splines, is discretized into several panels, with source points defined on each panel. The source intensities can then be derived using boundary integral equations. This method calculates the integrals on the panels based on the Gaussian orthogonal method, as follows: In the formula, , , For the i-th collocation point, and These are the collocation point and the source point normal, respectively. Let g be the Jacobian matrix of point g on the e-th unit. The weights are Gaussian; the assembly is performed over all collocation points i, thus obtaining the weights for the source strength. A system of linear equations with coefficients, where the unknowns are all Gaussian points. and matching points The source strength is obtained by solving for the source strength. The velocity potential can then be further solved. ; Introducing an equipotential auxiliary function, the homogeneous boundary integral equation it satisfies is: In the formula, It is obtained through independent, singular iterative equations. This represents the distance from the source point to the selected origin, which only appears as a scale reference to eliminate odd kernels; the panel method based on T-splines yields the velocity potential reconstruction formula as follows: In the formula, . The velocity expression for any point P in the flow field is: This equation solves for the velocity distribution of the fluid on the surface of the thin-shell structure, including the velocity components in all directions; after obtaining the velocity potential energy, its velocity can be further calculated. S4.

2. The velocity potential at the volume surface is described using the T-spline mixture function extracted by Bessel, in order to solve the singular integral in the boundary integral calculation; the higher-order continuity of the Bernstein polynomial allows for the direct calculation of the derivative, and the velocity potential at a point on the volume surface is constructed as follows: In the formula, The Bernstein bivariate basis functions can be calculated using the tensor product of the Bernstein univariate basis functions, and the calculation formula is as follows: In the formula, p defines the order of the polynomial. This represents the velocity potential energy at the control point. The velocity U of the flow field can be directly obtained as: Therefore, the pressure on the surface of the thin-shell structure can be calculated as follows: In the formula, For free-flow static pressure, For fluid density, For the free flow velocity, according to Bernoulli's equation, the pressure coefficient... It can be calculated as: Projecting the pressure along the normal to the shell surface into surface force density And map this onto the control points to obtain the expression for the aerodynamic load: For a thin shell with uniform geometry, the dynamic equation can be expressed as: In the formula, This is the aerodynamic load vector. For other mechanical loads, for low-speed, lift-dominated problems with relatively small geometric deflections, a unidirectional weakly coupled approximation can be used, where the aerodynamic forces are evaluated with the initial configuration and treated as constants within the iterations: In the formula, the global stiffness matrix K is updated with the design variables, while It does not update with displacement feedback.

6. The isogeometric topology optimization method based on dynamic response driving according to claim 5, characterized in that, Step S5 consists of two parts: S5.

1. Applying geometric analysis methods such as linear elastic materials to solve the structural response of thin-shell structures, the equilibrium equation of the structural response can be expressed as: In the formula, Let be the excitation angular frequency, M, C, K be the global mass, damping, and stiffness matrices respectively, d be the reset displacement, and F be the global external load vector. Let be... Let be the dynamic stiffness matrix, then the above equation can be expressed as: In frequency domain dynamics optimization, when design variables (such as topology density) are iteratively updated, the solution is directly repeated on the fully free system. The computational cost is extremely high; the MQSRV model order reduction strategy is adopted to reduce the computational cost while maintaining accuracy. Calculation cost will give Target frequency band Dividing the interval into N sub-intervals, the center frequency of each interval is defined as follows: In the formula, each The representative frequency of this interval is used to generate a local quasi-static basis; at the nth center frequency Define the frequency shift stiffness operator The basis derivation is as follows: The recursive process of the basis involves linear equations To avoid high-frequency local mode distortion, and to maintain the orthogonality of M, the basis is modified and normalized. Thus, the MQSRV basis... The calculation is as follows: This yields a set of M orthonormal bases, namely: By concatenating the basis vectors obtained from each frequency band, we can obtain the global reduced-order basis of MQSRV, which is: In the formula, Let the total reduced-order dimension be the approximate solution. Combined with the above formula and multiplied on the left The reduced dynamic control equations of MSQRV can be obtained, which significantly reduce the dimensionality of the original system: And it is equivalent to the full-order system in the sense of energy inner product. Among them, and Let the reduced-order mass, damping, stiffness matrix, and reduced-order external load vector be respectively, which can be derived as follows: Reduced-order dynamic stiffness matrix It can be deduced as: The reduced-order response is then solved as follows: according to The structural response can then be solved; S5.2 Construct a global density function based on Bézier elements to meet the requirements of topological description of thin shell structures. Specifically, given a vector containing the initial Bézier control density... The smoothing mechanism improves the smoothness of the Bezier control density by using the Shepard function on the Bezier cells, as follows: In the formula, The Shepard function represents the Bessel control point; the smooth, continuous local density distribution function. Density can be controlled by smoothing. and Bessel unit basis functions Represented as: Topological boundary The isomorphic profile of a local DDF can be used. This can be represented as: The common control points between adjacent Bessel elements ensure that local DDFs can be connected into a topological description model to describe the overall structure; therefore, local DDFs can be assembled into a global DDF through the natural connections between adjacent elements to describe the entire thin-shell structure design domain, and its expression is: In the formula, It is the number of Bézier elements in the design domain, symbol […]. Represents the local DDF of the i-th Bessel element; the structural topology of the entire design domain. It can be composed of the local structural topology of each Bessel unit, as follows: 。 7. The isogeometric topology optimization method based on dynamic response driving according to claim 6, characterized in that, The implementation methods of S6 include: Substitute the design variables into the S3 analysis model and calculate the reduced-order dynamic stiffness matrix using the MSQRV reduced-order model. Then solve for the reduced-order response. And then according to Calculate the displacement field.

8. The isogeometric topology optimization method based on dynamic response driving according to claim 7, characterized in that, The implementation methods of S7 include: S7.1 The dynamic flexibility minimization topology optimization model is expressed as: In the formula, This refers to the initial node density of the Bézier control points; the design variable values ​​must be between their minimum values. The minimum value between 1 and 0 is to avoid singularities; This represents the total number of design variables. The objective function for calculating global dynamic compliance is, where Calculate the size of a complex number; The range of excitation frequencies; Represent the global density distribution function; consider the excitation frequency as... Under the action of a simple harmonic load, and These are the complex forms of load and displacement, respectively; i is the imaginary unit; G is the volume constraint. Represents the volume fraction of a solid material. This represents the maximum value of material consumption. Represents a virtual displacement field, belonging to the kinematically acceptable space. ; Indicates Dirichlet boundary The specified displacement vector at the location; k, c, and m are the strain energy, damping energy, and mass bilinear energy, respectively: For a semi-linear load type l, the calculation formula is: In the formula, Represents body forces in thin-shell structures. This is the boundary traction force, belonging to the Neumann boundary. ; S7.

2. The sensitivity analysis of the optimization model requires calculating the first derivative of the objective function relative to the density distribution function. Based on the MSQRV order reduction strategy, the objective function can be approximated by the frequency band numerical integral form: objective function If a variable is composed of a real part and an imaginary part, then its absolute value is defined as: In the formula, and Let these represent the real and imaginary parts of the functional, respectively; the first derivative of the objective function with respect to the design variables and the dynamic compliance. related: The sensitivity analysis of the objective function is transformed into the derivation of dynamic compliance. The first derivative with respect to the design variables; dynamic compliance is used to derive the calculation results. It can be represented as: Dynamic flexibility The first derivative of the relative global density distribution function can be expressed as: In the formula, It is the derivative of the displacement field with respect to the global DDF. It is the derivative of the virtual displacement field with respect to the DDF; Combining the optimization model with the dynamic response analysis of the thin-shell structure, we can conclude that: The first derivative of the bilinear expression for strain energy k and mass m with respect to the global density distribution function is: because Then the following equations hold: Considering that the problem of minimizing dynamic compliance is self-adjoint, the virtual displacement field in the above equation can be eliminated: Then the flexibility The first derivative of the relative global density distribution function is: Further calculation of the first derivatives of the membrane stiffness A and bending stiffness D matrices with respect to the global DDF yields the dynamic compliance. The expression for the derivative of the global density distribution function is: Based on the construction principle of the global density distribution function, the first derivative of the density distribution function with respect to the design variables can be obtained as follows: Therefore, dynamic flexibility The final expression for sensitivity is: The sensitivity of the volume constraint function is: S7.3 Update the control density using the optimal criterion method. That is, design variables, repeat the process until the iteration termination condition is met, and obtain the enhanced structure topology optimization configuration based on the objective function and sensitivity calculated in the last iteration S7.

2.

9. The application of a dynamic response-driven isogeometric topology optimization method, applied to any one of claims 1 to 8, characterized in that: It is applied to the design of swept wing skin.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium includes a stored computer program, wherein the computer program, when executed by a processor, controls the device containing the storage medium to perform a dynamic response-driven isogeometric topology optimization method applied to any one of claims 1 to 8.

Images