A topology optimization method based on aerodynamic driven particle flow for flexible airfoils

Abstract

The application discloses a kind of topological optimization methods for flexible airfoil based on aerodynamic driven particle flow, it is related to aircraft structure design technical field, including: S1, based on NURBS model construction and initialization, generate uniformly distributed particles in parameter space, and it is distributed on the wing geometry model in physical space and completes initialization;Definition design parameter, including material attribute, structure size, aerodynamic parameter and boundary condition;Using smooth particle hydrodynamics kernel function, particle information is projected and mapped to control point and gauss point, establish the link between design variable and geometric control point.This method not only greatly enhances the accuracy and stability of flexible wing large deformation structure topology reconstruction, but also can directly output absolute smooth sawtooth-free CAD continuous boundary, significantly improve the safety, reliability and manufacturability of advanced slender aircraft structure design.

Description

Technical Field

[0001] This invention relates to the field of aircraft structural design technology, specifically to a topology optimization method for flexible airfoils based on aerodynamically driven particle flow. Background Technology

[0002] Currently, modern high-performance, long-endurance aircraft typically employ slender wing structures to achieve optimal aerodynamic efficiency. However, this significant flexible structural characteristic makes them highly susceptible to large geometrical nonlinear deformations under aerodynamic loads. Unlike the small deformation assumptions of traditional rigid wings, these large deformations not only alter the wing's aerodynamic shape and angle-of-attack distribution in real time, triggering complex nonlinear fluid-structure interaction effects, but also cause aerodynamic loads to exhibit significant "dynamic dependence," rendering traditional linear theory-based analysis methods completely ineffective. With increasing demands for structural safety, mechanical stability, and buckling / flutter resistance under extreme aircraft conditions, a systematic numerical-driven topology optimization method is urgently needed to accurately capture these geometrical nonlinear behaviors and achieve optimal aeroelastic design of the wing structure.

[0003] Topology optimization provides a systematic design approach for addressing nonlinear fluid-structure interaction challenges and improving wing performance in modern aircraft design. Currently, topology optimization methods are primarily evolving from implicit to explicit and multi-scale approaches. Mesh-based methods, such as the Variable Density Method (SIMP), perform aerodynamic-structural co-optimization by distributing material density across a fixed mesh. Geometry-based topology optimization methods, such as the Level Set Method (LSM) and Moving Deformable Component (MMC), flexibly update the design domain boundary through parametric geometry, overcoming the shortcomings of gray elements and serrated boundaries in SIMP, and demonstrating excellent boundary description capabilities when dealing with large deformation nonlinear problems in flexible wings. However, these traditional finite element method (FEM) topology optimization methods often face the inherent contradiction of separating the geometric model from the analysis model when applied to slender wings. Especially under strong geometric nonlinearity and complex aerodynamic-structural coupling, low-density regions are prone to severe mesh distortion or even element flipping due to large loads, leading to non-physical divergence and iterative oscillations in the finite element solution.

[0004] To address the numerical computation challenges posed by nonlinear large deformations and mesh distortion, various methods have been proposed. For example, element connectivity parameterization methods based on elastic connections and strategies based on second-order adjoint differentials have alleviated the ill-conditioned stiffness matrix and stability degradation problems to some extent. Meanwhile, to eliminate the mesh dependency between the fluid and solid domains, researchers have deeply integrated non-uniform rational B-splines (NURBS) with the aerodynamic surface element method, constructing an explicit coupling paradigm that greatly simplifies the derivation of the adjoint equation and improves the sensitivity and efficiency of the solution. Nevertheless, these methods only introduce improvements on the aerodynamic surface and cannot fundamentally eliminate the numerical limitations imposed on the internal topological evolution of the structure. Under extreme large deformation conditions, limitations such as high computational cost, mesh dependency, and artificial checkerboard patterns still exist, and their theoretical framework requires further refinement.

[0005] In recent years, the application of isogeometric analysis (IGA) methods in topology optimization has gradually gained attention. It directly uses NURBS basis functions from CAD systems as shape functions, achieving mathematical unification of geometric modeling and structural analysis, and demonstrating significant advantages in generating absolutely smooth boundaries and eliminating mesh dependencies. To further overcome the limitations of fixed parameter spaces on large deformation topology reconstruction, a particle flow topology optimization (PFTO) framework based on isogeometric material point method (I-MPM) is proposed, deeply integrating IGA with meshless particle methods such as the material point method (MPM). This method effectively combines the high-fidelity geometric characteristics of IGA with the numerical robustness of particle methods, fundamentally overcoming the solution divergence problem caused by mesh distortion under extreme nonlinear deformation. Nevertheless, the application of the PFTO method in flexible airfoils is still in the exploratory stage, and its theoretical framework and numerical implementation require further research and improvement. Summary of the Invention

[0006] To address the aforementioned technical challenges, a topology optimization method based on aerodynamically driven particle flow is provided for flexible airfoils. By introducing a movable particle flow into the design domain within the airfoil, the design variables of traditional fixed control points are transformed into a dynamic evolution problem of particle position and material properties, enabling continuous and controllable adjustment of the material distribution of the load-bearing skeleton inside the airfoil. The Lagrange-Euler hybrid description (matter point method) is used to simulate particle motion, which greatly enhances the response capability to extreme nonlinear topological changes in the internal structure while maintaining high-fidelity geometric continuity of the airfoil's outer surface.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A topology optimization method based on aerodynamically driven particle flow for flexible airfoils, comprising: S1. NURBS-based model construction and initialization: Generate uniformly distributed particles in the parameter space and distribute them on the wing geometric model in the physical space to complete the initialization; Define design parameters, including material properties, structural dimensions, aerodynamic parameters and boundary conditions; Utilize the smooth particle hydrodynamic kernel function to project and map particle information to control points and Gaussian points, and establish the relationship between design variables and geometric control points. S2. Aerodynamic load calculation: The source panel method based on NURBS is adopted. The numerical integration method based on NURBS uses NURBS basis functions to describe the surface velocity potential, calculates the aerodynamic load on the NURBS control points, and applies the corresponding boundary conditions and loads. The calculated aerodynamic load is introduced into the subsequent structural analysis as a constant external load. S3. Based on the nonlinear solution of the uniformly geometric particle flow, perform geometric nonlinear finite element analysis, calculate the tangent stiffness matrix and residual force for iterative solution, and use the Newton-Raphson iterative method to update the displacement field in the design domain to capture the nonlinear response of the structure. S4. Calculation and optimization update of the objective function: Calculate the objective function and perform sensitivity analysis. Use the KD tree to search for the nearest neighbor region and transfer the sensitivity information of the control points back to the particles. Use the moving asymptote method to update the design variables and determine whether the convergence condition is met. If it is met, proceed to the next step; otherwise, return to step S3 to enter the next iteration. S5. Construction of density distribution function: Introduce density distribution function and use relaxation factor to map density information at control points to the entire design domain, transform from discrete control points to continuous density field, describe the material distribution topology of the wing's internal structure, and output the final optimization result.

[0008] Preferably, the specific operations for NURBS-based model construction and initialization in step S1 are as follows: Uniformly distributed particles are generated in the parameter space and distributed onto the wing geometry model in physical space through mapping; design parameters including material properties, structural dimensions, aerodynamic parameters and boundary conditions are defined. A smooth particle hydrodynamic kernel function is used to determine the neighborhood with the grid point as the center. The relative distance between the particle position in the neighborhood and the center point is calculated and the kernel function is called to calculate the weight. The material information of the particle is projected and mapped to the control point and Gaussian point respectively.

[0009] Preferably, the specific operation for calculating the aerodynamic load in step S2 is as follows: The aerodynamic field is solved using the NURBS-based source panel method; In the calculation process, the numerical integration method based on NURBS is used to describe the velocity potential of the wing surface using NURBS basis functions; Solve the surface element equations to calculate the aerodynamic load distribution at the NURBS control points and apply the corresponding aerodynamic boundary conditions; the extracted aerodynamic loads will be directly introduced as constant external excitation loads into the geometric nonlinear structural analysis in step S3.

[0010] Preferably, the specific process of nonlinear solution based on isogeometric particle flow in step S3 is as follows: For the geometric nonlinear large deformation characteristics of airfoils, structural analysis is performed based on the isogeometric finite element analysis framework; the current tangential stiffness matrix of the structure and the residual force between internal forces and external aerodynamic loads are calculated; based on the Newton-Raphson iteration method, the displacement field in the design domain is continuously solved and updated by the residual force until the equations meet the convergence accuracy, thus capturing the nonlinear mechanical response of flexible airfoils under large deformation.

[0011] Preferably, the specific process of calculating the objective function and transferring sensitivity in step S4 is as follows: The objective function for topology optimization based on nonlinear displacement field calculation and its sensitivity to Gaussian point density are introduced. A sensitivity analysis transfer mechanism is introduced, and the KD tree algorithm is used to search the neighboring particle region around the control point. Based on the chain rule, the sensitivity information of the control point is transferred back to each particle to obtain the sensitivity of the objective function to the design variables of the bottom particle.

[0012] Preferably, the specific process of optimization update and convergence judgment in step S4 is as follows: Based on the acquired particle sensitivity, the position and attribute design variables of the particles are iteratively updated using the moving asymptote method; the rate of change of the objective function value is monitored and the volume constraint is detected, and the system determines whether the set convergence condition is met: if the convergence condition is met, the optimization is determined to end; otherwise, based on the updated particle flow state, the system returns to step S3 to perform the next nonlinear calculation and optimization iteration.

[0013] Preferably, the construction process of the density distribution function in step S5 is as follows: Density information at the converged control points is extracted, and a density distribution function is introduced. Using a set relaxation factor, the density information at the discrete control points is smoothly mapped to the entire airfoil design domain. By utilizing the high-order analytical properties of NURBS basis functions, a seamless transformation from the discrete control point set to the continuous physical density field is achieved, describing and extracting the high-fidelity topological boundary of the material distribution of the internal load-bearing skeleton of the wing.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention constructs an aerodynamically driven particle flow topology optimization system by using movable particles representing the material distribution in the physical domain of a flexible airfoil as design variables. It fundamentally overcomes the solution divergence problem caused by mesh distortion under extreme large deformations by employing a NURBS-based source panel method and Newton-Raphson nonlinear iterative solution. Forward mapping of control points and Gaussian points is achieved through a smooth particle hydrodynamic kernel function, and the sensitivity of ions is reversed using the KD tree algorithm, ensuring high-fidelity information interaction between particles, geometric meshes, and the physical field. Finally, a density distribution function is introduced to construct a continuous density field, and combined with the moving asymptote method optimization algorithm, accurate conversion from discrete particles to continuous topology is achieved. This method not only greatly enhances the accuracy and stability of topology reconstruction for large deformation flexible airfoil structures but also directly outputs absolutely smooth, serrated CAD continuous boundaries, significantly improving the safety, reliability, and manufacturability of advanced slender aircraft structure design. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the optimized method flow of the present invention; Figure 2 A schematic diagram of an undeformed and deformed configuration considering geometric nonlinearity is provided for this invention; Figure 3 This invention provides a schematic diagram of a nonlinear analysis of a classic cantilever beam using a particle flow topology optimization method. Figure 4 This invention provides a schematic diagram of a nonlinear analysis of an airfoil using a particle flow topology optimization method. Figure 5 A schematic diagram of NACA 2415 wing geometry parameter setting and surface aerodynamic load distribution provided by the present invention; Figure 6 A schematic diagram illustrating an iterative process for the topological distribution of all particles, control points, and structures within a design domain, provided by this invention; Figure 7 A schematic diagram of the evolution of wing rib structure considering geometric nonlinear effects (taking the tenth rib as an example) provided for the present invention. Figure 8 A comparative schematic diagram of linear and nonlinear optimized wing structures under the same nonlinear aerodynamic load provided by the present invention; Figure 9 This invention provides a complete schematic diagram of an optimized wing structure; Figure 10 This is a schematic diagram showing the structural analysis comparison between the original wing structure and the optimized wing structure under the same aerodynamic load, provided by the present invention. Detailed Implementation

[0016] The following description is intended to disclose the invention and enable those skilled in the art to implement it. The preferred embodiments described below are merely examples, and other obvious variations will occur to those skilled in the art. Example

[0017] This embodiment provides, for example Figure 1 The method shown is a topology optimization method based on aerodynamically driven particle flow for flexible airfoils, which specifically includes the following steps: S1. NURBS-based model construction and initialization: Generate uniformly distributed particles in the parameter space and distribute them on the wing geometric model in the physical space to complete the initialization; Define design parameters, including material properties, structural dimensions, aerodynamic parameters and boundary conditions; Utilize the smooth particle hydrodynamic kernel function to project and map particle information to control points and Gaussian points, and establish the relationship between design variables and geometric control points. S2. Aerodynamic load calculation: The source panel method based on NURBS is adopted. The numerical integration method based on NURBS uses NURBS basis functions to describe the surface velocity potential, calculates the aerodynamic load on the NURBS control points, and applies the corresponding boundary conditions and loads. The calculated aerodynamic load is introduced into the subsequent structural analysis as a constant external load. S3. Based on the nonlinear solution of the uniformly geometric particle flow, perform geometric nonlinear finite element analysis, calculate the tangent stiffness matrix and residual force for iterative solution, and use the Newton-Raphson iterative method to update the displacement field in the design domain to capture the nonlinear response of the structure. S4. Calculation and optimization update of the objective function: Calculate the objective function and perform sensitivity analysis. Use the KD tree to search for the nearest neighbor region and transfer the sensitivity information of the control points back to the particles. Use the moving asymptote method to update the design variables and determine whether the convergence condition is met. If it is met, proceed to the next step; otherwise, return to step S3 to enter the next iteration. S5. Construction of density distribution function: Introduce density distribution function and use relaxation factor to map density information at control points to the entire design domain, transform from discrete control points to continuous density field, describe the material distribution topology of the wing's internal structure, and output the final optimization result.

[0018] Panel analysis, due to its combination of computational efficiency and engineering accuracy, has been widely used in the aerodynamic performance evaluation of aircraft in the preliminary design stage. For the aerodynamic analysis needs of complex geometries, we employ NURBS surfaces to accurately describe the object surface, thus avoiding the approximation errors introduced at the geometric representation level by traditional low-order discretization. In the proposed method, the velocity potential distribution on the object surface is described by a source panel model, and the source intensity is discretized using NURBS basis functions consistent with the geometric representation. This geometric modeling strategy achieves a high-order continuous representation of the source intensity field while ensuring geometric accuracy, resulting in smoother and more stable boundary integral solutions. This section focuses on the NURBS-based source panel method for calculating aerodynamic loads on the surface of objects of arbitrary shapes. This method improves the accuracy and numerical stability of aerodynamic load prediction under complex geometric conditions while maintaining the computational efficiency advantages of panel methods. Consider a moving object in an inviscid incompressible fluid... An object exists in a flow field with a disturbed velocity potential. It can be represented as the source distribution on the flow field boundary:

[0019] in, Represents the boundary surface. It is any field point in the fluid domain. Indicates that it is located at Source point on; function This represents the source intensity distributed along the boundary. It is the distance from the field point to the source point; when the field point Located on the surface of an object At this point, by differentiating along the outward normal direction and applying an impenetrable boundary condition, the boundary integral governing equation is obtained, expressed as:

[0020] in, Points on a curved surface The unit normal vector, This is the normal fluid velocity, given by the boundary conditions of the object's surface. This integral equation belongs to the second kind of Fredholm integral equations and can be solved using numerical methods such as the panel method. When point... With the source When the surfaces coincide, the integral is singular. This special case can be handled appropriately when the object's surface has a closed boundary. For points on the surface where the normals are continuous, the integral can be expressed using Gauss's formula as follows:

[0021] in Let represent the unit normal vector of the source point. The nonsingular expression for the normal velocity can be obtained, written in the following form: ; Source density distribution Satisfy the following homogeneous integral equation:

[0022] The integral on the right-hand side of the homogeneous integral equation is nonsingular, and when the source point... Approaching the field point When, the integrand It tends to zero; therefore, the function It can be determined by homogeneous integral. (The point is missing from the original text.) When the equation also lies on the surface of the object, the singularities that appear in the homogeneous integral equation are correctly eliminated. Therefore, it can be restated as: ; Differentiating the velocity potential with respect to spatial coordinates yields the gradient, which corresponds to the perturbation velocity at any point within the fluid region. ; The source strength can be obtained by using the panel method. Furthermore, its velocity potential is calculated. The boundary geometry of the proposed method is accurately represented by a NURBS surface, and the expression for any point on the object surface can be derived: ; in, These are coordinates in parameter space. These are the coordinates of the A-th control point, and These are NURBS basis functions, defined as follows:

[0023] in and They are order and B-order spline basis functions, It is the weight, where These are B-spline basis functions:

[0024] The definition can be directly derived from The definition derives from this that on the surface of an object;

[0025] in

[0026] The symbols appearing in the above expression have the same definition. Therefore, it can be rewritten as: ; Numerical solutions are obtained by discretizing the boundary surface into a series of panels. The surface integral on each panel is calculated using the Gauss-Legendary quadrature method, where the continuous source density on the boundary is approximated by a finite number of discrete source points. Following Kouh's method, the Gaussian quadrature points serve as both source points and collocation points. Due to the parametric representation of the NURBS surface, the surface tangent vector can be directly obtained, allowing the construction of the surface normal vector and the corresponding area Jacobian matrix. The surface tangent vector is defined as follows:

[0027] in, and They represent the surfaces in and Tangent basis vectors in the direction. These tangent vectors can be used to further construct the geometric quantities of the surface. Unit normal vector. Surfaces in physical space and their corresponding area elements It can be represented as:

[0028] in Represents the cross product of vectors. The Euclidean norm of a vector. The area corresponding to the Jacobian matrix from the parameter space to the physical space is used to transform the integral in the parameter domain into the integral on the physical surface. The above relationship provides the geometric quantities required for subsequent numerical calculations of the boundary integral. Based on the above geometric description, the source intensity distribution on the boundary surface needs further discretization. To maintain consistency between the geometric representation and the physical field approximation, this invention uses the same NURBS basis functions as the geometric modeling to approximate the source intensity function, whose expression is:

[0029] Substitute the formula into the boundary integral control equation and place the points on the boundary. By applying impenetrable boundary conditions, the continuous boundary integral equation can be transformed into a discrete algebraic equation, as follows:

[0030] in ; in It is the Euclidean distance between the placement point and the integration point. This represents the gradient operator with respect to the coordinates of the placement point. Applying the equation to all placement points constructs a linear algebraic system concerning the unknown source intensity: ; in It is composed of integral coefficients The coefficient matrix consists of the source intensity distribution and the free terms. After obtaining the source intensity distribution, the perturbation velocity at any boundary point in the flow field can be directly calculated from the velocity potential gradient, and its expression is:

[0031] in, Point The velocity vector at that point. Further removing its normal component yields the tangential velocity distribution on the boundary surface: ; in This represents the projected component of the velocity along the normal direction. After obtaining the velocity field at the boundary surface, based on the assumptions of incompressible, inviscid potential flow, the aerodynamic pressure distribution can be further calculated using Bernoulli's equations. Boundary point aerodynamic pressure at the location and the corresponding pressure coefficient It can be represented as follows:

[0032] in, and These represent the static pressure and velocity of the free flow, respectively. For fluid density, For point The magnitude of the velocity vector at that point, and the aerodynamic pressure along... Directional action. Perpendicular to the structural surface, forming distributed aerodynamic loads. By integrating the pressure on the boundary surface, it can be converted into an equivalent nodal load vector acting on the structural degrees of freedom, expressed as: ; By substituting the aforementioned aerodynamic load vectors into the structural equations, a weak coupling analysis between aerodynamic forces and structural response can be performed, and the impact of aerodynamic loads on structural deformation can be evaluated. This method effectively alleviates the numerical instability and integral singularity problems commonly found in traditional panel methods by ensuring the continuity of the velocity potential field, thereby significantly improving the accuracy and numerical stability of the calculation results.

[0033] The Isogeometric Material Point Method (I-MPM) is a hybrid Lagrangian-Eulerian numerical method that integrates isogeometric analysis (IGA) into the traditional material point method (MPM). This method uses material points as Lagrangian carriers to describe the evolution of the material state, while introducing a NURBS-based Eulerian background mesh to solve the governing equations. This avoids mesh distortion problems and improves the continuity and numerical stability of the spatial approximation. In I-MPM, the continuum is first discretized into a finite number of material points under a reference configuration. Each material point contains complete physical information, including position, volume, density, stress, and internal variables related to material behavior. Let the position vector of the material points in the initial undeformed configuration be... During loading or iteration, its corresponding spatial location can be determined by a mapping function. This is represented as follows: ; in, This represents the spatial coordinates under the current configuration. This is a pseudo-time or loading step parameter used to describe the gradual application of external loads or boundary conditions. Based on the above mapping relationship, the velocity and acceleration of a particle are defined as its position relative to the loading parameters. First and second derivatives:

[0034] in This indicates the rate of evolution of a matter point in the current configuration. This corresponds to its evolutionary acceleration. To uniformly describe the evolution of different physical quantities during loading, an arbitrary field quantity defined throughout the entire object domain is introduced. In the Lagrange description, the field quantity It is represented by a single material point and varies only with the loading parameters. Its mass derivative is defined as: ; This form emphasizes the evolution of the physical state of a single matter point, embodying the fundamental idea of ​​the matter point method. In the Eulerian description, field quantities are defined at spatial locations, and the effects of time variation and spatial migration must be considered simultaneously. Its matter derivative can be expressed as: ; in Represents the gradient operator with respect to spatial coordinates. This is the evolutionary velocity field under the current configuration, used to describe the migration effects of matter points in space. This form is crucial when mapping particle information to an Eulerian grid. In I-MPM, the governing equations are not solved directly at the particle level, but rather spatially discretized using an Eulerian background grid. To address this issue, a finite-dimensional approximation space is introduced. And use particle integrals to numerically approximate continuous volume integrals: ; in, Indicates the reference configuration. This represents the position of the p-th material point. This indicates its corresponding volume weight. This represents the total number of matter points. This integration method allows the evolution of all physical quantities to be described through the information carried by the particles. However, traditional MPMs typically use low-order finite element basis functions for the mapping between particles and the mesh, whose gradients are discontinuous at element boundaries. When matter points cross element boundaries, element crossing errors are easily generated, leading to abrupt changes in internal forces and instability in sensitivity calculations. To overcome this problem, I-MPM introduces a smooth radial kernel function in the information transfer stage from particles to control points to improve the continuity of the mapping process.

[0035] In the particle-to-control point (P2C) stage, the aggregated variable at the control point The weighted sum of contributions from its neighboring particles is given as follows: ; in, and Let represent the spatial position of the q-th particle relative to the control point. It is the physical quantity carried by the particle. It is a cubic spline smoothing kernel function. This represents the radius of influence of the kernel function. Volume weights correspond to control points. This weighted mapping effectively reduces numerical oscillations caused by cross-domain particle interference. After P2C mapping, the governing equations are solved on a NURBS-based IGA grid. Due to the high-order continuity of NURBS basis functions, control point variables can be smoothly propagated to Gaussian integration points and mapped back to particles after the field equations are solved, thus forming a stable and consistent two-way information transfer mechanism. The approximate expression for any physical field in the physical domain can be uniformly written in the following form: ; in, This represents the degrees of freedom of the A-th control point. These are the corresponding physical domain NURBS basis functions, obtained through geometric mapping from the parametric domain basis functions. This refers to the number of control points in the influence domain. Through the above modeling and discretization process, I-MPM organically combines the advantages of the material point method in large deformations and historical variable tracking with the characteristics of isogeometric analysis in terms of high continuity and geometrically accurate expression, providing a stable, smooth, and physically consistent numerical basis.

[0036] Figure 2 This invention provides a schematic diagram considering both undeformed and deformed configurations with geometric nonlinearity. Within the aforementioned framework of equal geometric material points, the geometric nonlinearity of the structure is modeled using a general Lagrangian description. Let... This represents the material coordinates of the structure in the reference configuration. This represents its corresponding spatial coordinates under the current configuration. The motion of the structure from the reference configuration to the current configuration under external loads can be represented by the Lagrange displacement field. To describe, the two satisfy the following mapping relationship: ; in, This represents the displacement vector defined on the reference configuration, used to characterize the displacement change of a material point from its initial position to its current spatial position. This mapping relationship constitutes the fundamental starting point for describing structural deformation in geometric nonlinear analysis. Based on the above motion description, the deformation gradient tensor can be further defined as follows:

[0037] in It is a second-order unit tensor. This represents the gradient operator with respect to material coordinates. Deformation gradient tensor. This describes the local mapping relationship within the neighborhood of a point in the material. To reasonably measure the strain state under large deformation conditions, the Lagrange-Green strain tensor is introduced, defined as follows: ; This strain metric remains invariant to rigid body rotation and is applicable to finite deformation problems. Under the Lagrangian description, the second type of Piola-Kirchhoff stress... and strain The following relationship must be satisfied: ; in It is the tangential elastic tensor of the material, the magnitude of which is related to the point density variable of the material. Related. To ensure numerical stability, an interpolation method is introduced. ; in, Let be the initial elastic tensor of the solid material. To prevent the minimum density parameter from having zero stiffness, The penalty exponent is used. This interpolation model provides a reasonable tangent modulus for nonlinear iteration while ensuring numerical stability. In the reference configuration... Under these conditions, the virtual work equilibrium condition of the structure can be expressed as:

[0038] in Represents the volume force density vector. It is applied to the boundary. Surface force, This is a virtual displacement field. Based on this, the residual vector corresponding to the nonlinear equilibrium problem can be defined as follows:

[0039] in, and Let represent the internal force vector and the external force vector, respectively. Within the IMPM framework, the internal force vector is obtained by numerical integration over material points, and its expression is:

[0040] in, This indicates the total number of material points. This represents the integration domain corresponding to the p-th material point. Representation and control points Related isogeometric basis functions, This represents the stress tensor at a point in the material. This form fully embodies the basic idea of ​​I-MPM: using particles as integration points and control points as degrees of freedom. To solve the above nonlinear equilibrium equations, the Newton-Raphson method is used to linearize the residuals, thus obtaining the iterative incremental formula: ; in The global tangent stiffness matrix consists of material stiffness terms and geometric stiffness terms.

[0041] In the isogeometric material point method, the construction of the tangent stiffness matrix combines control point interpolation and particle integration, thereby improving the convergence and numerical stability of the nonlinear solution process while ensuring geometric accuracy.

[0042] Based on the aforementioned topology optimization framework for isogeometric particle flow, this study further considers the topology optimization problem of continuum structures subjected to aerodynamic loads, where only the unidirectional influence of aerodynamic loads on the structural response is considered. That is, the aerodynamic pressure field is first calculated based on the given flow conditions and structural geometry, and then treated as an external load acting on the structure, while neglecting the feedback effect of structural deformation on the flow field. This assumption is applicable to structures with stiffness dominance and small deformation, and significantly reduces the computational cost of coupled analysis.

[0043] The geometric model of the design domain is precisely represented by NURBS basis functions, which enables a unified description of geometry, analysis process, and material distribution. These represent the Neumann boundary and the Dirichlet boundary, respectively. Region The internal material layout is determined by the global density distribution function (DDF). The function is described as being constructed from the physical densities of NURBS control points, which are determined by the particle's design variables through a particle control point projection scheme. Under aerodynamic pressure loads, the topology optimization problem aiming to equivalently minimize compliance can be formulated as follows:

[0044] in This represents a set of design variables based on particles, determined by particle positions. and its material consumption variables Composition. The compliance functional is derived from the strain tensor. Implicitly dependent on displacement field And through density-dependent constitutive tensors Depends on design variables. and These represent the internal force vector and the external force vector, respectively. Representation and control points Related NURBS basis functions, It is its gradient relative to the material coordinates. Represents particles The stress it bears. This formula highlights a key feature of the I-MPM: the particles act as integration points, while the control points provide the global degrees of freedom.

[0045] External force virtual work It consists solely of aerodynamic pressure acting on the surface of the structure.

[0046] in, Boundary point pneumatic pressure at the location, The outward unit normal vector is defined in the reference configuration. The pressure field is obtained independently by the flow solver and treated as specified in the structural analysis, reflecting the unidirectional aerodynamic-structural coupling assumption adopted in this invention. The material layout within the design domain is represented by the density distribution function (DDF).

[0047] in, These are the parametric coordinates of NURBS voxels. Is the order of NURBS basis functions, This refers to the physical density of the control points. These control point densities are determined by particle variables through a numerical projection scheme, thus ensuring consistency between the particle distribution and the density representation of the continuous medium. The volume constraints are as follows:

[0048] in This represents the specified maximum permissible material volume fraction. The internal force vectors in the residual equations are calculated using a particle-based numerical integration method within the isogeometric point method of materials (I-MPM) framework, as follows:

[0049] in, Represents particles Stress at the point, It is its reference position. This corresponds to the particle volume. The formula highlights the core feature of I-MPM, which is using particles as numerical integration points while retaining NURBS control points as global degrees of freedom. This is to solve the nonlinear equilibrium equations. Using the Newton-Raphson method, the following incremental form is obtained:

[0050] in The global tangent stiffness matrix consists of material and geometric contributions. Within the PFTO-I-MPM framework, this matrix is ​​assembled using NURBS basis functions and particle-based numerical integration to ensure robustness and quadratic convergence of the nonlinear solution process.

[0051] The proposed isogeometric particle flow topology optimization method requires calculating the first derivatives of the objective function and volume constraint function with respect to design variables to drive the evolution of particle physics variables, and further characterizes the structural topology using the NURBS density distribution function (DDF). Therefore, this section starts from the first derivative of the objective function under geometric nonlinear conditions and systematically derives its sensitivity expression with respect to design variables. In the optimization process, it is assumed that external loads are independent of design variables. Under large deformation conditions, the strain energy of the structure is defined in the reference configuration. The objective function can be written as follows:

[0052] Objective function for global DDF The first derivative can be written as:

[0053] The derivative of the objective function with respect to the displacement field distribution function (DDF) contains explicit terms and implicit terms introduced through the displacement field, the latter being the core challenge in sensitivity analysis of geometrically nonlinear problems. The weak-form reference configuration of the geometrically nonlinear equilibrium equations under the displacement field distribution function can be written as:

[0054] in

[0055] For nonlinear weak forms regarding Taking the derivative, we get:

[0056] Similarly, since the external loads are independent of the design variables, we have: ; Within the PFTO / I-MPM framework, the global DDF is obtained by interpolating the control point density:

[0057] in It is the control point physical density, which is obtained by aggregating particle variables:

[0058] Therefore, the first derivative of the global distribution function (DDF) with respect to the particle design variables is:

[0059] The first derivative of the objective function with respect to the particle design variables can ultimately be written as:

[0060] In numerical implementation, this sensitivity can be effectively evaluated using the particle-Gaussian point-control point mapping of the I-MPM.

[0061] Numerical examples validated the effectiveness of the proposed aerodynamic topology optimization method for wing structures considering geometric nonlinearity. The content is mainly divided into two parts: firstly, verifying the accuracy of geometric nonlinear finite element analysis, and secondly, ensuring the reliability of structural response calculations under large deformations through comparative analysis; thirdly, applying the optimization framework to the design of a complex 14-layer wing rib structure, verifying the method's aerodynamic topology optimization capability in handling large engineering shell structures. The wing geometric model involved was modeled using Matlab, and then post-processed using Rhinoceros 3D software. The focus was on validating the effectiveness of the proposed aerodynamic nonlinear topology optimization method in the conceptual design stage; therefore, all data in the numerical examples were dimensionless. Given that the aerodynamic load calculation used the panel method, its theoretical assumptions limit its applicability in low-speed, potential flow environments, and lift-dominant conditions; therefore, all optimization designs were performed within this applicable range. In all numerical examples, except for specifically defined geometric thickness and angle of attack, the basic physical parameters remained consistent: the Young's modulus of the solid material... Young's modulus of porous materials Poisson's ratio Penalty parameters Free flow velocity air density The convergence criterion for aerodynamic topology optimization is set as follows: when the difference in the 1-norm of the design variables in two consecutive iterations is less than... The optimization process terminates when the number of iterations reaches the maximum value of 150 steps.

[0062] Figure 3 The aim was to verify the ability of the developed PFTO method to capture the large rotational and large deformation behavior of structures. A cantilever beam model with an aspect ratio of 5:1 and physical dimensions of 50 cm × 10 cm was constructed. Material properties were set using dimensionless parameters or standard units, with the elastic modulus being [value missing]. Poisson's ratio is The boundary conditions are strictly set according to the classical cantilever beam theory: a full displacement constraint is applied to the left boundary, and a concentrated downward force is applied to the right end face. Figure 3 (b) shows the deformed mesh in the final equilibrium state simulated by this method. The results show that the mesh maintains good quality under large deformation and no distortion occurs. Figure 3 (c) Further details the mechanical response characteristics of the structure: the gray dashed line in the figure represents the force-displacement path predicted by the linear stiffness, exhibiting a linear distribution; the blue solid line represents the actual response after considering geometric nonlinear effects. It can be seen that in the small deformation stage, the two essentially coincide. However, with increasing load, due to the large rotation of the structure, the geometric stiffness matrix changes, and the nonlinear result (blue) gradually deviates from the linear prediction. The accurate reproduction of the nonlinear hardening / softening behavior strongly confirms the accuracy of the PFTO method in solving nonlinear continuum mechanics problems.

[0063] Reference Figure 4 As shown, to further verify the applicability of the PFTO method in the design of complex engineering structures, Figure 4 (a) The mesh convergence test process is illustrated in detail. As the mesh density gradually increases from 40×250 to 120×750, the objective function value initially fluctuates before stabilizing. Data shows that when the mesh size reaches 80×500, the objective function value stabilizes around 40.069, and the corresponding topological features are clearly visible. Considering both computational accuracy and cost, it can be determined that an 80×500 mesh has met the convergence requirements. Therefore, to ensure consistency and computational efficiency, this mesh size will be used in all subsequent examples. Furthermore, Figure 4 (b) The initial geometry and optimized deformable configuration of the wing model are visually demonstrated, further illustrating the stability of the method in handling geometric nonlinear problems involving large structural deformations.

[0064] Reference Figure 5 As shown, a geometric model with 14 ribs was established based on classic NACA airfoil data, and nonlinear topology optimization design was performed for verification. The naming convention clearly defines the geometric characteristics of the airfoil: "2" indicates that the maximum camber is 2% of the chord length, "4" indicates that the maximum camber is located at 40% of the chord length, and "15" indicates that the maximum relative thickness is 15%. This study selected this airfoil as the design object for the 14th rib. The specific design domain definition and boundary conditions are as follows... Figure 5 As shown in (b), in the model, the sparse support locations are simplified to fixed constraints. To maintain the external aerodynamic shape of the wing and ensure assembly accuracy, both the fixed constraint region and the external profile of the airfoil are set as undesignable regions.

[0065] Reference Figure 4 As shown, to illustrate the convergence properties, Figure 6 Details of the two design phases are presented. Analysis shows that the particle behavior exhibits clear Lagrangian characteristics. Starting from an initial uniform distribution (equal material consumption), during the iteration process, the particles "flow" and redistribute within the design domain, driven by sensitivity analysis, until they converge to the target position. Notably, at the end of optimization, some particles moved outside the IGA mesh boundary. This is due to the "reverse correction sensitivity" mechanism preventing excessive particle aggregation. Since their displacements remain within the support domain of the SPH kernel function, these particles can still effectively participate in boundary optimization without negatively impacting the final topology. Therefore, by comprehensively considering the physical location and material consumption of the particles and mapping them to control points to construct the density distribution function, the PFTO method successfully achieves topology optimization for the ribs. This result fully demonstrates the superior effectiveness, excellence, and irreplaceable engineering value of this method in handling curved surface design domains.

[0066] Reference Figure 7 The second-layer IGA mesh reveals the optimal material distribution for the ribs, forming a continuous material distribution structure that matches the aerodynamic load distribution. The optimized structure exhibits higher density in the leading and central regions, corresponding to the lift regions where negative pressure loads are concentrated. Material is effectively removed from the trailing edge and regions with lower loads, achieving a lightweight design. This validates the effectiveness of the PFTO method in the surface design domain. Thanks to the powerful geometric modeling capabilities of NURBS, the optimized structure has smooth and continuous boundaries, with clear interfaces between the real and void regions. The high consistency between the geometric model and the numerical analysis model significantly improves numerical accuracy and optimization stability.

[0067] Depend on Figure 7 and Figure 8 It is evident that introducing geometric nonlinear modeling not only corrects the kinematic errors that linear theory cannot capture in large deformations and avoids unreasonable stress concentration, but also significantly reduces the peak stress and deformation of the structure through a material layout that is more in line with the physical reality of fluid-structure interaction, resulting in a reliable design with higher structural integrity, safety margin and load-bearing efficiency.

[0068] Depend on Figure 9 The corresponding optimized embodiments illustrate that the optimized structure not only maintains excellent geometric continuity, but also can adaptively match complex aerodynamic loads distributed along the spanwise direction, fully demonstrating the versatility and application potential of the NURBS-based nonlinear aerodynamic particle flow topology optimization method in the design of complex aircraft components.

[0069] Depend on Figure 10Numerical comparisons show that the stiffness of the optimized structure has increased by more than an order of magnitude, while the maximum displacement and peak stress have both decreased sharply. This fully verifies the significant advantages of the high-order surface element aerodynamic topology optimization strategy of this invention in generating high-load-bearing, manufacturable rib structures and improving their aerodynamic environmental adaptability.

[0070] This invention effectively enhances the structural response and engineering applicability of flexible airfoils under extreme nonlinear aerodynamic conditions by integrating Lagrange particle description and Eulerian mesh description. This method combines the numerical stability of isogeometric analysis and the material point method with the geometric tunability of particle flow, achieving continuous evolution of material distribution and seamless integration of geometric analysis within a unified framework. Specifically, on the one hand, the information coupling mechanism between Lagrange particles and the Eulerian background mesh (combined with SPH kernel functions and NURBS basis functions) employed by I-MPM successfully overcomes the limitations of traditional finite element methods in handling mesh distortion, significantly improving the physical stability and solution accuracy of flexible structures during large aerodynamic deformation processes. On the other hand, the intelligent "flow" of particles breaks through the limitations of traditional fixed control points, enabling the structural topology to dynamically respond to aerodynamic load paths and high strain energy distributions during nonlinear fluid-structure interaction processes, thereby generating more rational, clearly defined, and numerically independent flexible skeleton structures. By optimizing the aggregation and sparsity behavior of particles under aerodynamic drive during the process, the material is adaptively reinforced in the critical force transmission area and eliminated in the non-critical area, effectively meeting the dual requirements of flexible wings to resist aerodynamic bending and torsional deformation and maintain high structural integrity.

[0071] Through in-depth comparative analysis of flexible airfoils and related benchmark structures, the versatility and robustness of the nonlinear PFTO method were verified. Results show that this method not only effectively avoids the mesh distortion and evaluation failure problems that traditional linear analysis easily encounters when dealing with large deformations, but also generates high-quality aerodynamic structural configurations with excellent buckling resistance, extremely low maximum displacement, and equivalent stress (von Mises stress). Based on its high-fidelity structural representation and robust nonlinear optimization process, this computational strategy is expected to be further extended to the design of more advanced and complex aerospace equipment, including dynamic aeroelastic constraints, complex three-dimensional flow field multiphysics coupling, and deformable aircraft and metamaterial aerodynamic configurations.

[0072] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the claimed invention. The scope of protection claimed by the appended claims and their equivalents is defined.

Claims

1. A topology optimization method based on aerodynamically driven particle flow for flexible airfoils, characterized in that, include: S1. NURBS-based model construction and initialization: Generate uniformly distributed particles in the parameter space and distribute them on the wing geometric model in the physical space to complete the initialization; Define design parameters, including material properties, structural dimensions, aerodynamic parameters and boundary conditions; Utilize the smooth particle hydrodynamic kernel function to project and map particle information to control points and Gaussian points, and establish the relationship between design variables and geometric control points. S2. Aerodynamic load calculation: The source panel method based on NURBS is adopted. The numerical integration method based on NURBS uses NURBS basis functions to describe the surface velocity potential, calculates the aerodynamic load on the NURBS control points, and applies the corresponding boundary conditions and loads. The calculated aerodynamic load is introduced into the subsequent structural analysis as a constant external load. S3. Based on the nonlinear solution of the uniformly geometric particle flow, perform geometric nonlinear finite element analysis, calculate the tangent stiffness matrix and residual force for iterative solution, and use the Newton-Raphson iterative method to update the displacement field in the design domain to capture the nonlinear response of the structure. S4. Calculation and optimization update of the objective function: Calculate the objective function and perform sensitivity analysis. Use KD tree to search for the neighborhood and transfer the sensitivity information of the control point back to the particle. The design variables are updated using the moving asymptote method, and it is determined whether the convergence condition is met. If it is met, proceed to the next step; otherwise, return to step S3 and proceed to the next iteration. S5. Construction of density distribution function: Introduce density distribution function and use relaxation factor to map density information at control points to the entire design domain, transform from discrete control points to continuous density field, describe the material distribution topology of the wing's internal structure, and output the final optimization result.

2. The topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The specific operations for NURBS-based model construction and initialization in step S1 are as follows: Uniformly distributed particles are generated in the parameter space and distributed onto the wing geometry model in physical space through mapping; design parameters including material properties, structural dimensions, aerodynamic parameters and boundary conditions are defined. A smooth particle hydrodynamic kernel function is used to determine the neighborhood with the grid point as the center. The relative distance between the particle position in the neighborhood and the center point is calculated and the kernel function is called to calculate the weight. The material information of the particle is projected and mapped to the control point and Gaussian point respectively.

3. The topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The specific steps for calculating aerodynamic loads in step S2 are as follows: The aerodynamic field is solved using the NURBS-based source panel method; In the calculation process, the numerical integration method based on NURBS is used to describe the velocity potential of the wing surface using NURBS basis functions; Solve the surface element equations to calculate the aerodynamic load distribution at the NURBS control points and apply the corresponding aerodynamic boundary conditions; the extracted aerodynamic loads will be directly introduced as constant external excitation loads into the geometric nonlinear structural analysis in step S3.

4. The topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The specific process of the nonlinear solution based on the isogeometric particle flow in step S3 is as follows: For the geometric nonlinear large deformation characteristics of the wing, structural analysis is performed based on the isogeometric finite element analysis framework; the current tangential stiffness matrix of the structure and the residual forces between internal forces and external aerodynamic loads are calculated. Based on the Newton-Raphson iterative method, the displacement field in the design domain is continuously solved and updated by residual force until the equations meet the convergence accuracy, thus capturing the nonlinear mechanical response of flexible airfoils under large deformation.

5. A topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The specific process of calculating the objective function and transferring sensitivity in step S4 is as follows: The objective function for topology optimization based on nonlinear displacement field calculation and its sensitivity to Gaussian point density are introduced. A sensitivity analysis transfer mechanism is introduced, and the KD tree algorithm is used to search the neighboring particle region around the control point. Based on the chain rule, the sensitivity information of the control point is transferred back to each particle to obtain the sensitivity of the objective function to the design variables of the bottom particle.

6. The topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The specific process of optimization update and convergence judgment in step S4 is as follows: Based on the acquired particle sensitivity, the moving asymptote method is used to iteratively update the particle position and attribute design variables; the rate of change of the objective function value is monitored and the volume constraint is detected, and the system determines whether the set convergence condition is met: if the convergence condition is met, the optimization is determined to end; Otherwise, based on the updated particle flow state, return to step S3 for the next nonlinear calculation and optimization iteration.

7. A topology optimization method for flexible airfoils based on aerodynamically driven particle flow according to claim 1, characterized in that, The construction process of the density distribution function in step S5 is as follows: Density information at the converged control points is extracted, and a density distribution function is introduced. Using a set relaxation factor, the density information at the discrete control points is smoothly mapped to the entire airfoil design domain. By utilizing the high-order analytical properties of NURBS basis functions, a seamless transformation from the discrete control point set to the continuous physical density field is achieved, describing and extracting the high-fidelity topological boundary of the material distribution of the internal load-bearing skeleton of the wing.

Images