High aspect ratio wing topological optimization method considering pneumatic-structural coupling

Abstract

The invention provides a high aspect ratio wing topological optimization method considering pneumatic-structure coupling. The method comprises the following steps: establishing a pneumatic-structure strong coupling topological optimization model; performing finite element analysis to obtain an original variable field of the high aspect ratio wing; generating a filtering variable field and a projection variable field; deformation displacement of each finite element unit is obtained based on a UMMI material interpolation model; establishing information transmission, structure displacement response, aerodynamic load transmission and aerodynamic load updating strategies under pneumatic-structure coupling, and analyzing and updating the aerodynamic load; according to the method, manufacturing constraints such as geometric nonlinearity, multi-material layout and wing rib orientation size are comprehensively considered, design variables in an original variable field are subjected to iterative optimization solution, a wing optimization configuration is obtained, the three core problems of calculation distortion, performance deviation and manufacturing difficulty in lightweight design of the high aspect ratio wing are successfully solved, and the method is suitable for the lightweight design of the high aspect ratio wing. And the comprehensive performance and economical efficiency of the aircraft are remarkably improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of structural optimization, and particularly relates to a large-aspect-ratio wing topology optimization method considering aerodynamic-structural coupling. BACKGROUND

[0002] In modern engineering design, especially in the field of aerospace, pursuing lightweight and high performance of structures has become a general trend. Topology optimization technology, as an advanced and effective design method, can automatically explore the optimal material layout scheme through mathematical models and algorithms to achieve the purpose of reducing weight while meeting specific performance requirements. In recent years, with the development of material science, the optimization design of multi-material structures has attracted widespread attention because it can further expand the design space, making the structure not only lightweight but also more adaptable to complex working environments.

[0003] Currently, multi-material topology optimization mainly relies on three mainstream methods: density method, structural progressive method and level set method. Among them, the density method is more widely used in solving practical engineering problems due to its flexibility and ease of implementation. However, for large-aspect-ratio wings, the traditional constant load topology optimization method faces severe challenges. Large-aspect-ratio wings will experience significant deformation during flight, which directly affects the distribution of unmanned aerial vehicle aerodynamic loads and changes the stress state of the structure. Therefore, the traditional constant load optimization strategy cannot accurately reflect this dynamic process, and the reliability of its results is questionable.

[0004] In recent years, researchers have begun to realize the importance of aerodynamic-structural coupling effects on optimization design. Aerodynamic-structural coupling refers to the interaction between aerodynamic forces and structural responses (such as deformation and vibration) during the flight of an aircraft. Currently, most research on wing aerodynamic-structural coupling focuses on local issues such as wing aerodynamic profile optimization and flutter analysis, without systematically integrating aerodynamic-structural coupling effects into the overall topology optimization design of large-aspect-ratio wings. SUMMARY

[0005] The present application provides a large-aspect-ratio wing topology optimization method considering aerodynamic-structural coupling to solve the problem of multi-material topology optimization of large-aspect-ratio wings under the condition of real-time changes in aerodynamic loads.

[0006] To achieve the above-mentioned purpose, the technical solution provided by the present application is as follows:

[0007] A large-aspect-ratio wing topology optimization method considering aerodynamic-structural coupling, comprising the following steps:

[0008] Step 1: Determine the optimization objective and constraints for the topology optimization problem of a high aspect ratio wing multi-material structure. Based on the optimization objective and constraints, establish a high aspect ratio wing multi-material topology optimization model considering aerodynamic-structural coupling. The optimization objective is to minimize the wing's terminal compliance. The optimization constraints include rib orientation dimension control constraints, local bending deformation control, overall wing mass constraints, and single material volume fraction constraints.

[0009] The multi-material topology optimization model for a high aspect ratio airfoil considering aerodynamic-structural coupling is as follows:

[0010]

[0011] In the formula: Represents the set of design variables. For design variables, representing the first... The finite element and the first To avoid singularities in the structural stiffness matrix caused by design variables at 0 during finite element analysis, a lower limit for design variables is introduced. ;when When, it indicates that the unit is entirely composed of the j-th phase material, when When the value is between the lower limit of the design variable and 1, it indicates that the unit is composed of the first... The degree of phase material composition; in order to better determine the first The UMMI multi-material interpolation model uses an exponential function penalty to determine which phase of material each finite element is ultimately composed of, ensuring that the final variable field converges to the lower limit of the design variables or 1; the initial values ​​of all design variables are set to 0.1. For the first Projected variable field of phase material, For the first Phase solid material filtering variable field, For the first The original variable field of the phase material, Indicates the first Phase solid materials in the first Projection variables of the unit; Represents the total number of design domain units; To optimize the objective, This refers to the structural displacement vector in the Nth step of the iterative solution. For the internal forces at the unit nodes, To employ an update strategy based on wingtip deformation to participate in the updating of aerodynamic loads, The residual force vector; For the Heaviside projection function, This is a density filtering matrix; For the overall structural quality. is the first density of the phase material, is the preset upper limit of the overall structure mass, is the first volume of the finite element; two materials are used in the model to simulate the main beam and the rib in the wing, wherein represents the hole size control constraint of the rib material; represents the volume fraction constraint of the material j, is the set upper limit of the volume fraction, is the overall structure volume, is the bending deformation control constraint;

[0012] Step 2: according to the large-aspect-ratio-wing multi-material topology optimization model, finite element mesh is divided for the large-aspect-ratio-wing to generate a plurality of finite element units; according to the boundary constraint and the initial load condition of the large-aspect-ratio-wing structure, finite element analysis is performed on the large-aspect-ratio-wing multi-material topology optimization model to obtain the original variable field of the large-aspect-ratio-wing;

[0013] wherein the large-aspect-ratio-wing is composed of two-phase materials; different variables in the original variable field correspond to different finite element units of the wing, and the variable is used to describe whether any phase material of the two-phase materials exists in the corresponding finite element unit;

[0014] Step 3: the original variable field is sequentially subjected to density filtering and Heaviside projection operation to generate a filtered variable field and a projected variable field; wherein:

[0015] the filtered variable field of each phase material is obtained from the original variable field of the phase material and a set filtering matrix; the projected variable field is obtained by projecting the filtered variable field through a projection function, and the projection variable in the projected variable field of each phase material corresponds to the filtered variable in the filtered variable field of the phase material one by one;

[0016] Step 4: a UMMI material interpolation model is established according to the projected variable field obtained in step 3 to characterize the mapping relationship between the projected variable field and the unit material property; based on the UMMI material interpolation model, a finite element analysis considering nonlinear geometry is performed by calling a finite element commercial solver to obtain the deformation displacement of each finite element unit in the large-aspect-ratio-wing;

[0017] Step 5: according to the deformation displacement of each finite element unit, a data information transmission strategy under aerodynamic-structure coupling, a structure displacement transmission strategy, an aerodynamic load transmission strategy and an aerodynamic load updating strategy are established; and aerodynamic-structure analysis is performed on the wing, and whether the aerodynamic load of each unit needs to be updated is judged according to the analysis result, if yes, the aerodynamic load is updated and the step 2 is returned; otherwise, step 6 is executed;

[0018] Step 6: Establishing material hole directional size control constraints by using local cylindrical search and introducing local material rate method; and establishing local bending deformation control constraints of the high-aspect-ratio wing by using unit node displacement;

[0019] Step 7: Determining constraint functions according to the material hole directional size control constraints and the local bending deformation control constraints of the high-aspect-ratio wing;

[0020] According to the deformation displacement of each finite element unit, the sensitivity of the objective function and the constraint function to the design variable is calculated respectively; and according to the sensitivity calculation result, the design variable in the original variable field is iteratively optimized and solved by using a gradient-based optimization algorithm under multiple constraint conditions, until the objective function reaches the optimal under the constraint conditions, and then the optimized configuration of the high-aspect-ratio wing is obtained;

[0021] The multiple design constraint conditions include mass constraint, volume constraint, material hole directional size control constraint, wing bending deformation control constraint and design variable constraint.

[0022] Step 8: Lightweight reconstruction design is performed on the optimized configuration of the high-aspect-ratio wing by using a truss structure.

[0023] Further, the step 1 includes the following sub-steps:

[0024] Step 1.1: Establishing a wing end point compliance expression, and taking the minimum wing end point compliance as an optimization objective; wherein the wing end point compliance expression is:

[0025]

[0026] In the formula, is an external force of a unit node force, is a final displacement obtained by nonlinear iterative solution;

[0027] Step 1.2: Establishing a wing finite element analysis model considering nonlinear geometry, performing finite element meshing on the wing to obtain a plurality of finite element units; wherein a shell unit is used to simulate the wing skin, and a solid unit is used to simulate the design domain, and each finite element unit is assigned with two candidate materials, respectively representing the wing main beam and the wing rib.

[0028] Step 1.3: Establishing constraint conditions, including: wing overall mass constraint, single material volume fraction constraint, directional size control constraint and bending deformation control constraint;

[0029] Step 1.4: Constructing a high-aspect-ratio wing multi-material topology optimization model considering aerodynamic-structure coupling according to the optimization objective, the wing original design variable field and the constraint conditions.

[0030] Furthermore, step 3 includes the following sub-steps:

[0031] Step 3.1: Pre-set the filter radius, calculate the filter variable for any j-th phase material in the original variable field, and output the filter variable field. :

[0032]

[0033] In the formula, , These represent the number of materials involved in the design and the number of designable units, respectively. This includes the main beam material and the wing rib material; For the original variable field, ; Indicates the first Unit and the Design variables related to phase materials; For the filtering matrix, ; Indicates the k-th unit about the th The weight coefficient of the unit is expressed as follows:

[0034]

[0035] In the formula, For the first The volume of the unit For the first The volume of the unit Indicates the k-th unit about the th Distance weight of the cell Indicate the ξth unit about the ξth unit. Distance weight of the cell;

[0036] definition: ,in Given a pre-defined filter radius, For the first Unit and the

[0037] Center distance of units;

[0038] Step 3.2: Use the Heaviside function to filter the variable field. Perform a projection operation to generate a projected variable field. Among them, the projected variable field The projection variable of the i-th element with respect to the j-th phase solid material for:

[0039]

[0040] In the formula, These are projection parameters used to control the steepness of the projection. During the optimization process, an incremental adjustment strategy is used to dynamically increase the value. The projection threshold is 0.5.

[0041] Furthermore, in step S3.2, an incremental adjustment strategy is adopted to dynamically adjust... The process is as follows:

[0042] In iterations 1–30, the projection parameters Set to a fixed value of 1;

[0043] Starting from step 31, the projection steepness control parameter is adjusted every 10 iterations. Multiply by 1.5 and repeat this process until the maximum number of iterations is reached.

[0044] Furthermore, in the UMMI multi-material interpolation model established in step 4, the UMMI interpolation of Young's modulus is expressed as the sum of the material weights of each phase:

[0045]

[0046] In the formula, For projection variables The weight function, For the first Young's modulus of the phase material, where m is the number of materials involved in the optimization. ;

[0047] Using the SIMP penalty function to adjust the weight function Defined as:

[0048]

[0049] in, As a penalty factor; For the first Unit about the first Projected variables of phase solid materials, for Unit about the first Projected variables of phase solid materials.

[0050] Furthermore, in step 4, the specific process of performing finite element analysis based on the UMMI material interpolation model is as follows:

[0051] Using the UMMI material interpolation model, the material properties of each finite element are calculated based on the variables of each finite element.

[0052] The material properties obtained are applied to the corresponding finite element units to form a finite element analysis model reflecting the optimized design;

[0053] The finite element analysis model is imported into a finite element commercial solver to perform finite element iterative optimization solving, and a finite element analysis result file is outputted, from which the deformation displacement of each finite element unit is extracted.

[0054] Further, in step 5, the data information transmission strategy under the aerodynamic-structure coupling is:

[0055] Finite element solving: on the basis of the original aerodynamic load, a commercial solver computer wing structure deformation is called to calculate, and the displacement of each node of the wing skin in each direction is outputted;

[0056] Deformation point cloud calculation: the skin node coordinates before the wing structure deformation are added to the corresponding displacement to obtain the appearance point cloud after the skin deformation;

[0057] Reverse modeling: based on the point cloud spatial coordinate information after the skin deformation, the convergence body interpolation, creative modeling, inclusive body and split body functions of the UG commercial software are utilized to realize the reverse modeling of the deformed wing;

[0058] Deformed aerodynamic load calculation: the deformed wing geometric model is imported into Fluent to perform CFD simulation, the aerodynamic load of the deformed wing is calculated, and a load file is exported;

[0059] Aerodynamic load transmission: the deformed aerodynamic load is mapped to the optimized finite element model to realize the update of the aerodynamic load;

[0060] The structure displacement transmission strategy is:

[0061] The initial coordinates of the wing skin and the deformation in each direction are extracted, and the two are added to obtain the point cloud distribution after the skin deformation; the point cloud is fitted into a convergence body composed of numerous small planes by using the UG commercial software, and the plane convergence body is reversely generated into a geometric framework in the UG creative modeling application module; the obtained geometric framework is sequentially subjected to inclusive body generation and split body processing to obtain the geometric model after the wing deformation; finally, the geometric model is imported into the Fluent module in Workbench to calculate the aerodynamic load of the deformed wing skin appearance;

[0062] The aerodynamic load transmission strategy is:

[0063] Aerodynamic load spatial search mapping: A spatial search method is used to traverse the output aerodynamic load file one by one to find the aerodynamic load point closest to each point in the point cloud model, and the pressure value at the aerodynamic load point is transferred to the corresponding point in the point cloud model; in addition, each interpolated load is multiplied by a lift coefficient, which is the ratio of the total lift in the aerodynamic load file to the total lift after interpolation;

[0064] Transfer of point cloud loads to the wing finite element model: Subtract the deformation of the corresponding node on the skin from the coordinate values ​​in the deformed point cloud load file to obtain the aerodynamic load file of the wing finite element model in optimization.

[0065] Update and optimize the loads on the wing finite element model:

[0066] The aerodynamic load update strategy is as follows: through the wingtip section center node Change in displacement in direction The formula for determining whether to update the aerodynamic load is:

[0067]

[0068] In the formula, This represents the displacement in the y-direction of the center node of the wingtip section in the current iteration step; Center node of wingtip section The reference object for directional displacement is updated every time the aerodynamic load is updated during the optimization process. Then update to the current iteration step. Values, where, before the first update of the aerodynamic loads, This represents the displacement value in the y-direction of the center node of the wingtip section at step 0; The upper limit of the pre-introduced displacement change is set to 30%.

[0069] Furthermore, in step 6, the process of establishing the material pore orientation size control constraint is as follows:

[0070] Based on the characteristics of the traditional wing spars-rib structure, and considering the rib material, in the first... Introducing material rate oriented towards rib material within the local search region of the element :

[0071] ;

[0072] In the formula, Indicates the first The volume of the unit; As a unit The local search area centered on the search center. For unit Projected variables regarding the rib material; This is the penalty factor, and its value is 3.

[0073] Introducing a lower limit for material yield The maximum pore size control constraint expression for the rib material is obtained as follows:

[0074]

[0075] In the formula, In the first A material rate oriented towards the rib material is introduced within the local search region of the element;

[0076] Introducing the p-means condensation function, all local constraints are condensed into a single global constraint. Based on the definition of the condensation function, a constructor is constructed. Its expression is defined as:

[0077]

[0078] Then, establish the global expression for the directional size control constraint:

[0079]

[0080] in, The expression is

[0081]

[0082] In the formula, It is the number of units. It is the penalty factor for the condensation function, with a value of 80;

[0083] The process of establishing the bending deformation control constraint is as follows:

[0084] Multiple nodes evenly positioned on the wingtip section are selected as observation nodes. The displacements of each observation node are extracted and a system is established.

[0085] Maximum nodal displacement control constraint model in the direction:

[0086]

[0087] In the formula, Indicates the node number is The observation nodes are at Displacement in direction, This is the preset upper limit of the constraint;

[0088] Each observation point is Displacement in direction Expressed as:

[0089]

[0090] wherein, is the global displacement vector of the wing structure, is the displacement extraction vector;

[0091] Definition:

[0092] According to the definition of the p-mean agglomeration function, the variable is constructed, and its expression is:

[0093]

[0094] The maximum nodal displacement control constraint model in the direction of is rewritten as:

[0095]

[0096] wherein, is the number of observation nodes, is the index of the agglomeration function;

[0097] The bending deformation control constraint is obtained.

[0098] Further, in the step 7,

[0099] The mass constraint: the mass of the high-aspect-ratio wing calculated based on the projection variable field is less than or equal to the pre-set upper limit of the mass;

[0100] The volume constraint: the volume of the same phase material calculated based on the projection variable of each finite element unit corresponding to the same phase material is less than or equal to the pre-set upper limit of the volume;

[0101] The design variable constraint: the variables in the projection variable field are limited in the open interval (0, 1].

[0102] The advantages of the present application are:

[0103] The present application successfully solves the three core problems of calculation distortion, performance deviation and manufacturing difficulty in the lightweight design of a high-aspect-ratio wing by establishing an aerodynamic-structure strong coupling topology optimization model and comprehensively considering geometric nonlinearity, multi-material layout and wing rib directional size manufacturing constraints. The final effect is to automatically design an ultra-lightweight wing configuration that meets the requirements of aerodynamic performance, structural strength and manufacturability, significantly improving the comprehensive performance and economy of the aircraft. Specifically embodied in:

[0104] 1、Traditional topology optimization is often carried out under fixed load, ignoring the feedback influence of wing deformation on aerodynamic load, leading to design deviation. The method of the present application realizes real aerodynamic-structure strong coupling analysis by establishing aerodynamic load updating strategy (step 5) and bidirectional data transmission strategy, and updating the aerodynamic load changed due to structural deformation in real time in optimization iteration. This makes the finally optimized wing configuration achieve global optimization of stiffness and aerodynamic efficiency under real flight working conditions, significantly reduces the performance loss caused by ignoring coupling effect, and improves the comprehensive performance of the wing.

[0105] 2、Since the high-aspect-ratio wing will produce large deformation under load, it shows strong geometric nonlinearity. The present application actively considers the nonlinear geometric effect in finite element analysis (step 4), rather than using the traditional linear assumption, which greatly improves the calculation accuracy of structural response (stress, strain, displacement), thereby ensuring the correctness of the sensitivity analysis and optimization direction based on the response, and fundamentally avoiding the risk of distorted optimization results and even structural failure caused by linear assumption.

[0106] 3、The traditional topology optimization result often produces complex and porous "truss" structure, which is difficult to meet the requirements of aviation structure on fatigue life and processability. The present application creatively introduces the wing rib directional size control and local cylindrical search method (step 6), which can actively control the minimum size and directionality of the holes in the optimization result, directly generating a topology configuration with clear structure, reasonable layout, easy engineering interpretation and manufacturing, greatly reducing the difficulty and cost of subsequent detailed design and production. BRIEF DESCRIPTION OF DRAWINGS

[0107] The above and / or additional aspects and advantages of the present application will become apparent and more readily appreciated from the following description of the embodiments, taken in conjunction with the accompanying drawings, in which:

[0108] Figure 1 is a flowchart of a high-aspect-ratio wing topology optimization method considering aerodynamic-structure coupling of the present application;

[0109] Figure 2 is a schematic diagram of the wing structure to be optimized provided by the embodiment of the present application;

[0110] Figure 3 is a load schematic diagram of the wing to be optimized provided by the embodiment of the present application;

[0111] Figure 4 is a schematic diagram of the structure displacement and aerodynamic load updating process in the embodiment of the present application;

[0112] Figure 5 is a schematic diagram of the reverse modeling process in the embodiment of the present application;

[0113] Figure 6Fig. 1 is a schematic diagram of the process of delivering aerodynamic load based on spatial search in the embodiment of the present application;

[0114] Figure 7 Fig. 2 is a schematic diagram of directional size control using cylindrical search in the embodiment of the present application;

[0115] Figure 8 Fig. 3 is a schematic diagram of wing bending deformation control in the embodiment of the present application;

[0116] Figure 9 Fig. 4 is a schematic diagram of the comparison between the topological optimization effect of the high-aspect-ratio wing considering aerodynamic-structure coupling and the topological optimization effect of the high-aspect-ratio wing not considering aerodynamic-structure coupling, wherein 9a is a topological optimization effect diagram not considering aerodynamic-structure coupling, and 9b is a topological optimization effect diagram considering aerodynamic-structure coupling;

[0117] Figure 10 Fig. 5 is a diagram of the overall reconstruction effect of the topological optimization result of the high-aspect-ratio wing considering aerodynamic-structure coupling by using a truss structure in the embodiment of the present application;

[0118] Figure 11 Fig. 6 is a diagram of the top view of the high-aspect-ratio wing considering aerodynamic-structure coupling in the embodiment of the present application; Figure 10 Fig. 7 is an enlarged view of position A in Fig. 6;

[0119] Figure 12 Fig. 8 is a diagram of the side view of the high-aspect-ratio wing considering aerodynamic-structure coupling in the embodiment of the present application; Figure 10 Fig. 9 is a diagram of the side view of the high-aspect-ratio wing considering aerodynamic-structure coupling in the embodiment of the present application.

[0120] In the figure: 1-wing root, 2-wing tip, 3-front main beam, 4-rear main beam, 5-secondary beam, 6-wing rib, 7-main beam diagonal cable, 8-main beam reinforcing rod, 9-wing rib diagonal cable, 10-wing rib reinforcing rod, 11-main beam secondary beam. DETAILED DESCRIPTION

[0121] The embodiments of the present application are described in detail below, which are exemplary and intended to explain the present application, and cannot be understood as a limitation of the present application.

[0122] With reference to Figure 1 , the embodiment of the present application provides a topological optimization method of a high-aspect-ratio wing considering aerodynamic-structure coupling, comprising the following steps:

[0123] Step 1: determining the optimization target, wing design variable field and constraint condition of the topological optimization problem of the high-aspect-ratio wing multi-material structure, and establishing a high-aspect-ratio wing multi-material topological optimization model considering aerodynamic-structure coupling according to the optimization target and optimization constraint. Including the following steps:

[0124] Step 1.1: taking the wing end point compliance optimization target to establish the expression of the optimization target. Wherein, the wing end point compliance is defined as:

[0125] (1)

[0126] in For element nodal forces, This is the displacement in the last step of the nonlinear iterative solution process, which is the final displacement obtained by solving the equilibrium equations using the finite element method.

[0127] Step 1.2: Establish a finite element analysis model of the wing considering nonlinear geometry, perform finite element meshing on the wing to obtain several finite element elements; among them, the wing skin is simulated using shell elements, and the design domain is simulated using solid elements. Each finite element element is assigned two candidate materials, representing the wing spars and wing ribs respectively.

[0128] Step 1.3: Establish constraints, including: overall wing mass constraints, single material volume ratio constraints, directional size control constraints, and bending deformation control constraints.

[0129] Step 1.4: Based on the optimization objective, wing design variable field, and optimization constraints, construct a multi-material topology optimization model for a high aspect ratio wing considering aerodynamic-structural coupling. The optimization objective is to minimize the wing's terminal compliance. Optimization constraints include overall structural mass, single-material volume fraction, directional dimension control, and wing bending deformation control. The model expression is as follows:

[0130] (2)

[0131] In the formula, Represents the set of design variables. For design variables, representing the first... The finite element element and the first To avoid singularities in the structural stiffness matrix caused by design variables at 0 during finite element analysis, a lower limit for design variables is introduced. ;when When, it indicates that the unit is entirely composed of the j-th phase material, when When the value is between the lower limit of the design variable and 1, it indicates that the unit is composed of the first... The degree of phase material composition; in order to better determine the first The UMMI multi-material interpolation model uses an exponential function penalty to determine which phase of material each finite element is ultimately composed of, ensuring that the final variable field converges to the lower limit of the design variables or 1; the initial values ​​of all design variables are set to 0.1. For the first Projected variable field of phase material, For the first Phase solid material filtering variable field, For the first The original variable field of the phase material, Indicates the first Phase solid materials in the first Projection variables of the unit; Represents the total number of design domain units; To optimize the objective, This refers to the structural displacement vector in the Nth step of the iterative solution. For the internal forces at the unit nodes, To employ an update strategy based on wingtip deformation to participate in the updating of aerodynamic loads, The residual force vector; For the Heaviside projection function, This is a density filtering matrix; For the overall structural quality. For the first The density of the phase material, This is a pre-set upper limit for the overall structural mass. For the first The model represents the volume of a finite element; two materials are used in the model to simulate the main spars and ribs of the wing, respectively. The constraint represents the control of the hole size for the rib material, with the superscript 1 indicating the rib material; the constraint represents the volume fraction of material j. The upper limit of the set body fraction ratio, For the overall structural volume, Constraints for controlling bending deformation.

[0132] Step 2: First, based on the multi-material topology optimization model of the high aspect ratio airfoil, the airfoil is meshed using the finite element method (FEM) to generate several finite element elements. The original variable field of the high aspect ratio airfoil is then established based on the number of finite element elements and the topology optimization model. The principles and processes of finite element meshing and establishing the original variable field are well-known in this field and will not be elaborated here.

[0133] In this design, each finite element element of the high aspect ratio wing contains two candidate materials, representing the wing spars and ribs, respectively. Different variables in the original variable field correspond to different finite element elements of the wing, and these variables describe whether either of the two material phases exists in the corresponding finite element element.

[0134] Step 3: Generate filtered and projected variable fields based on the original variable field and density filtering and Heaviside projection. Based on the original variable field obtained in Step 2, the filtered and projected variable fields are established sequentially. The filtered variable field for each phase material is obtained by multiplying the original variable field of that phase material by the filtering matrix. The projected variables in the projected variable field of each phase material correspond one-to-one with the filtered variables in the filtered variable field of that phase material. The projected variable field is obtained by processing the filtered variable field using a projection function. Specifically, this includes:

[0135] Step 3.1: Establishing the filtered variable field

[0136] The filter radius is set in advance, and the filtered variable of any element in the original variable field is calculated The filtered variable field is output The original variable field and the filtered variable field The relationship is:

[0137] (3)

[0138] Wherein, , are the number of materials involved in the design and the number of units that can be designed, respectively, wherein , including the main beam material and the wing rib material; is the original variable field, ; represents the design variable related to the kth unit and the jth material; is the filtered matrix, ; represents the distance weight of the kth unit with respect to the jth unit. The expression of

[0139] is:

[0140] (4)

[0141] Wherein, is the volume of the kth unit, is the volume of the jth unit, represents the distance weight of the kth unit with respect to the jth unit. The definition of

[0142] (5)

[0143] Wherein, is the filter radius given in advance, is the distance between the center of the kth unit and the jth unit.

[0144] Step 3.2: Establishing the projection variable field

[0145] The filtered field is operated by using the Heaviside function to obtain the projection variable field , and its expression is:​​​​​​

[0146] (6)

[0147] wherein, is the projection parameter, which is used to control the steepness of the projection. In the optimization process, an incremental adjustment strategy is adopted to dynamically increase, in the present invention, the projection parameter is taken as a fixed value of 1 in the first 30 iteration steps, and thereafter the projection parameter is increased by 1.5 times every 10 steps until it reaches the maximum iteration step of 128. is the projection threshold value, which is taken as 0.5. is the projection variable of the first phase solid material of the first element, is the projection variable of the first phase solid material of the first element.

[0148] Step 4: Establish a UMMI material interpolation model according to the projection variable field obtained in step 3 to characterize the mapping relationship between the final projection variable value and the material properties of the element. And based on the best effort UMMI material interpolation model, a finite element commercial solver is called to perform a finite element analysis considering nonlinear geometry to obtain the deformation displacement of each finite element in the high-aspect-ratio wing. In the established UMMI multi-material interpolation model, the UMMI interpolation of Young's modulus is expressed in the form of the weight sum of each phase material:

[0149] (7)

[0150] wherein, is the weight function of the projection variable , is the Young's modulus of the first phase material, is the number of materials involved in the optimization, in the present invention the value of is 2, including the main beam material and the rib material.

[0151] The weight function is defined as:

[0152] (8)

[0153] wherein p is the penalty factor, in the present invention the value is taken as 3. is the projection variable of the first phase solid material of the first element, is the projection variable of the first phase solid material of the first element.

[0154] In step 4, the specific process of finite element analysis based on the UMMI material interpolation model is as follows:

[0155] Firstly, the material properties of each element are calculated according to the variables of each finite element unit by using the UMMI material interpolation model. Then, the obtained material properties are applied to the corresponding finite element unit to form a finite element analysis model reflecting the optimization design. Finally, the finite element analysis model is imported into a commercial finite element solver to perform finite element iterative optimization solving, and the deformation displacement of each finite element unit is extracted from the finite element analysis result file.

[0156] Step 5: According to the deformation displacement of each finite element unit, the data information transmission strategy under aerodynamic-structure coupling, the structure displacement transmission strategy, the aerodynamic load transmission strategy and the aerodynamic load updating strategy are established; and the aerodynamic-structure analysis is performed on the wing, and according to the analysis result, it is judged whether the aerodynamic load of each unit needs to be updated, if yes, the aerodynamic load is updated and returns to step 2; otherwise, step 6 is executed. The specific construction process of each strategy is as follows:

[0157] (1) Data information transmission strategy under aerodynamic-structure coupling

[0158] The weak coupling data transmission method is adopted in the present application, which mainly includes the following five stages, as shown in the figure. Figure 4

[0159] First stage: finite element solving. On the basis of the original aerodynamic load, the commercial software solver is called to calculate the deformation of the wing structure, and the displacement of each node of the wing skin in each direction is exported.

[0160] Second stage: calculate the deformed point cloud. The coordinates of the skin nodes before the deformation of the wing structure are added to the corresponding displacement to obtain the point cloud of the deformed skin.

[0161] Third stage: reverse modeling. Based on the spatial coordinate information of the deformed point cloud of the skin, the convergence body interpolation, creative modeling, inclusive body and split body series functions in the UG commercial software are used to realize the reverse modeling of the deformed wing.

[0162] Fourth stage: calculate the aerodynamic load after deformation. The deformed wing geometry model is imported into Fluent for CFD simulation, the aerodynamic load after deformation of the wing is calculated, and the load file is exported.

[0163] Fifth stage: the aerodynamic load calculated after deformation is mapped to the finite element model of the wing in optimization to realize the updating of the aerodynamic load.

[0164] (2) Structure displacement transmission strategy ​

[0165] After determining the need to update the aerodynamic load of the wing structure, the initial coordinates of the wing structure skin and the deformation in each direction are first extracted, and the two are added to obtain the point cloud distribution after the deformation of the skin. Then, the point cloud is fitted into a convergent body composed of numerous small planes using UG, and the geometric framework is generated in the UG creative modeling application module. The fitting accuracy of the geometric framework and the deformed finite element model is controlled by the reverse plane element size, and the smaller the value, the closer the reverse generated geometric framework to the deformed finite element model. Then, the obtained geometric framework is further processed by the containing body generation and split body to obtain the geometric model of the deformed wing. Finally, the geometric model is imported into the Fluent module in Workbench to calculate the aerodynamic load of the deformed wing skin, as shown in Figure 5 .

[0166] (3) Aerodynamic load transfer strategy

[0167] In order to transfer the aerodynamic load of the deformed wing calculated by the Fluent module in the commercial software Workbench to the structural finite element model in topology optimization, an aerodynamic load transfer module based on spatial search is established. Since the aerodynamic mesh used for aerodynamic load calculation and the nodes and meshes of the structural mesh in topology optimization do not have a one-to-one mapping relationship, and their numbering does not match, and the structural mesh in optimization is always in the initial state of deformation, the deformed point cloud model is used for transition processing of aerodynamic load, which mainly includes the following three steps, as shown in Figure 6 .

[0168] 1) Aerodynamic load spatial search mapping. The spatial search method is used to traverse the output aerodynamic load file one by one to find the nearest aerodynamic load point to each point in the point cloud model, and the pressure value at the aerodynamic load point is transferred to the corresponding point in the point cloud model. In addition, since the number of aerodynamic meshes divided each time when calculating the aerodynamic load is inconsistent, the total amount of data in the output aerodynamic load file is not consistent. In order to keep the total lift in the calculated aerodynamic load consistent, each interpolated load is multiplied by a lift coefficient, which is defined as the ratio of the total lift in the aerodynamic load file to the interpolated total lift.

[0169] 2) Point cloud load transfer to structural finite element model. Since the deformed point cloud model is composed of the deformation of the skin nodes in each direction in the structural finite element model, and the point number in the point cloud model and the node number in the structural skin model keep one-to-one correspondence, the coordinate value in the deformed point cloud load file is subtracted from the deformation of the corresponding node on the skin, that is, the aerodynamic load file of the finite element model in optimization is obtained.

[0170] 3) Update and optimize the loads on the wing finite element model. Since pressure loads are volume loads, they need to be applied to the elements. However, the aerodynamic load file of the finite element model obtained above contains point load information. Therefore, the mean value method is used, i.e., the load is applied to the elements. The pressure value is the average of the pressure values ​​at its four nodes.

[0171] (4) Aerodynamic load update strategy

[0172] In the process of optimizing wing topology considering aerodynamic-structural coupling, if the aerodynamic loads are updated in real time for each iteration, the optimization problem may be frequently changed, making convergence difficult. However, if the number of iterations in between is too large, the load span before and after the update will be large, making the optimization prone to abrupt changes. Therefore, establishing a suitable aerodynamic load update strategy is crucial.

[0173] This invention employs an aerodynamic load update strategy based on deformation control, specifically, determining whether to update the aerodynamic load by measuring the displacement change VD in the y-direction of the wingtip section center node. This is defined as...

[0174] (9)

[0175] In the formula This represents the displacement in the y-direction of the center node of the wingtip section in the current iteration step. The reference object for the y-direction displacement of the center node of the wingtip section is used. During the optimization process, the aerodynamic load is updated each time. Then update to the current iteration step. Values, where the aerodynamic loads are updated for the first time. This represents the displacement value in the y-direction of the center node of the wingtip section at step 0. To pre-introduce an upper limit for displacement change, this invention, after numerous numerical trials in the calculation examples, [details omitted]. Set it to 30%.

[0176] Step 6: Establish material void orientation size control constraints using local cylinder search and the introduction of local material rate method; establish local bending deformation control constraints for high aspect ratio airfoils using element nodal displacement. Among these,

[0177] The process for establishing constraints for the directional size control of material pores is as follows:

[0178] Based on the characteristics of traditional wing spars-rib structures, a maximum aperture member size control constraint is adopted using a local cylindrical search along the wingspan direction, such as... Figure 7 As shown.

[0179] Regarding the rib material (material 1), in the first Introducing material rate oriented towards rib material within the local search region of the element , where

[0180] (10)

[0181] where, is the volume of the element; is a penalty factor, which is set to 3. is the local search region centered at the element, is the element with respect to the projected variable of the rib material.

[0182] Equation (10) indicates that there is a certain material ratio only when the variables related to the rib material in the element are not equal to 0. Introducing the material ratio lower limit , the maximum hole size control constraint expression of the material 1 (rib material) is:

[0183] (11)

[0184] where, is the local search region in the element.

[0185] Introducing the p-mean aggregation function to aggregate all local constraints into a single global constraint, according to the definition of the aggregation function, the function is constructed, whose expression is defined as

[0186] (12)

[0187] That is, the global expression of the directional size control constraint is:

[0188] (13)

[0189] where, the expression is:

[0190] (14)

[0191] where, is the number of elements, is the penalty factor of the aggregation function, which is set to 80.

[0192] The process of local bending deformation control constraint of high aspect ratio wing is:

[0193] According to the characteristics of the traditional wing beam-rib structure form, the maximum hole member size control constraint of local cylindrical search along the length direction of the wingspan is adopted, such as​Figure 8 As shown. To precisely control wing bending deformation, multiple nodes evenly positioned on the wingtip section are selected as observation nodes, such as... Figure 8 As shown, the displacement is extracted and the maximum nodal displacement control constraint model in the y-direction is established as follows:

[0194] (15)

[0195] in, Indicates the node number is The observation nodes are at Displacement in the direction, This is a preset upper limit constraint. Furthermore, the displacements of the above observation nodes can be derived from the overall displacement vector of the wing structure. Extracted from, i.e. It can be written in the following form:

[0196] (16)

[0197] In the formula, the displacement extraction vector Defined as:

[0198] (17)

[0199] That is, only the node numbered in this vector is 0. The observation nodes are at The directional component is 1, and all other terms are 0. Since the max function in equation (15) is discontinuous, it is non-differentiable when calculating the sensitivity of this constraint. Therefore, in this invention, a continuous p-mean condensation function is used to approximate the max function. Based on the definition of the p-mean condensation function, variables are constructed... Its expression is

[0200] (18)

[0201] Right now

[0202] (19)

[0203] In the formula, The number of observation nodes, Let be the exponent of the condensation function. Therefore, the final control constraint for wing bending deformation can be written in the following form:

[0204] (20)

[0205] Step 7: First, determine the constraint function based on the material hole orientation size control constraint and the local bending deformation control constraint of the large aspect ratio wing obtained in Step 6.

[0206] According to the deformation displacement of each finite element unit calculated in step 4, the sensitivity of the objective function and the constraint function to the design variable is calculated respectively. According to the sensitivity calculation result, the design variable of the original variable field is iteratively optimized and solved under multiple constraint conditions by using a gradient-based optimization algorithm, so that the objective function reaches the optimal value under the premise of meeting each constraint condition. It should be noted that the optimal value in the application refers to the minimum value of the objective function (end-point compliance of the high-aspect-ratio wing) that reaches stable convergence under the premise of meeting other constraints, wherein the convergence criterion is that the optimization target value of the previous and next two iteration steps changes by not more than 0.1%.

[0207] The plurality of design constraints include mass constraint, volume constraint, material hole directional size control constraint, wing bending deformation control constraint and design variable constraint. The mass constraint is that the mass of the wing calculated based on the projection variable field is less than or equal to the pre-set upper limit of the mass. The volume constraint is that the volume of the same phase material calculated based on the projection variable of each unit corresponding to the same phase material is less than or equal to the pre-set upper limit of the volume. The design variable constraint is that the variable in the projection variable field is greater than zero and less than or equal to 1. The specific process is as follows:

[0208] The sensitivity of the optimization response in the optimization model is derived, and then the derived sensitivity is transmitted to the optimization algorithm. Since the derivation process between different constraints is similar, only the sensitivity of the optimization target end-point compliance is derived here.

[0209] The present application uses the adjoint method and Taylor formula to derive the approximate expression of the end-point compliance of the structure based on the node displacement increment. According to the definition, the adjoint vector is introduced The end-point compliance of the wing is obtained

[0210] (21)

[0211] In the formula, is the end-point compliance of the wing, is the external force of the unit node, is the final displacement obtained by solving the balance equation of the finite element, is the residual force vector. The aerodynamic load of the present application is not related to the design variable, therefore, the sensitivity of the end-point compliance to the projection variable can be expressed as:

[0212] (22)

[0213] In the formula, is the tangent stiffness matrix. can be further expressed as

[0214] (23)

[0215] Eliminating unknown terms by adjoint method , so set

[0216] (24)

[0217] To get the approximate expression of sensitivity, here we add a new case to the original case 1, which is case 2. The superposition of case 1 and case 2 is case 3, which is where is usually a small value, which is taken as in this invention. Therefore, the equilibrium equations of the two cases are respectively

[0218] (25)

[0219] (26)

[0220] where is the displacement vector of case 1, is the displacement vector of case 2. The above equation is rearranged as

[0221] (27)

[0222] The Taylor expansion of the internal force of the second case is obtained as follows:

[0223] (28)

[0224] where is the high-order term of the displacement difference of the two cases, which is omitted. The above equation is substituted into equation (3-34) to obtain the approximate expression of

[0225] (29)

[0226] where is the displacement difference of the two cases. In addition, the internal node force of the element can be approximately expressed as:

[0227] (30)

[0228] where is the difference of strain energy of the element in the two cases, which can be directly obtained in the commercial software solver. Combined with the UMMI interpolation model, we finally get

[0229] (31)

[0230] where​​

[0231] (32)

[0232] Figure 2 and Figure 3 The constraints and load conditions of the wing structure to be optimized in this embodiment are shown, wherein the wing span is 8m and the wing chord length is 1m. The topology optimization of this wing is performed using the method of this invention, and the optimized configuration is shown below. Figure 9 As shown in Figure 9b, and compared with the optimization results without considering aerodynamic / structural coupling ( Figure 9 a) For comparison, material 1 in the figure represents the wing rib material, and material 2 represents the main sparsity material. As can be seen from the figure, after considering aerodynamic / structural coupling, the optimized configuration shows that the front main sparsity gradually tapers along the wingspan and its position moves from the center of the leading edge camera wing. Some staggered ribs are distributed at the end of the rear main sparsity. However, materials with higher Young's modulus are more concentrated at the wing root, and the front main sparsity at the wingtip also becomes thinner. Figure 9 The features within the red circle in section b are shown. Furthermore, several inclined ribs are distributed near the wingtip of the main sparsity, and the spatial arrangement of these ribs is staggered from the optimized configuration without considering aerodynamic-structural coupling. Under the influence of large deformation, areas with larger aerodynamic loads also shift towards the wing root as the deformation increases. However, in the optimized configuration considering aerodynamic-structural coupling in this invention, the movement of the strong material towards the wing root exhibits the same trend as the change in aerodynamic load distribution. This sufficiently demonstrates that the optimization results considering aerodynamic-structural coupling are more accurate and have higher practical engineering value.

[0233] It should be noted that in the method of this invention, the gradient-based GCMMA optimization algorithm is selected, and a self-developed topology optimization program is used for optimization iteration. Specifically, other gradient-based optimization algorithms such as ConLin, GCM, MDPA, SLP, and QP can also achieve optimization iteration. Furthermore, some other optimization methods such as optimization criteria methods, mathematical programming methods, and asymptotic structure optimization methods can also be used for the optimization design of the method of this invention.

[0234] Step 8: Redesign the optimized configuration using a truss structure.

[0235] The topology optimization results considering aerodynamic coupling are redesigned using a truss structure scheme with greater lightweight potential. Figures 10-11The reconstruction design scheme is given, in which the structure mode of the front main beam 3 and the rear main beam 4 envelope, the auxiliary beam 5 and the intermediate several cable-stayed rods is used to represent the material distribution of the full entity main beam in the optimization result, wherein the cable-stayed rods include several main beam cable-stayed rods 7 and several wing rib cable-stayed rods 9, the main beam cable-stayed rods 9 are close to the wing tip; several main beam reinforcing rods 8 and main beam auxiliary beams 11 are also arranged on the main beam. The wing ribs 6 refer to the spatial distribution in the optimization result, the interval between adjacent wing ribs is set to 0.4 meters, and a wing rib 6 is added at the cross section of the wing root 1 and the wing tip 2 respectively. The structural form of the wing rib 6 is selected as the conventional layout, and the stability of the wing rib is enhanced through the wing rib cable-stayed rod 9 and the vertical wing rib reinforcing rod 10. In addition, the length of the auxiliary beam 5 at the trailing edge is given according to the optimization result. Finally, the overall mass of the structure of the reconstruction design scheme according to the optimization configuration considering the aerodynamic-structural coupling and the bending control constraint is 16.948 kg. The specific size parameters in the reconstruction scheme are given in detail in Table 1, wherein the main beam and the main beam reinforcing rod are divided into S1-S10 for a total of 10 segments from the wing root to the wing tip.

[0236] Table 1

[0237]

[0238] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited to this, any skilled person in the art can easily think of various equivalent modifications or replacements within the technical range disclosed by the present application, and these modifications or replacements should be covered within the protection scope of the present application.

Claims

1. A method for topology optimization of a high-aspect-ratio wing considering aero-structural coupling, characterized in that, The method comprises the following steps: Step 1: determining an optimization objective and an optimization constraint of a high-aspect-ratio wing multi-material structure topology optimization problem, and establishing a high-aspect-ratio wing multi-material topology optimization model considering aerodynamic-structure coupling according to the optimization objective and the optimization constraint; wherein the optimization objective is to minimize the wing end compliance; the optimization constraint comprises wing rib directional size control constraint, local bending deformation control, wing overall mass constraint and single material volume fraction constraint; The high-aspect-ratio wing multi-material topology optimization model considering aerodynamic-structure coupling is as follows: In the formula: Represents the set of design variables. For design variables, representing the first... The finite element and the first To avoid singularities in the structural stiffness matrix caused by design variables at 0 during finite element analysis, a lower limit for design variables is introduced. ;when When, it indicates that the unit is entirely composed of the j-th phase material, when When the value is between the lower limit of the design variable and 1, it indicates that the unit is composed of the first... The degree of phase material composition; in order to better determine the first The UMMI multi-material interpolation model uses an exponential function penalty to determine which phase of material each finite element is ultimately composed of, ensuring that the final variable field converges to the lower limit of the design variables or 1; the initial values ​​of all design variables are set to 0.

1. For the first Projected variable field of phase material, For the first Phase solid material filtering variable field, For the first The original variable field of the phase material, Indicates the first Phase solid materials in the first Projection variables of the unit; Represents the total number of design domain units; To optimize the objective, This refers to the structural displacement vector in the Nth step of the iterative solution. For the internal forces at the unit nodes, To employ an update strategy based on wingtip deformation to participate in the updating of aerodynamic loads, The residual force vector; For the Heaviside projection function, This is a density filtering matrix; For the overall structural quality. For the first The density of the phase material, This is a pre-set upper limit for the overall structural mass. For the first The model represents the volume of a finite element; two materials are used in the model to simulate the main spars and ribs of the wing, respectively. Represents the constraint on the hole size control of the rib material; represents the volume fraction constraint of material j. The upper limit of the set body fraction ratio, For the overall structural volume, Constraints for controlling bending deformation; Step 2: according to the high-aspect-ratio wing multi-material topology optimization model, finite element meshing is performed on the high-aspect-ratio wing to generate a plurality of finite element units; finite element analysis is performed on the high-aspect-ratio wing multi-material topology optimization model according to the high-aspect-ratio wing structure boundary constraint and initial load condition, and an original variable field of the high-aspect-ratio wing is obtained; Wherein the high-aspect-ratio wing is composed of two-phase materials; different variables in the original variable field correspond to different finite element units of the wing, and the variables are used to describe whether any phase material of the two-phase materials exists in the corresponding finite element unit; Step 3: the original variable field is sequentially subjected to density filtering and Heaviside projection operation to generate a filtered variable field and a projected variable field; wherein: The filtered variable field of each phase material is obtained from the original variable field of the phase material and a set filtering matrix; the projected variable field is obtained by projecting the filtered variable field through a projection function, and the projection variable in the projected variable field of each phase material corresponds to the filtered variable in the filtered variable field of the phase material one by one; Step 4: a UMMI material interpolation model is established according to the projected variable field obtained in step 3 to represent the mapping relationship between the projected variable field and the unit material property; based on the UMMI material interpolation model, a finite element analysis considering nonlinear geometry is performed by calling a commercial finite element solver to obtain the deformation displacement of each finite element unit in the high-aspect-ratio wing; Step 5: according to the deformation displacement of each finite element unit, a data information transmission strategy, a structure displacement transmission strategy, an aerodynamic load transmission strategy and an aerodynamic load updating strategy under aerodynamic-structure coupling are established; and aerodynamic-structure analysis is performed on the wing, and whether the aerodynamic load of each unit needs to be updated is judged according to the analysis result; if yes, the aerodynamic load is updated and the step 2 is returned; otherwise, step 6 is executed; Step 6: a local cylindrical search and a local material rate method are used to establish a material hole directional size control constraint; and a unit node displacement is used to establish a local bending deformation control constraint of the high-aspect-ratio wing; Step 7: according to the material hole directional size control constraint and the local bending deformation control constraint of the high-aspect-ratio wing, a constraint function is determined; According to the deformation displacement of each finite element unit, the sensitivity of the objective function and the constraint function to the design variable is calculated respectively; according to the sensitivity calculation result, the design variable in the original variable field is iteratively optimized and solved by using the gradient-based optimization algorithm under multiple constraint conditions, until the objective function reaches the optimal under the constraint condition, and then the optimal configuration of the high-aspect-ratio wing is obtained; Wherein, the multiple design constraints include mass constraint, volume constraint, material hole directional size control constraint, wing bending deformation control constraint and design variable constraint; Step 8: The optimal configuration of the high-aspect-ratio wing is reconfigured by using the truss structure for lightweight design.

2. The method of claim 1, wherein, The step 1 includes the following sub-steps: Step 1.1: Establishing a wing end point compliance expression, and taking the minimum wing end point compliance as the optimization objective; wherein, the wing end point compliance expression is: wherein Fextis an external force to the unit node force, u is the final displacement obtained by the nonlinear iterative solution; Step 1.2: Establishing a wing finite element analysis model considering nonlinear geometry, and performing finite element meshing on the wing to obtain a plurality of finite element units; wherein, the wing skin is simulated by shell elements, and the design domain is simulated by solid elements, each finite element unit is assigned with two candidate materials, representing the wing main beam and wing rib respectively; Step 1.3: Establishing constraint conditions, including: wing overall mass constraint, single material volume fraction constraint, directional size control constraint and bending deformation control constraint; Step 1.4: According to the optimization objective, the original design variable field of the wing and the constraint conditions, a high-aspect-ratio wing multi-material topology optimization model considering aerodynamic-structure coupling is constructed.

3. The method of topology optimization of a high-aspect-ratio wing according to claim 1 or 2, characterized in that, The step 3 includes the following sub-steps: Step 3.1: Set the filter radius in advance, calculate the filtered variable of any j-phase material in the original variable field, and output the filtered variable field : In the formula, , These represent the number of materials involved in the design and the number of designable units, respectively. The value of 2 represents the material of the main beam and the rib. For the original variable field, ; Indicates the first Unit and the Design variables related to phase materials; For the filtering matrix, ; Indicates the k-th unit about the th The weight coefficient of a unit is expressed as follows: wherein Vkis the volume of the kth cell, Vkis the volume of the kth cell, Vkis the volume of the kth cell, Dk,ξrepresents the distance weight of the kth cell with respect to the Dk,ξrepresents the distance weight of the kth cell with respect to the Dk,ξrepresents the distance weight of the kth cell with respect to the Dk,ξrepresents the distance weight of the kth cell with respect to the Definitions: wherein is a predefined filter radius, is the first unit and the first center distance of the cells; Step 3.2: Filtering the variable field with the Heaviside function Performing a projection operation to generate a projected variable field ; where the projected variable field is given by: where the projected variable of the i-th cell with respect to the j-th phase entity material is given by: In the formula, is a projection parameter, used to control the steepness of the projection, An incremental adjustment strategy is adopted to dynamically increase during the optimization process; is a projection threshold, which is 0.

5.

4. The method of claim 3, wherein, In the step S3.2, the increment adjustment strategy is adopted to dynamically adjust The process is as follows: At 1-30 iteration steps, projection parameters is taken as a fixed value 1; From step 31, every 10 iteration steps, the projection steepness control parameter is multiplied by 1.5, and this process is repeated until the maximum number of iteration steps is reached.

5. The method of topology optimization of a high-aspect-ratio wing according to claim 1 or 2, wherein, In the UMMI multi-material interpolation model established in the step 4, the UMMI interpolation expression of Young's modulus is expressed as the sum of the weights of each phase material: wherein is a weight function of the projection variable , is the Young's modulus of the th material, is the number of materials involved in the optimization, ; The weight function is defined as with the SIMP penalty function wherein is a penalty factor; is the unit with respect to the projected variable of the is the unit with respect to the projected variable of the 6. The method of claim 5, wherein, In the step 4, the specific process of finite element analysis based on the UMMI material interpolation model is: Using the UMMI material interpolation model, the material properties of each finite element unit are calculated according to the variables of each finite element unit; The obtained material properties are applied to the corresponding finite element units to form a finite element analysis model reflecting the optimization design; The finite element analysis model is imported into a commercial finite element solver to perform finite element iterative optimization solving, and a finite element analysis result file is output, from which the deformation displacement of each finite element unit is extracted.

7. The method of topology optimization of a high-aspect-ratio wing according to claim 1 or 2, wherein In the step 5, The data information transmission strategy under the aerodynamic-structure coupling is: Finite element solving: on the basis of the original aerodynamic load, a commercial solver is called to calculate the wing structure deformation, and the displacement of each node of the wing skin in each direction is output; Deformation point cloud calculation: the node coordinates of the wing skin before deformation are added to the corresponding displacement to obtain the point cloud of the deformed skin; Reverse modeling: based on the point cloud spatial coordinate information of the deformed skin, the convergence body interpolation, creative modeling, inclusive body and split body functions of UG commercial software are used to realize the reverse modeling of the deformed wing; Deformed aerodynamic load calculation: the deformed wing geometry model is imported into Fluent for CFD simulation to calculate the aerodynamic load of the deformed wing, and the load file is exported; Aerodynamic load transfer: mapping the deformed aerodynamic load to the optimized finite element model to realize the update of aerodynamic load; The structural displacement transfer strategy is: The initial coordinates of the wing skin and the deformation in each direction are extracted, and the two are added to obtain the point cloud distribution of the deformed skin; the point cloud is fitted into a convergent body composed of a large number of small planes by using UG commercial software, and the plane convergent body is reversely generated into a geometric framework in the UG creative modeling application module; the obtained geometric framework is sequentially subjected to an inclusive body generation and a split body processing to obtain a deformed wing geometric model; finally, the geometric model is imported into the Fluent module in Workbench to calculate the aerodynamic load of the deformed wing skin shape; The aerodynamic load transfer strategy is: Aerodynamic load space search mapping: a space search method is used to traverse the output aerodynamic load file one by one to find the nearest aerodynamic load point for each point in the point cloud model, and the pressure value at the aerodynamic load point is transferred to the corresponding point in the point cloud model; and each load after interpolation is multiplied by the lift coefficient, and the lift coefficient is the ratio of the total lift in the aerodynamic load file to the total lift after interpolation; Point cloud load to wing finite element model transfer: subtracting the deformation of the corresponding node on the skin from the coordinate value in the deformed point cloud load file to obtain the aerodynamic load file of the optimized wing finite element model; Updating the load of the optimized wing finite element model: The aerodynamic load updating strategy is: through the center node of the wing tip section The displacement change amount in the direction Judging whether to update the aerodynamic load, and the updating formula is: In the formula, is the displacement of the center node of the wing tip section in the y direction in the current iteration step; is the displacement of the center node of the wing tip section in the y direction in the first iteration step; is the reference object of the displacement in the y direction, and the aerodynamic load is updated every time in the optimization process, is updated to the displacement of the center node of the wing tip section in the y direction in the current iteration step, and the aerodynamic load is updated every time in the optimization process, is the value of the displacement of the center node of the wing tip section in the y direction in the first iteration step, and the aerodynamic load is updated every time in the optimization process, is the displacement of the center node of the wing tip section in the y direction in the first iteration step; is the upper limit of the pre-introduced displacement change, and the value is set to 30%.

8. The method of claim 7, wherein, In step 6, the material hole directional size control constraint establishment process is: According to the characteristics of the traditional wing beam-rib structure form, the wing rib material is introduced in the local search area of the unit The material rate of the wing rib material is introduced in the local search area of the unit : ； In the formula, Indicates the first The volume of the unit; As a unit The local search area centered on the search center. Let k be the projection variable of the element with respect to the rib material; This is the penalty factor, and its value is 3. Introducing a material rate lower bound The maximum hole size control constraint expression for the wing rib material is obtained as In the formula, is the first Introducing a material rate oriented to the rib material within the local search region of the cell. The p-mean condensing function is introduced to condense all local constraints into a single global constraint. According to the definition of the condensing function, the function is defined as Then, the global expression for establishing the directional size control constraint is: wherein The expression is In the formula, is the number of units, is a penalty factor of the cohesion function, and has a value of 80; The establishment process of the bending deformation control constraint is: Selecting multiple nodes with uniform positions on the wing tip section as observation nodes, extracting the displacement of each observation node and establishing Maximum nodal displacement control constraint model in the direction: In the formula, represents the displacement of the observation node numbered in the direction , and is a preset upper limit of the constraint. The displacement of each observation point in the direction of the line of sight is expressed as: the displacement of each observation point in the direction of the line of sight is expressed as: the displacement of each observation point in the wherein is the overall displacement vector of the wing structure, is the displacement extraction vector; Definitions: According to the definition of the p-mean cohesion function, the variable is constructed, whose expression is The expression of the maximum nodal displacement control constraint model in the direction is rewritten as: The expression of the maximum nodal displacement control constraint model in the direction is rewritten as: wherein is the number of observation nodes, is the index of the cohesion function; then the bending deformation control constraint is obtained .

9. The method of claim 8, wherein, In step 7, The mass constraint: the mass of the high-aspect-ratio wing calculated based on the projection variable field is less than or equal to the pre-set upper limit of the mass; The volume constraint: the volume of the same phase material calculated based on the projection variable of each finite element unit corresponding to the same phase material is less than or equal to the pre-set upper limit of the volume; The design variable constraint: the variables in the projection variable field are limited in the open interval (0, 1].

Images