Multidisciplinary Coupling Optimization Design Method for Integrated Joined-Wing Layout

Abstract

The present invention provides a multidisciplinary coupling optimization design method for an integrated joined-wing layout. First, determine the design objectives, design conditions, and design variables for the multidisciplinary optimization of the joined-wing layout. Then, the design variables are divided into global design variables and local design variables. Next, with high lift-to-drag ratio, low structural mass, and low electromagnetic shielding rate as the system-level design objectives, and the global design variables as the design variables, fix the local design variables and conduct system-level optimization. After that, with high lift-to-drag ratio, low structural mass, and low electromagnetic shielding rate as the subsystem-level design objectives respectively, and the local design variables of the corresponding disciplines as the design variables, fix the global design variables and the local design variables of other disciplines and conduct subsystem-level optimization for each discipline. Finally, determine whether the convergence condition is reached. If not, continue with system-level optimization and subsystem-level optimization until convergence. The present invention realizes the multidisciplinary coupling optimization design of aerodynamics / structure / electromagnetics for the integrated joined-wing layout.

Description

Technical Field

[0001] The present invention relates to the field of aircraft design, and particularly to a multidisciplinary coupling optimization design method for an integrated joined-wing layout. Background Technique

[0002] For the integrated design of an ultra-long-endurance loitering unmanned aerial vehicle platform with a joined-wing layout, it is often required to have a high lift-to-drag ratio, a small structural mass, and a small electromagnetic shielding ratio, which is a typical multidisciplinary coupling optimization design problem of aerodynamics / structure / electromagnetics. At the same time, the design objectives of multiple disciplines bring high-dimensional design variables. Therefore, the integrated design of an ultra-long-endurance loitering unmanned aerial vehicle platform with a joined-wing layout is a typical high-dimensional multi-objective problem.

[0003] Currently, conventional multidisciplinary optimization methods are divided into two categories: monolithic and distributed. The monolithic method regards the optimization problem as a whole for solution, such as the multidisciplinary feasible method and the single-discipline feasible method. These methods have a high dimension of design variables and are prone to the problem of dimensional disaster, resulting in low design efficiency. The distributed method decomposes the optimization problem into multiple levels for analysis, including the collaborative optimization method and the parallel subspace optimization method. Among them, the collaborative optimization method decouples the subsystems by introducing auxiliary design variables and consistency constraint conditions, but it has the problems of insufficient convergence and easy premature convergence. The parallel subspace optimization method maintains multidisciplinary consistency through the global sensitivity equation, but it needs to solve gradient information in the optimization, and there are problems of difficult gradient solution and large computational amount.

[0004] In addition, when performing multidisciplinary optimization of aerodynamics / structure / electromagnetics for a joined-wing layout, it is necessary to continuously calculate the structural characteristics of the current joined-wing layout scheme. However, the structural characteristic analysis methods of traditional configuration aircraft, including the piecewise linear stiffness method and the solid element finite element method, are not applicable to the integrated joined-wing configuration design. The main reason is that the piecewise linear stiffness method cannot well adapt to the special stiffness configuration of the joined wing and has problems of insufficient accuracy, while the solid element finite element method requires a large number of elements and nodes, with a long analysis time and a large computational cost, resulting in low optimization design efficiency. Summary of the Invention

[0005] To solve the above problems, the present invention proposes a multidisciplinary coupling optimization design method for an integrated joined-wing layout. By classifying design variables and adopting a hierarchical design method of multidisciplinary system-level optimization and subsystem-level optimization, it can reasonably optimize global design variables and local design variables and significantly reduce the subsystem analysis time, realizing the efficient multidisciplinary optimization design of the joined-wing layout. In addition, for the integrated joined-wing layout, the structural target calculation adopts the finite element equivalent beam method based on engineering beam theory. On the one hand, it combines the classical structural force transmission model of engineering beam theory, and on the other hand, it adds a finite element beam model, which can be used to analyze the special stiffness configuration of the joined wing. In the electromagnetic discipline optimization, an improved rapid electromagnetic shielding evaluation method is adopted, which can quickly calculate the electromagnetic shielding rate of the aircraft antenna array in the integrated joined-wing layout, so as to provide an efficient and accurate evaluation method for the structural design and aerodynamic / structural / electromagnetic multidisciplinary design of the integrated joined-wing layout.

[0006] The technical solution of the present invention is specifically as follows and includes the steps:

[0007] Step 1: Establish the initial configuration of the aircraft with an integrated joined-wing layout and parameterize the configuration; determine the design objectives, design conditions, and design variables for the multidisciplinary optimization of the aircraft with an integrated joined-wing layout, as well as the design space range of each design variable.

[0008] The design objectives include a high lift-to-drag ratio, low structural mass, and low electromagnetic shielding rate; the design variables include aerodynamic discipline design variables, structural discipline design variables, and electromagnetic discipline design variables.

[0009] Step 2: Divide the design variables determined in Step 1 into global design variables XS and local design variables XL.

[0010] Step 3: Conduct system-level optimization.

[0011] Taking a high lift-to-drag ratio, low structural mass, and low electromagnetic shielding rate as the system-level design objectives, taking the global design variables as the system-level optimization design variables, and fixing the local design variables as constants; using the CFD method for aerodynamic target analysis, using the equivalent beam method based on engineering beam theory for structural target analysis, and using the improved rapid electromagnetic shielding evaluation method for electromagnetic target analysis; using the genetic algorithm for optimization to obtain the global optimal design result.

[0012] Step 4: Conduct subsystem-level optimization respectively, including aerodynamic discipline optimization, structural discipline optimization, and electromagnetic discipline optimization.

[0013] The optimization of the aerodynamic discipline takes the high lift-to-drag ratio as the optimization design goal of the aerodynamic discipline, uses the local design variables in the aerodynamic discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the structural discipline design variables, and the local design variables in the electromagnetic discipline design variables as constants; uses the CFD method for aerodynamic target analysis, and uses the genetic algorithm for optimization to obtain the local optimal design result of the aerodynamic discipline;

[0014] The optimization of the structural discipline takes the low structural mass as the optimization design goal of the structural discipline, uses the local design variables in the structural discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic discipline design variables after aerodynamic discipline optimization, and the local design variables in the electromagnetic discipline design variables as constants; combines the aerodynamic loads obtained by using the CFD method for aerodynamic target analysis in the aerodynamic discipline optimization process, uses the equivalent beam method based on the engineering beam theory for structural target analysis, and uses the genetic algorithm for optimization to obtain the local optimal design result of the structural discipline;

[0015] The optimization of the electromagnetic discipline takes the low electromagnetic shielding rate as the design goal of the electromagnetic discipline, uses the local design variables in the electromagnetic discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic discipline design variables after aerodynamic discipline optimization, and the local design variables in the structural discipline design variables after structural discipline optimization as constants; uses the improved rapid electromagnetic shielding evaluation method for electromagnetic target analysis, and uses the genetic algorithm for optimization to obtain the local optimal design result of the electromagnetic discipline;

[0016] Step 5: Judge whether the convergence end condition is reached. If not, return to step 3 to continue the system-level optimization and subsystem-level optimization until the convergence end condition is reached, end the optimization, and obtain the multidisciplinary coupling optimization design result of the integrated joined-wing layout.

[0017] Furthermore, in step 1, the configuration is parameterized, and the determined design variables are as follows: The aerodynamic discipline design variables include the two-dimensional airfoil design parameters of the front wing and the rear wing and the three-dimensional shape layout parameters of the joined wing;

[0018] The two-dimensional airfoil adopts the CST parameterization method, and the airfoil function is expressed as the product of the class function and the shape function, where the coefficients of the shape function are the two-dimensional airfoil design parameters;

[0019] The three-dimensional shape layout parameters include:

[0020] The span of the front wing 、The chord length of the front wing root 、The taper ratio of the front wing root 、The twist angle of the front wing tip 、The sweep angle of the front wing 、Anhedral angle of the front wing ; Span of the rear wing 、Chord length at the root of the rear wing 、Tip-to-root ratio of the rear wing 、Twist angle at the tip of the rear wing 、Swept angle of the rear wing 、Dihedral angle of the rear wing ; Position of the leading edge of the front wing ,Position of the leading edge of the rear wing ;Height ratio of the end plate ;

[0021] Among them, is the height difference between the front wing and the rear wing at the wing root, is the height difference between the front wing and the rear wing at the end plate;

[0022] The design variables of the structural discipline include structural dimension parameters and structural layout parameters;

[0023] Among them, the structural dimension parameters include: Thickness of the web of the first front wing beam: ,Thickness of the web of the second front wing beam: ;Area of the upper flange of the first front wing beam: ,Area of the lower flange of the first front wing beam: ;Area of the upper flange of the second front wing beam: ,Area of the lower flange of the second front wing beam: ;Thickness of the front wing rib: ;Thickness of the front wing skin: ;Thickness of the web of the first rear wing beam: ,Thickness of the web of the second rear wing beam: ;Area of the upper flange of the first rear wing beam: ,Area of the lower flange of the first rear wing beam: ;Area of the upper flange of the second rear wing beam: ,Area of the lower flange of the second rear wing beam: ;Thickness of the rear wing rib: ;Thickness of the rear wing skin: ;

[0024] The structural layout parameters include: Position of the first front wing beam: ,Position of the second front wing beam: ;Spacing of the front wing ribs: ;Position of the first rear wing beam: ,Position of the second rear wing beam: ;Spacing of the rear wing ribs: ;

[0025] The design variables of the electromagnetic discipline include antenna spacing ,Suppose there are m antennas arranged on the side of the fuselage, then ; among them, is the distance between the leading edge of the front wing and the first antenna, is the distance between the first antenna and the second antenna, and so on, is the distance between the (m - 1)-th antenna and the m-th antenna.

[0026] Furthermore, the global design variable XS in step 2 includes all the design parameters of the two-dimensional airfoil in the aerodynamic discipline design variables and the following parameters in the three-dimensional layout parameters of the joined wing:

[0027] Front wingspan , front wing root chord length , front wing root taper ratio , front wing sweep angle ; rear wingspan , rear wing root chord length , rear wing root taper ratio , rear wing forward sweep angle ; front wing leading edge position , rear wing leading edge position ; end plate height ratio ;

[0028] The local design variable XL in step 2 includes the local design variables of the aerodynamic discipline , the local design variables of the structural discipline , the local design variables of the electromagnetic discipline ;

[0029] Among them, the local design variables of the aerodynamic discipline include the front wing tip twist angle , front wing dihedral angle in the three-dimensional layout parameters of the joined wing; rear wing tip twist angle , rear wing anhedral angle

[0030] The local design variables of the structural discipline include all the structural dimension parameters and structural layout parameters;

[0031] The local design variables of the electromagnetic discipline include all the antenna spacings .

[0032] Furthermore, the mathematical model of the system-level optimization in step 3 is specifically as follows:

[0033]

[0034]

[0035]

[0036]

[0037] Among them, is the weight coefficient of the lift-drag ratio optimization target, is the weight coefficient of the structural mass optimization target, is the weight coefficient of the electromagnetic shielding rate target; is the lift-drag ratio, is the structural mass, is the electromagnetic shielding rate, is the lift-drag ratio of the initial configuration, is the structural mass of the initial configuration, is the electromagnetic shielding rate of the initial configuration; The constraint conditions include aerodynamic discipline constraints, structural discipline constraints, and electromagnetic discipline constraints. Specifically:

[0038] Aerodynamic discipline constraints: The cruise lift coefficient is not less than the set value of the lift coefficient , and the change in the pitching moment coefficient is not greater than the set value of the moment change ;

[0039] Structural discipline constraints: The cruise structural equivalent stress is not greater than the allowable stress , the maximum wing deformation does not exceed the set value of the deformation amount , and the wing torsional deformation amount in the yaw control state does not exceed the set value of the torsion angle ;

[0040] Electromagnetic discipline constraints: The electromagnetic shielding rate is not greater than the electromagnetic shielding rate of the initial configuration .

[0041] Furthermore, the specific mathematical model for aerodynamic discipline optimization in step 4 is as follows:

[0042]

[0043]

[0044] Among them, is the lift-drag ratio; The constraint conditions are specifically: The cruise lift coefficient is not less than the set value of the lift coefficient , and the change in the pitching moment coefficient is not greater than the set value of the moment change ;

[0045] The specific mathematical model for structural discipline optimization is as follows:

[0046]

[0047]

[0048] Among them, is the structural quality; the constraint conditions are specifically: the equivalent stress of the cruise structure shall not be greater than the allowable stress , the maximum deformation of the wing shall not exceed the set value of the deformation amount , and the torsional deformation amount of the wing under the yaw control state shall not exceed the set value of the torsional angle ;

[0049] The optimization mathematical model of the electromagnetic discipline is specifically as follows:

[0050]

[0051]

[0052] Among them, is the electromagnetic shielding rate; the constraint conditions are specifically: the electromagnetic shielding rate shall not be greater than the electromagnetic shielding rate of the initial configuration .

[0053] Furthermore, in steps 3 and 4, an improved rapid electromagnetic shielding evaluation method is adopted for electromagnetic target analysis, which specifically includes the following steps:

[0054] Step a2: In the finite element software, for the joined-wing aircraft with the antenna array layout to be calculated, ignoring the thickness of the patch antenna, divide the aircraft surface into a patch antenna area and other areas, generate triangular mesh elements according to the areas, and perform mesh encryption on the leading and trailing edges of the wing, and store the vertex coordinates of the mesh elements;

[0055] Step b2: Import the joined-wing aircraft in step a2 into the electromagnetic analysis software. In the electromagnetic analysis software, set the field source on the ground and set the receiving end on the joined-wing aircraft; generate a plane wave as the incident wave in the direction from the field source to the receiving end; the incident wave consists of several rays carrying the set electric field, and the number of rays is equal to the number of mesh elements in step a2; perform ray and mesh element intersection judgment. If the ray intersects with the mesh element, it is determined that the ray shoots towards the bright area;

[0056] Step c2: Perform field strength tracking on the rays shooting towards the bright area, and use the improved bounce ray method to calculate the field strength distribution of each patch antenna on the joined-wing aircraft;

[0057] Step d2: Calculate the electromagnetic shielding rate of the integrated joined-wing aircraft based on the field strength distribution and the field source electric field strength in all patch antenna regions .

[0058] Furthermore, in step b2, the process of determining whether the ray is directed towards the bright area is as follows: For each ray, traverse each grid cell on the surface of the integrated joined-wing aircraft model:

[0059] For each grid cell, first perform a single occlusion judgment:

[0060] Obtain the ray direction vector based on the field source position and the integrated joined-wing aircraft model position , and obtain the normal vector of the grid cell based on the vertex coordinates of the grid cell ;

[0061] If , it is determined that there is no intersection, and this grid cell is marked as the dark area, and this ray is directed towards the dark area; if , it is determined that there is an intersection, and this grid cell is marked as the bright area, and this ray is directed towards the bright area;

[0062] Then consider the mutual occlusion between grid cells and perform a multiple occlusion judgment:

[0063] For a grid cell marked as the bright area, draw a line segment from its center point to the field source and determine whether this line segment intersects with other grid cells. If this line segment intersects with other grid cells, then this grid cell is occluded by other grid cells, and this grid cell is re-marked as the dark area;

[0064] Finally, obtain all the rays directed towards the bright area.

[0065] Furthermore, step c2 specifically includes the following sub-steps:

[0066] Step c2.1: For the rays directed towards the bright area, assume that the rays propagate in a straight line in a homogeneous medium, and calculate the electric field strength at the intersection point where the rays propagate in a straight line through the following formula :

[0067]

[0068] Among them, is the differential of the grid cell at the field source, is the differential of the grid cell at the intersection point; the complex expression is used to describe the propagation characteristics of plane electromagnetic waves in space, is the imaginary unit, is the wave number, is the distance from the intersection point to the field source, is the wavelength of the ray; is the electric field strength of the field source;

[0069] Step c2.2: Determine the intersection type between the ray and the integrated joined-wing aircraft model. The intersection types include intersection in the plane and intersection at the edge. If the ray intersects with the grid cell in the plane, use the geometric optics method to simulate the ray reflection and calculate the direction vector of the reflected ray. and the electric field strength of the reflected ray. ; If the intersection is at the edge, use the improved uniform theory of diffraction to simulate the ray diffraction and calculate the direction vector of the diffracted ray. and the electric field strength of the diffracted ray. ;

[0070] Step c2.3: Trace each reflected and diffracted ray, and perform the intersection judgment between the ray and the grid cell. If the ray intersects with the grid cell, the ray is directed towards the bright area.

[0071] Step c2.4: Repeat steps c2.1 to c2.3 until the electric field carried by the traced ray is less than the set value or the ray is directed towards the dark area.

[0072] Step c2.5: For each patch antenna in the integrated joined-wing aircraft patch antenna area, obtain the field strength distribution of the patch antenna through the following process:

[0073] For each grid cell corresponding to the patch antenna, calculate the sum of the electric fields carried by the direct, reflected, and diffracted rays on each grid cell to obtain the total electric field of each grid cell on the patch antenna:

[0074] is the total electric field of the th grid cell in the patch antenna, so as to obtain the field strength distribution of the patch antenna; traverse each patch antenna to obtain the field strength distribution of each patch antenna area on the integrated joined-wing aircraft.

[0075] Furthermore, the method for determining the intersection type in step c2.2 is that if the angle between the normal of a bright area grid cell and the normal of the adjacent grid cell is less than 10 degrees, the intersection type is judged as intersection in the plane; if the normal angle is greater than or equal to 10 degrees, the intersection type is judged as intersection at the edge, and the bright area grid cell is the edge grid cell.

[0076] Furthermore, in step c2.2, when the intersection type is intersection at the edge, the process of using the improved uniform theory of diffraction to simulate the ray diffraction and calculate the direction vector of the diffracted ray and the electric field strength of the diffracted ray is specifically as follows:

[0077] First, solve for the ray direction vector after diffraction :

[0078] Obtain the vertex coordinates of the edge grid cell and the grid cells adjacent to it, and obtain the normal vector of the edge grid cell based on the vertex coordinates of the edge grid cell ; Obtain the edge vector based on the common edge of the edge grid cell and the grid cells adjacent to it ;

[0079] Let the incident direction vector of the ray be , and the angle between the incident direction vector and the normal vector be the incident angle ; The angle between the incident direction vector of the ray and the edge vector is denoted as ; Consider the diffraction phenomenon as the incident ray deflecting by an angle around the edge vector . Let the ray direction vector after diffraction be , and the angle between the ray direction vector after diffraction and the normal vector be the diffraction angle ;

[0080] Calculate the incident angle based on the incident direction vector of the ray and the normal vector :

[0081] Calculate the deflection angle using the formula , where is the diffraction coefficient;

[0082] Assume , is the triaxial component of the normal vector ; Obtain the transformation matrix for rotating around the edge vector through the following expression:

[0083]

[0084] Decompose into a vector parallel to and a vector perpendicular to , and the specific expressions are:

[0085]

[0086]

[0087] Then the ray direction vector after diffraction is:

[0088] Then solve for the ray electric field strength after diffraction :

[0089] According to the incident direction vector of the ray and the edge vector calculate :

[0090] According to the diffraction direction vector of the ray and the normal vector calculate the diffraction angle :

[0091] The ray electric field strength after diffraction is obtained by calculating the following formula :

[0092] In the formula, is the electric field strength when the ray propagates linearly to the intersection point, and at this time the intersection point is the diffraction point, is the diffraction coefficient tensor, is the defocusing distance, is the distance from the observation point to the diffraction point;

[0093] The diffraction coefficient tensor is calculated by the following formula :

[0094]

[0095] D s = - e -jπ / 4 2π 2πbsin β 0 [ F( b La( θ s - θ' i ) ) cos (( θ s - θ' i ) / 2) - F( bLa( θ s + θ' i ) ) cos (( θ s + θ' i ) / 2) ]

[0096] D h = - e -jπ / 4 2π 2πbsin β 0 [ F bLa θ s - θ' i cos θ s - θ' i 2 + F( bLa( θ s + θ' i ) ) cos (( θ s + θ' i ) / 2) ]

[0097] The functions in the above expressions , , and intermediate variables The explanation is as follows:

[0098] When : F x =[ πx -2x e j π 4 -2 x 2 e -j π 4 ] e j( π 4 +x)

[0099] When :

[0100]

[0101]

[0102] 。

[0103] Beneficial effects:

[0104] The present invention provides a multidisciplinary coupling optimization design method for an integrated joined-wing layout, mainly realizing the multidisciplinary coupling optimization of the joined-wing aerodynamics / structure / electromagnetics. Aiming at the problems of low efficiency, insufficient convergence, difficult gradient solution, and low accuracy in the existing technology, the present invention uses a multidisciplinary secondary optimization method to divide the design variables into global design variables and local design variables of each sub-discipline, greatly reducing the dimension of the optimization design and solving the problem of low efficiency; the design spaces of the system level and the subsystem level do not overlap, and there is no need to introduce coupling auxiliary design variables and state coupling design variables, so the convergence is good; when classifying the design variables, the classification basis is the coupling relationship of each discipline of the joined-wing, and there is no need to solve the gradient, so the calculation amount is small; in addition, for the integrated joined-wing layout, the structural target calculation adopts the finite element equivalent beam method based on the engineering beam theory. On the one hand, it combines the classical structural force transmission model of the engineering beam theory, and on the other hand, it adds a finite element beam model, which can be used to analyze the special stiffness configuration of the joined-wing; in addition, in the electromagnetic discipline optimization, an improved fast evaluation method for electromagnetic shielding is adopted, which is improved on the basis of the shooting and bouncing ray method. On the one hand, the diffraction coefficient is directly defined to quickly solve the diffraction angle; on the other hand, the three-order diffraction tensor coefficient in the traditional shooting and bouncing ray method is reduced in dimension to avoid magnetic field calculation, and the electric field of the diffracted ray can be quickly solved, which can provide an efficient and accurate evaluation method for the structural design and the multidisciplinary design of aerodynamics / structure / electromagnetics of the integrated joined-wing layout. Brief Description of the Drawings

[0105] Figure 1 is the flowchart of the implementation of the solution of the embodiment of the present invention;

[0106] Figure 2 This is the initial configuration of the joined wing in the embodiment of the present invention;

[0107] Figure 3 This is the three-dimensional model for electromagnetic analysis of the simplified joined wing in the embodiment of the present invention;

[0108] Figure 4 This is the detailed flowchart of the method for rapid equivalent analysis of the structural characteristics in the embodiment of the present invention;

[0109] Figure 5 This is the three-dimensional model for structural analysis of the simplified joined wing in the embodiment of the present invention;

[0110] Figure 6 This is the schematic diagram of the sectional division of the joined wing in the embodiment of the present invention;

[0111] Figure 7 This is the detailed flowchart of the method for analyzing the electromagnetic shielding rate in the embodiment of the present invention;

[0112] Figure 8 This is the grid schematic diagram of the electromagnetic shielding rate analysis in the embodiment of the present invention;

[0113] Figure 9 This is the comparison diagram of the aerodynamic / structural three-dimensional models of the joined wing before and after optimization in the embodiment of the present invention;

[0114] Figure 10 This is the comparison diagram of the electromagnetic three-dimensional models of the joined wing before and after optimization in the embodiment of the present invention;

[0115] Wherein: 1 - first patch antenna; 2 - second patch antenna; 3 - third patch antenna; 4 - fourth patch antenna; 5 - fuselage; 6 - connecting end plate; 7 - front wing; 8 - rear wing; 9 - first front wing beam; 10 - second front wing beam; 11 - front wing rib; 12 - front wing skin; 13 - first rear wing beam; 14 - second rear wing beam; 15 - rear wing rib; 16 - rear wing skin; 17 - single-sided connecting end plate; 18 - aerodynamic and structural model before optimization; 19 - aerodynamic and structural model after optimization; 20 - electromagnetic model before optimization; 21 - electromagnetic model after optimization; Specific Embodiments

[0116] In order to make the technical problems, technical solutions and beneficial effects solved by the present invention clearer and more understandable, and to enable those skilled in the art to better understand the solution of the present invention, the following further details and fully describe the present invention in conjunction with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0117] This embodiment takes a certain double-beam joined wing ultra-long endurance unmanned aerial vehicle as an example, and the initial configuration of this joined wing aircraft is as Figure 2As shown, this configuration mainly consists of a fuselage, a forward wing, a rear wing, connecting end plates, and a tail wing. Propellers for vertical takeoff and landing are arranged at both ends of the two connecting end plates, with a total of 4 propellers. For the sake of simplicity in the view, Figure 2 is not shown in the figure. Next, the multidisciplinary coupling optimization design method proposed by the present invention for the integrated joined-wing layout is adopted to conduct multidisciplinary optimization design of aerodynamics, structure, and electromagnetics for this joined-wing aircraft. To simplify the problem, a half-model design is carried out for this double-beam joined-wing ultra-long-endurance loitering unmanned platform. Taking the electromagnetic analysis model as an example, the simplified electromagnetic analysis model is as shown in Figure 3 : The simplified model consists of a fuselage, connecting end plates, a forward wing, and a rear wing. The antenna layout scheme is as follows: 4 patch antennas are evenly arranged on the side of the fuselage, and the thickness of the antenna is ignored.

[0118] As shown in Figure 1 , the flowchart of the implementation of the solution of the embodiment of the present invention specifically includes the following steps:

[0119] Step 1: Establish the initial configuration of the aircraft with an integrated joined-wing layout, and parameterize the configuration; determine the design objectives, design conditions, and design variables for the multidisciplinary optimization of the aircraft with an integrated joined-wing layout, as well as the design space range of each design variable;

[0120] The design objectives include a high lift-to-drag ratio, a low structural mass, and a low electromagnetic shielding rate; the design variables include aerodynamic discipline design variables, structural discipline design variables, and electromagnetic discipline design variables;

[0121] In this embodiment, the design conditions are as follows:

[0122] Cruise condition: angle of attack 0°, flight speed 40 m / s, flight altitude: H = 5 km;

[0123] Vertical takeoff and landing condition: including two states, a hover state and a yaw control state. Among them, in the hover state, the pulls of the four propellers are the same, that is, , where, is the pull of the left front propeller, is the pull of the right front propeller, is the pull of the left rear propeller, is the pull of the right rear propeller; in the yaw control state, the diagonal rotors increase to 1.5 times the lift, that is, , .

[0124] The design variables, design space, and parameterization method are specifically as follows:

[0125] The aerodynamic discipline design variables include the two-dimensional airfoil design parameters of the forward wing and the rear wing and the layout parameters of the three-dimensional shape of the joined wing;

[0126] The two-dimensional airfoil adopts the CST parameterization method, and the airfoil function is expressed as the product of a class function and a shape function, where the coefficients of the shape function are the design parameters of the two-dimensional airfoil;

[0127] In this embodiment, only the front airfoil is optimized; the front airfoil is selected as the airfoil to be optimized, and the CST parameterization method is adopted. The mathematical expression of the airfoil is:

[0128] In the formula, is the ordinate of the airfoil, is the abscissa of the airfoil, is the class function, is the shape function;

[0129]

[0130]

[0131]

[0132] In the formula, 、 are the category parameters for controlling the geometry. In this embodiment, 、 are taken; the coefficients of the shape function are the weight coefficients, is the Bernstein polynomial, and n is its order. In this embodiment, n = 6 is taken;

[0133] are the design parameters of the two-dimensional airfoil. In this embodiment, the two-dimensional airfoil is parameterized for the upper and lower surfaces, and the weight coefficients - represent the design parameters of the upper surface, and the weight coefficients - represent the design parameters of the lower surface;

[0134] The following table gives the initial values and upper and lower value ranges (design space ranges) of the design parameters of the two-dimensional airfoil in the design variables of the aerodynamic discipline:

[0135]

[0136] The three-dimensional external layout parameters include:

[0137] The front wingspan 、the front wing root chord length 、the front wing root taper ratio 、the front wing tip twist angle 、the front wing sweep angle 、the front wing dihedral angle ; the rear wingspan 、the rear wing root chord length , the trailing wing root ratio , the trailing wing tip twist angle , the trailing wing forward sweep angle , the trailing wing dihedral angle ; the leading edge position of the forward wing , the leading edge position of the trailing wing ; the endplate height ratio ;

[0138] Among them, is the height difference between the forward wing and the trailing wing at the wing root, is the height difference between the forward wing and the trailing wing at the endplate;

[0139] The following table gives the initial values and upper and lower value ranges (design space ranges) of the layout parameters of the three-dimensional shape in the design variables of the aerodynamic discipline:

[0140]

[0141] The described structural discipline design variables include structural dimension parameters and structural layout parameters;

[0142] Among them, the structural dimension parameters include:

[0143] The thickness of the web of the first forward wing beam: , the thickness of the web of the second forward wing beam: ; the area of the upper flange of the first forward wing beam: , the area of the lower flange of the first forward wing beam: ; the area of the upper flange of the second forward wing beam: , the area of the lower flange of the second forward wing beam: ; the thickness of the forward wing rib: ; the thickness of the forward wing skin: ; the thickness of the web of the first trailing wing beam: , the thickness of the web of the second trailing wing beam: ; the area of the upper flange of the first trailing wing beam: , the area of the lower flange of the first trailing wing beam: ; the area of the upper flange of the second trailing wing beam: , the area of the lower flange of the second trailing wing beam: ; the thickness of the trailing wing rib: ; the thickness of the trailing wing skin: ; the structural layout parameters include: the position of the first forward wing beam: , the position of the second forward wing beam: ; the forward wing rib spacing: ; the position of the first trailing wing beam: , the position of the second trailing wing beam: ; the trailing wing rib spacing: ;

[0144] The following table gives the initial values and upper and lower value ranges (design space ranges) of the structural dimension parameters among the design variables of the structural discipline:

[0145]

[0146] The following table gives the initial values and upper and lower value ranges (design space ranges) of the structural layout parameters among the design variables of the structural discipline:

[0147]

[0148] The design variables of the electromagnetic discipline include the antenna spacing , the number of this parameter is related to the number of antennas. Assuming that m antennas are arranged on the side of the fuselage, then ; among them, is the distance between the leading edge of the forewing and the first antenna, is the distance between the first antenna and the second antenna, and so on, is the distance between the (m - 1)-th antenna and the second antenna.

[0149] In this embodiment, a total of 4 antennas are arranged on the side of the fuselage, and the size specification of the antenna is 200 250 mm. The following table gives the initial values and upper and lower value ranges (design space ranges) of the design variables of the electromagnetic discipline:

[0150]

[0151] In this embodiment, based on the above design variables, the specific way to parameterize the configuration is:

[0152] Use secondary development of CATIA software to parameterize the three-dimensional joined-wing configuration, including parametric modeling of the aerodynamic / structural / electromagnetic disciplines; specifically: when using CATIA to establish the digital model of the three-dimensional joined-wing configuration, take the leading-edge position coordinates of the wing root and wing tip, the ratio of the wing root and wing tip, and the rotation angle to change the wing tip profile as the control parameters of the aerodynamic discipline design variables; take the thick surface size of the skin, the sketch size of the beam section, the coordinates of the beam section, the thick surface size of the rib, and the translation matrix size of the rib as the control parameters of the structural discipline design variables; take the coordinates of the leading edge of the antenna as the control parameters of the electromagnetic discipline design variables; and record the script file, and control the change of the three-dimensional joined-wing configuration by modifying the above control parameters in the script file.

[0153] Among them, in the aerodynamic discipline, the span lengths of the front and rear wings, the sweep angles of the front wing, the dihedral angles of the front wing, the sweep angles of the rear wing, the anhedral angles of the rear wing, the leading-edge positions of the front and rear wings, and the end-plate height ratio are controlled by changing the leading-edge position coordinates of the wing root and wing tip. The root chord lengths and root-to-tip ratios of the front and rear wings are controlled by changing the ratios of the wing root and wing tip. The tip twist angles of the front and rear wings are controlled by changing the rotation angles of the tip profiles.

[0154] In the structural discipline, the skin thickness of the front and rear wings is controlled by changing the thick-curved surface dimensions of the skin. The web thickness of the beams and the areas of the upper and lower flanges of the front and rear wings are controlled by changing the sketch dimensions of the beam profiles. The positions of the front and rear wing beams are controlled by changing the coordinates of the beam profiles. The thicknesses of the ribs of the front and rear wings are controlled by changing the thick-curved surface dimensions of the ribs. The rib spacings of the front and rear wings are controlled by changing the translation matrix dimensions of the ribs.

[0155] In the electromagnetic discipline, the antenna spacing is controlled by changing the coordinates of the antenna leading edge.

[0156] Step 2: Divide the design variables determined in Step 1 into global design variables XS and local design variables XL;

[0157] The global design variables XS include all the design parameters of the two-dimensional airfoil in the design variables of the aerodynamic discipline and the following parameters in the three-dimensional outer shape layout parameters of the connected wing:

[0158] Span length of the front wing 、Root chord length of the front wing 、Root-to-tip ratio of the front wing 、Sweep angle of the front wing ; Span length of the rear wing 、Root chord length of the rear wing 、Root-to-tip ratio of the rear wing 、Sweep angle of the rear wing ; Leading-edge position of the front wing ,Leading-edge position of the rear wing ;End-plate height ratio ;

[0159] The local design variables XL include the local design variables of the aerodynamic discipline 、The local design variables of the structural discipline 、The local design variables of the electromagnetic discipline ;

[0160] Among them, the local design variables of the aerodynamic discipline include the tip twist angle of the front wing 、Dihedral angle of the front wing ;Tip twist angle of the rear wing 、Anhedral angle of the rear wing ;

[0161] The local design variables of the structural discipline Include all structural dimension parameters and structural layout parameters;

[0162] Local design variables in the electromagnetic discipline Include all antenna spacings .

[0163] Step 3: Perform system-level optimization;

[0164] Specifically: Taking high lift-to-drag ratio, small structural mass, and small electromagnetic shielding ratio as system-level design objectives, using global design variables as system-level optimization design variables, and fixing local design variables as constants; Using the CFD method for aerodynamic target analysis, using the equivalent beam method based on engineering beam theory for structural target analysis, and using an improved fast electromagnetic shielding evaluation method for electromagnetic target analysis; Using the genetic algorithm for optimization to obtain the global optimal design result;

[0165] In this embodiment, the mathematical model of system-level optimization is specifically as follows:

[0166]

[0167]

[0168]

[0169]

[0170] Among them, is the weight coefficient of the lift-to-drag ratio optimization target, is the weight coefficient of the structural mass optimization target, is the weight coefficient of the electromagnetic shielding ratio target, and all are taken as 1 in this embodiment; are the lift-to-drag ratio, structural mass, and electromagnetic shielding ratio respectively, are the lift-to-drag ratio, structural mass, and electromagnetic shielding ratio of the initial configuration respectively; The constraint conditions include aerodynamic discipline constraints, structural discipline constraints, and electromagnetic discipline constraints, specifically:

[0171] Aerodynamic discipline constraint: The cruise lift coefficient is not lower than the lift set value , and in this embodiment, it is taken as , and the change in pitch moment coefficient is not greater than the moment change set value , and in this embodiment, it is taken as ;

[0172] Structural discipline constraint: The cruise structural equivalent stress is not greater than the allowable stress , and the maximum wing deformation does not exceed the deformation set value , in this embodiment, take , unit: meter, the wing twist deformation amount under the yaw control state shall not exceed the set value of the twist angle , in this embodiment, take ;

[0173] Electromagnetic discipline constraint: the electromagnetic shielding rate shall not be greater than the electromagnetic shielding rate of the initial configuration .

[0174] Step 4: Conduct subsystem-level optimization respectively, including aerodynamic discipline optimization, structural discipline optimization, and electromagnetic discipline optimization;

[0175] Specifically: The aerodynamic discipline optimization takes the high lift-to-drag ratio as the aerodynamic discipline optimization design goal, takes the local design variables in the aerodynamic discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the structural discipline design variables, and the local design variables in the electromagnetic discipline design variables as constants; uses the CFD method to conduct aerodynamic target analysis, and uses the genetic algorithm to conduct optimization to obtain the local optimal design result of the aerodynamic discipline of the joined-wing layout;

[0176] In this embodiment, the mathematical model of the aerodynamic discipline optimization is specifically as follows:

[0177]

[0178]

[0179] Among them, is the lift-to-drag ratio; the constraint conditions are specifically: the cruise lift coefficient shall not be lower than the set value of the lift , in this embodiment, take , the change amount of the pitching moment coefficient shall not be greater than the set value of the moment change , in this embodiment, take ;

[0180] The structural discipline optimization takes the low structural mass as the structural discipline optimization design goal, takes the local design variables in the structural discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic discipline design variables after aerodynamic discipline optimization, and the local design variables in the electromagnetic discipline design variables as constants; combines the aerodynamic loads obtained by using the CFD method to conduct aerodynamic target analysis in the aerodynamic discipline optimization process, uses the equivalent beam method based on the engineering beam theory to conduct structural target analysis, and uses the genetic algorithm to conduct optimization to obtain the local optimal design result of the structural discipline of the joined-wing layout;

[0181] In this embodiment, the optimization mathematical model of the structural discipline is as follows:

[0182]

[0183]

[0184] Among them, is the structural mass; the constraint conditions are specifically: the equivalent stress of the cruise structure is not greater than the allowable stress , the maximum deformation of the wing does not exceed the set value of the deformation amount . In this embodiment, is taken. Under the yaw control state, the torsional deformation of the wing does not exceed the set value of the torsional angle . In this embodiment, ;

[0185] The optimization of the electromagnetic discipline takes the low electromagnetic shielding rate as the design goal of the electromagnetic discipline, uses the local design variables in the electromagnetic discipline design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic discipline design variables after aerodynamic discipline optimization, and the local design variables in the structural discipline design variables after structural discipline optimization as constants; adopts an improved rapid evaluation method for electromagnetic shielding to conduct electromagnetic target analysis, and uses the genetic algorithm for optimization to obtain the local optimal design result of the electromagnetic discipline of the joined-wing layout;

[0186] In this embodiment, the optimization mathematical model of the electromagnetic discipline is as follows:

[0187]

[0188]

[0189] Among them, is the electromagnetic shielding rate; the constraint conditions are specifically: the electromagnetic shielding rate is not greater than the electromagnetic shielding rate of the initial configuration .

[0190] In this embodiment, as Figure 4 shown, in steps 3 and 4, the equivalent beam method based on the engineering beam theory is used for structural target analysis to obtain the equivalent stress of the cruise structure , the maximum deformation of the wing , the torsional deformation of the wing under the yaw control state and the specific steps of the overall structural mass are as follows. The equivalent beam method based on the engineering beam theory can refer to the publication number Chinese Patent Application "Fast Equivalent Analysis Method for Structural Characteristics Oriented to Integrated Joined-Wing Configuration Design":

[0191] First, establish a global coordinate system , and the origin of the global coordinate system is established at the design center of gravity of the integrated joined-wing layout aircraft, the axis points backward along the chord line of the longitudinal symmetric section of the aircraft, the axis is located in the section and perpendicular to the axis upward, the axis is perpendicular to the section and inward to the aircraft, satisfying the right-hand system;

[0192] Step a1: Obtain the design variables of the joined-wing configuration to be analyzed. Regard the joined-wing to be analyzed as a spatial rigid frame structure formed by the combination of the front wing, the rear wing, and the connecting end plate, and simplify it into a spatial six-degree-of-freedom rigid frame structure model; the "rigid frame" here refers to a rigid frame structure, rather than a flexible structure;

[0193] In this embodiment, regard the joined-wing shown in Figure 2 as a spatial rigid frame structure formed by the combination of the front wing, the rear wing, and the connecting end plate, and simplify it into a spatial six-degree-of-freedom rigid frame structure model; the simplified three-dimensional structural analysis model is as shown in Figure 5 .

[0194] Step b1: In the spatial six-degree-of-freedom rigid frame structure model, divide the front wing, the rear wing, and the connecting end plate into several sections at equal intervals according to the set interval, and obtain the leading-edge point coordinates of each section in the global coordinate system; among them, the sections of the front wing and the rear wing are parallel to the longitudinal symmetric plane of the aircraft where the joined-wing is located, and the section of the connecting end plate is perpendicular to the longitudinal symmetric plane of the aircraft where the joined-wing is located;

[0195] In this embodiment, divide the spatial six-degree-of-freedom rigid frame structure model according to the following two principles, and the division schematic diagram is as shown in Figure 6 :

[0196] 1) Enough sections need to be divided for the front wing and the rear wing to accurately describe the global stiffness. The front wing takes 9 sections at an interval of 0.3 m, and the rear wing takes 5 sections at an interval of 0.275 m;

[0197] 2) The stiffness of the connecting end plate contributes little to the global stiffness, and 4 sections are taken at equal intervals;

[0198] The method for obtaining the leading-edge point coordinates of each section is to first obtain the leading-edge point coordinates of the two end sections and the connecting section of the front wing, the rear wing, and the connecting end plate as follows:

[0199] #1 (-0.8505, 0, 0.1610); #5 (-0.4058, -1.1, 0.1347); #9 (0.1197, -2.4, 0.1035); #10 (0.5633, 0, 0.3579); #14 (-0.0569, -1.1, 0.2173); #15 (-1.158, -1.1, 0.263); #18 (0.842, -1.1, 0.263);

[0200] Then, the coordinates of the leading edge points of the remaining cross-sections are calculated through geometric relationships.

[0201] Step c1: Based on the structural discipline design variables and the coordinates of the leading edge points of each cross-section, calculate the centroid coordinates of each cross-section in the global coordinate system Take the centroid of each cross-section as the origin, establish the centroid coordinate system of each cross-section, and the three-axis directions of the centroid coordinate system are correspondingly parallel to the three-axis directions of the global coordinate system;

[0202] Taking the front wing root cross-section, i.e., cross-section #1, as an example, calculate its centroid coordinates and obtain its centroid coordinate system , which specifically includes the following sub-steps:

[0203] Step c1.1: Calculate the total structural area of the cross-section ;

[0204] In the formula, is the reduction coefficient or conversion coefficient, representing the ratio of the Young's modulus of the material used to the set reference material. In this embodiment, the materials used are all aluminum alloys, so the reduction coefficient ; is the skin thickness of the corresponding area of the cross-section. For example, in this embodiment, this cross-section is the front wing root cross-section, so is the front wing skin thickness , is the skin arc length micro-element; is the beam web thickness of the corresponding area of the cross-section. In this embodiment, is the front wing beam web thickness , is the beam web arc length micro-element; is the area of the k-th type of wing beam flange in the corresponding area of the cross-section. In this embodiment, there are 4 types of wing beam flanges in the corresponding area of the front wing root cross-section, namely the area of the first front wing beam upper flange , the area of the first front wing beam lower flange , the area of the second front wing beam upper flange and the area of the second front wing beam lower flange .

[0205] Step c1.2: Calculate the total static moment of the section about the X-axis in the global coordinate system and the total static moment about the Z-axis

[0206]

[0207]

[0208] wherein, is the Z-direction coordinate of the skin arc length microelement in the corresponding area of the section in the global coordinate system, is the Z-direction coordinate of the beam web arc length microelement in the corresponding area of the section in the global coordinate system, is the Z-direction coordinate of the k-th type of wing beam flange in the corresponding area of the section in the global coordinate system; is the X-direction coordinate of the skin arc length microelement in the corresponding area of the section in the global coordinate system, is the X-direction coordinate of the beam web arc length microelement in the corresponding area of the section in the global coordinate system, is the X-direction coordinate of the k-th type of wing beam flange in the corresponding area of the section in the global coordinate system.

[0209] Step c1.3: Calculate the centroid coordinates of the section ;

[0210]

[0211] Move the origin of the global coordinate system to the centroid coordinates , to obtain the centroid coordinate system corresponding to the section;

[0212] Step d1: Calculate the moment of inertia about the X-axis and the moment of inertia about the Z-axis of each section in its respective centroid coordinate system in step b1;

[0213] Specifically, since the angle between the centroid coordinate system and the section's principal inertia axis coordinate system is a small quantity, the centroid coordinate system is regarded as the principal inertia axis coordinate system. Thus, the moments of inertia and of the sections divided in step b1 about the principal inertia axis, which is also the centroid axis, are obtained through the following calculation formulas:

[0214]

[0215] ;

[0216] Step e1: Calculate the rigid centroid coordinates of each section in its respective centroid coordinate system in step b1. Taking the rigid centroid of each section as the origin, establish the rigid centroid coordinate system of each section, and the three-axis directions of the rigid centroid coordinate system are correspondingly parallel to the three-axis directions of the centroid coordinate system;

[0217] Taking section #1 as an example specifically, the description is as follows:

[0218] Using the engineering beam theory, the wing is regarded as a beam supported on the fuselage. According to the moment of inertia in the centroid coordinate system of the section obtained in step d1 and , the rigid centroid coordinates of the section are obtained through the following calculation formulas:

[0219]

[0220]

[0221] Among them, is the total static moment of the section on the X-axis of the centroid coordinate system, is the total static moment of the section on the Z-axis of the centroid coordinate system, is the perpendicular distance from the resultant shear flow force of the skin arc length microelement in the corresponding area of the section to the set moment reference point, is the perpendicular distance from the resultant shear flow force of the beam web arc length microelement in the corresponding area of the section to the set moment reference point is twice the area enclosed by the single closed chamber closed curve corresponding to the skin arc length microelement in the corresponding area of the section, is twice the area enclosed by the single closed chamber closed curve corresponding to the beam web arc length microelement in the corresponding area of the section, is the shear modulus of the set reference material;

[0222] Move the origin of the centroid coordinate system to the rigid centroid coordinates to obtain the corresponding rigid centroid coordinate system of the section;

[0223] Step f1: Calculate the moment of inertia about the X-axis , moment of inertia about the Z-axis and torsional stiffness of each section in its respective rigid centroid coordinate system in step b1;

[0224] Specifically, taking section #1 as an example, the description is as follows:

[0225] According to the rigid centroid coordinates of the section obtained in step e1, the moment of inertia about the X-axis and moment of inertia about the Z-axis of the section in the rigid centroid coordinate system are calculated using the following formula:

[0226]

[0227]

[0228] According to the formula

[0229] Obtain the torsional stiffness per unit length of the equivalent beam ;

[0230] Step g1: Based on the X-axis moment of inertia , Z-axis moment of inertia and torsional stiffness of each section in its own center-of-gravity coordinate system obtained in step f1, establish the element global stiffness matrix of each element and assemble to obtain the structural global stiffness matrix of the spatial six-degree-of-freedom rigid frame structure ; Specifically, it includes the following sub-steps:

[0231] Step g1.1: Consider each section as a node, and the position of the node is the center-of-gravity position of the section. Two adjacent nodes and node form an element; the local coordinate system of the element takes node as the origin, and the coordinate axes are respectively axis, axis, axis. Among them, axis is the connection line between node and node , axis coincides with the axis of the center-of-gravity coordinate system of the section where node is located, and axis forms a right-hand coordinate system with axis and axis;

[0232] The simplified configuration of the connecting wing is a spatial six-degree-of-freedom rigid frame. Therefore, each node has 6 degrees of freedom, each element has 12 degrees of freedom, and each element is considered as a beam element with a constant cross-section. Based on the moments of inertia and torsional stiffness of each section in its own center-of-gravity coordinate system obtained in step f1, obtain the local stiffness matrix of the element formed by node and node :

[0233]

[0234] In the formula, is Young's modulus, is the shear modulus, is the total cross-sectional area of the section, and are the moments of inertia of the center-of-gravity coordinate system of the section, , is the element length, i.e., the distance between the two sections; the 1-6 columns of the matrix respectively correspond to the displacements of node in the direction displacement ( ), Directional displacement ( ), Directional displacement ( ); Directional rotation angle ( ); Directional rotation angle ( ); Directional rotation angle ( ); Similarly, the 6 degrees of freedom of the nodes corresponding to columns 7 - 12 of the matrix ;

[0235] Step g1.2: According to the formula

[0236] Convert the element local stiffness matrix in the local coordinate system of the element to the element global stiffness matrix in the global coordinate system, where

[0237]

[0238]

[0239] where is the coordinate transformation matrix, is the Euler transformation matrix, is the element global stiffness matrix, is the three - axis rotation Euler angles from the local coordinate system of the element to the global coordinate system;

[0240] Step g1.3: Assemble the element global stiffness matrix of each element to obtain the structural global stiffness matrix ;

[0241] Specifically: Assemble the element global stiffness matrix of each element into the structural global stiffness matrix in sequence according to the node numbers of the elements. The node numbers determine their positions in the structural global stiffness matrix;

[0242] In this embodiment, combined with the structural material properties of the skin, wing ribs, and wing beams, obtain the structural global stiffness matrix The size of the matrix is orders, so it will not be shown;

[0243] Step h1: According to the structural global stiffness matrix obtained in step g1, establish and solve the structural equation to obtain the cruise structural equivalent stress of the connected - wing configuration, wing deformation , wing torsional deformation amount under yaw control state, and the overall structural mass, specifically including the following sub - steps:

[0244] Step h1.1: Establish the structural equation, specifically:

[0245] 1) Load processing:

[0246] Based on the aerodynamic loads obtained from the aerodynamic target analysis using the CFD method during the optimization process of the aerodynamic discipline, obtain the loads of the nodes and express the loads of the nodes as: In the formula,

[0247] are the force loads in three directions of the node in the global coordinate system, and are the moment loads in three directions of the node in the global coordinate system;

[0248] Assemble the load vector and express it as:

[0249] For a configuration with the number of nodes , its load vector is a column vector with a size of , which is composed of the loads of each node;

[0250] After assembling to obtain the load vector , combined with the overall stiffness matrix obtained in step g1, establish the structural equation , where is the displacement vector matrix to be solved, including displacements and rotations in three directions, namely direction displacement, direction displacement, direction displacement; direction rotation, direction rotation, direction rotation;

[0251] 2) Boundary condition processing:

[0252] Limit the positions of the corresponding nodes at the front wing root and the rear wing root in the displacement vector to 0; and limit the displacement vector of the corresponding nodes connecting the end plates to 0;

[0253] Step h1.2: Solve the structural equation, specifically:

[0254] 1) From the structural equation we can obtain ​​, the deflection and rotation angle of each node are obtained; the deflection and rotation angle of each node include the following six parameters:

[0255] Directional displacement: , Directional displacement: , Directional displacement: ; Directional rotation angle: , Directional rotation angle: , Directional rotation angle: ;

[0256] 2) According to the deflection and rotation angle of each obtained node, interpolate the shape functions of the deflection and rotation angle for each element, obtain the deflection and rotation angle of any point within the element, and further obtain the wing deformation of the joined-wing configuration and the wing torsional deformation under yaw control , specifically as follows:

[0257] X-direction displacement shape function:

[0258] In the formula is the displacement of the element in the direction along the X-axis, is the undetermined coefficient, obtained by interpolation, Y-direction displacement shape function:

[0259] In the formula is the displacement of the element in the direction along the Y-axis, is the undetermined coefficient, obtained by interpolation, Z-direction displacement shape function:

[0260] In the formula is the displacement of the element in the direction along the Z-axis, is the undetermined coefficient, obtained by interpolation, rotation angle shape function about the X-axis:

[0261] In the formula is the rotation angle of the element about the X-axis in the position, is the undetermined coefficient, obtained by interpolation, rotation angle shape function about the Y-axis:

[0262] In the formula is the rotation angle of the element about the Y-axis in the position, is a coefficient to be determined, obtained by interpolation. The shape function of the rotation angle about the Z-axis:

[0263] where is the rotation angle of the element about the Z-axis at the position ; is a coefficient to be determined, obtained by interpolation. The wing deformation of the joined-wing configuration is obtained ; the wing twist deformation under the yaw control state ;

[0264] 3) According to the obtained displacement shape functions in the X, Y, and Z directions, calculate the equivalent stress of the element structure of the joined-wing configuration , specifically as follows:

[0265] According to the displacement shape functions, obtain the linear strains in the X-axis, Y-axis, and Z-axis directions:

[0266]

[0267] where is the linear strain in the X-axis direction, is the linear strain in the Y-axis direction, is the linear strain in the Z-axis direction;

[0268] According to the displacement shape functions, obtain the shear strains in the XY plane, YZ plane, and ZX plane:

[0269] , ,

[0270] where is the shear strain in the XY plane, is the shear strain in the YZ plane, is the shear strain in the ZX plane;

[0271] Based on the obtained linear strains, calculate the normal stress as follows:

[0272]

[0273]

[0274]

[0275] where is the Poisson's ratio of the material;

[0276] Based on the obtained shear strains, calculate the shear stress as follows:

[0277]

[0278] Based on the obtained normal stress and shear stress, the stress state tensor is obtained

[0279]

[0280] Solve its characteristic equation: ; Obtain the eigenvalues , that is, the three principal stresses of the element 、 、 , is the identity matrix;

[0281] Based on the obtained three principal stresses 、 、 , the equivalent stress is calculated by using the fourth strength theory:

[0282] σ= 1 2 [ σ 1 - σ 2 2 + σ 2 - σ 3 2 + σ 3 - σ 1 2 ] ;

[0283] 4) Calculate the structural mass of the entire joined-wing configuration according to the parameters of each section , and the specific calculation formula is as follows:

[0284] m = ∑ i=1 n [Δ x i * A i+1 + A i 2 ] ρ mean

[0285] In the formula, 、 are the total structural areas of the sections where the nodes are located respectively, is the distance between the nodes , n is the number of nodes, is the average density of the section.

[0286] In this embodiment, as Figure 7 shown, in steps 3 and 4, an improved electromagnetic occlusion fast evaluation method is used for electromagnetic target analysis to obtain the electromagnetic occlusion rate of the joined-wing configuration, and the specific steps are as follows:

[0287] Step a2: In the finite element software, for the integrated joined-wing aircraft with the antenna array layout scheme to be calculated, ignoring the thickness of the patch antenna, divide the aircraft surface into a patch antenna area and other areas, generate triangular mesh elements according to the areas, and perform mesh encryption on the leading and trailing edges of the wing, and store the vertex coordinates of the mesh elements;

[0288] In this embodiment, the size of the patch antenna on the side of the fuselage is , and the operating frequency is , the corresponding wavelength size is ; The size of the triangular grid unit is selected according to the wavelength of the incident wave. Specifically, the grid unit size is taken as one-tenth of the wavelength, that is, 0.05 m; Due to the large curvature of the leading and trailing edges of the wing, the grid is appropriately refined to better fit the shape. There are a total of 7,134 grid units. The divided grid is as shown in Figure 8 ;

[0289] Step b2: Import the joined-wing aircraft in step a2 into the electromagnetic analysis software. In the electromagnetic analysis software, set the field source on the ground and set the receiving end on the joined-wing aircraft; Generate a plane wave as the incident wave in the direction from the field source to the receiving end; The incident wave is composed of several rays carrying the set electric field, and the number of rays is equal to the number of grid units in step a2; Perform ray-grid unit intersection judgment. If the ray intersects with the grid unit, it is determined that the ray is directed towards the bright area;

[0290] The process of determining whether the ray is directed towards the bright area is specifically as follows. For each ray, traverse each grid unit on the surface of the joined-wing aircraft;

[0291] For each grid unit, first perform single occlusion judgment:

[0292] Obtain the ray direction vector according to the field source position and the position of the joined-wing aircraft , and obtain the normal vector of the grid unit according to the vertex coordinates of the grid unit ;

[0293] In this embodiment, in the global coordinate system , assuming that the joined-wing aircraft is at an altitude of 1000 m and the field source is 2000 m to the left of the aircraft, then the receiving end coordinates , the field source coordinates , the ray direction vector After normalization, it is , and the intensity of the electric field initially carried by the ray (field source electric field intensity) is set to .

[0294] If , it is determined that there is no intersection, and this grid unit is marked as the dark area, and this ray is directed towards the dark area; If , it is determined that there is an intersection, and this grid unit is marked as the bright area, and this ray is directed towards the bright area;

[0295] Then consider the mutual occlusion between grid units and perform multiple occlusion judgment:

[0296] For a grid cell marked as a bright area, draw a line segment from the center point to the field source, and determine whether this line segment intersects with other grid cells. If the line segment intersects with other grid cells, then this grid cell is blocked by other grid cells, and this grid cell is re-marked as a dark area;

[0297] Finally, all the rays directed towards the bright area are obtained.

[0298] Step c2: Perform field strength tracking on the rays directed towards the bright area, and use the improved bounce ray method to calculate the field strength distribution of each patch antenna on the integrated joined-wing aircraft; specifically, it includes the following sub-steps:

[0299] Step c2.1: For the rays directed towards the bright area, assume that the rays propagate in a straight line in a homogeneous medium, and calculate the electric field strength at the intersection point where the rays propagate in a straight line through the following formula :

[0300]

[0301] Among them, is the differential of the grid cell at the field source, is the differential of the grid cell at the intersection point; the complex expression is used to describe the propagation characteristics of plane electromagnetic waves in space, is the imaginary unit, is the wave number, is the distance from the intersection point to the field source, is the wavelength of the ray;

[0302] Step c2.2: Determine the intersection type of the ray with the integrated joined-wing aircraft. The intersection types include intersection in the plane and intersection at the edge; if the ray intersects with the grid cell in the plane, use the geometric optics method to simulate the ray reflection, and calculate the direction vector of the reflected ray and the electric field strength of the reflected ray ; if the intersection is at the edge, use the improved uniform theory of diffraction to simulate the ray diffraction, and calculate the direction vector of the diffracted ray and the electric field strength of the diffracted ray ;

[0303] In this embodiment, the method for determining the intersection type is that if the included angle between the normal of a bright area grid cell and the normal of an adjacent grid cell is less than 10 degrees, it is determined as an intersection in the plane; if the included angle of the normal is greater than or equal to 10 degrees, it is determined as an intersection at the edge, and this bright area grid cell is an edge grid cell.

[0304] ① When the intersection is in the plane, use the geometric optics method (GO) to simulate the reflection of the ray, and calculate the direction vector of the reflected ray and the electric field intensity of the reflected ray The specific process is as follows:

[0305] First, solve the direction vector of the reflected ray :

[0306] Obtain the normal vector of the grid cell according to the vertex coordinates of the grid cell ;

[0307] Let the incident direction vector of the ray be , and the angle between the incident direction vector of the ray and the normal vector of the grid cell is the incident angle ; Let the direction vector of the reflected ray be , and the angle between the direction vector of the ray after reflection and the normal vector of the grid cell is the reflection angle ; According to the law of reflection, the incident angle is equal to the reflection angle ; In the plane formed by the incident direction vector and the normal vector , decompose into a vector parallel to the normal vector and a vector perpendicular to the normal vector ;

[0308] where is the projection vector of on

[0309] ;

[0310]

[0311] Because , so is symmetric to about the normal vector , that is

[0312] Obtain The expression of

[0313] Then solve the electric field intensity of the reflected ray through the following formula :

[0314] where is the electric field intensity when the ray propagates linearly to the intersection point obtained in step 31, and the intersection point is the reflection point at this time represents the reflection coefficient, is the differential of the grid cell at the incident point, is the differential of the grid cell after reflection, is the distance from the observation point to the incident point after reflection.

[0315] ② When the intersection is at the edge, use the improved uniform theory of diffraction to simulate the diffraction of the ray and calculate the direction vector of the ray after diffraction and the electric field strength of the ray after diffraction The specific process is as follows:

[0316] First, solve the direction vector of the ray after diffraction :

[0317] Obtain the vertex coordinates of the edge grid cell and the grid cells adjacent to it, and obtain the normal vector of the edge grid cell according to the vertex coordinates of the edge grid cell ; Obtain the edge vector according to the common edge of the edge grid cell and the grid cells adjacent to it ;

[0318] The main shielding components of the joined-wing aircraft are the connecting end plates and the front and rear wings. Simplify the diffraction problem to a flat-plate diffraction problem. Let the incident direction vector of the ray be , take the normal vector of the flat plate as , take the edge vector of the flat plate as , the angle between the incident direction vector of the ray and the normal vector is the incident angle ; The angle between the incident direction vector of the ray and the edge vector is denoted as ; Consider the diffraction phenomenon as the incident ray deflecting by an angle around the edge vector . Let the direction vector of the ray after diffraction be , and the angle between the direction vector of the ray after diffraction and the normal vector is the diffraction angle ;

[0319] According to the incident direction vector of the ray and the normal vector calculate the incident angle :

[0320] Define the diffraction coefficient , calculate the deflection angle :

[0321] Assume , is the three - axis component of the normal vector; The transformation matrix rotating around the edge vector is obtained by the following expression .

[0322]

[0323] Decompose into a vector parallel to and a vector perpendicular to . The specific expression is:

[0324]

[0325]

[0326] Then the diffracted ray direction vector is:

[0327] Then solve the diffracted ray electric field strength :

[0328] Calculate according to the incident direction vector of the ray and the edge vector :

[0329] Calculate the diffraction angle according to the diffracted direction vector of the ray and the normal vector :

[0330] According to Keller's cone, the diffracted ray electric field strength is obtained by the following formula

[0331] In the formula, is the electric field strength when the ray propagates linearly to the intersection point, and the intersection point is the diffraction point at this time. is the diffraction coefficient tensor, is the defocusing distance, is the distance from the observation point to the diffraction point; As mentioned before, the main shielding components of the integrated joined - wing aircraft are the connecting end plates and the front and rear wings. The diffraction problem is simplified to a flat - plate diffraction problem. For the flat - plate diffraction problem, there is , ;

[0332] For the diffraction coefficient tensor in the conventional bounce - ray method Perform dimensionality reduction. Without considering the spatial dimensional characteristics of the electromagnetic field, only consider the solution of the electric field, and calculate the diffraction coefficient tensor through the following formula :

[0333]

[0334] D s = - e -jπ / 4 2π 2πbsin β 0 [ F( b La( θ s - θ' i ) ) cos (( θ s - θ' i ) / 2) - F( bLa( θ s + θ' i ) ) cos (( θ s + θ' i ) / 2) ]

[0335] D h = - e -jπ / 4 2π 2πbsin β 0 [ F bLa θ s - θ' i cos θ s - θ' i 2 + F( bLa( θ s + θ' i ) ) cos (( θ s + θ' i ) / 2) ]

[0336] The functions and intermediate variables in the above expressions are explained as follows:

[0337]

[0338] When : F x =[ πx -2x e j π 4 -2 x 2 e -j π 4 ] e j( π 4 +x)

[0339] When :

[0340]

[0341]

[0342] 。

[0343] Step c2.3: Trace each ray after reflection and diffraction, and determine whether the ray intersects with the grid cell. If the ray intersects with the grid cell, then the ray is directed towards the bright area;

[0344] Step c2.4: Repeat Step c2.1 to Step c2.3 until the electric field carried by the traced ray is less than the set value or until it shoots into the dark area;

[0345] In this embodiment, the set value , that is, when the electric field strength of the ray after reflection or diffraction is less than , it is considered that the ray shoots into the dark area and no longer contributes to the superposition of the electric field strength of the patch antenna.

[0346] Step c2.5: For each patch antenna in the patch antenna area of the integrated joined-wing aircraft, obtain the field strength distribution of the patch antenna through the following process:

[0347] For each grid cell corresponding to the patch antenna, calculate the sum of the electric fields carried by the direct, reflected, and diffracted rays on each grid cell to obtain the total electric field of each grid cell on the patch antenna:

[0348] is the total electric field of the th grid cell in the patch antenna, so as to obtain the field strength distribution of the patch antenna; traverse each patch antenna to obtain the field strength distribution of each patch antenna on the integrated joined-wing aircraft.

[0349] Step d2: Based on the field strength distribution of each patch antenna and the field source electric field strength, calculate the electromagnetic shielding rate of the integrated joined-wing aircraft ; specifically:

[0350] Regard the four patch antennas as a region, and according to the field strength distribution of each patch antenna obtained in step c2, perform weighted averaging according to the area size of the grid cells to obtain the average electric field of the entire region , and the specific calculation formula is:

[0351] Among them, represents the total electric field of the th grid cell in the region, represents the area of the th grid cell in the region, represents the sum of the areas of all grid cells in the region; and are obtained by calculating the vertex coordinates of the grid cells;

[0352] Then calculate the electromagnetic shielding rate of the integrated joined-wing aircraft according to the following formula :

[0353] Among them, represents the field source electric field strength.

[0354] Step 5: Determine whether the convergence end condition is met. If not, return to Step 3 to continue the system-level optimization and subsystem-level optimization until the convergence end condition is met, and then end the optimization to obtain the multidisciplinary coupling optimization design result of the integrated joined-wing layout. The convergence end condition is that the change in the objective function value of five consecutive generations of populations is less than 1%, or the maximum number of iterations, 100, is reached.

[0355] In this embodiment, the convergence condition is reached after 57 iterations, and the convergence is good. The main trends of the optimization are as follows:

[0356] As Figure 9 and 10 shown in the comparison diagrams of the aerodynamic / structural and electromagnetic models before and after optimization in this embodiment, where the black color represents the initial configuration before optimization, and the white wireframe represents the configuration after optimization. For the aerodynamic discipline, the relative thickness of the forewing increases, the relative camber of the wing increases, and the chord length of the forewing increases, which is beneficial to improving the lift of the forewing; for the structural discipline, the thickness of the forewing ribs increases, while the thickness of the rear wing beam web, the thickness of the rear wing skin, and the thickness of the rear wing ribs decrease because the initial strength of the rear wing is redundant; for the electromagnetic discipline, the antenna spacing increases because the forewing is the main shielding component, and the forewing area should be avoided as much as possible.

[0357] The following table is a comparison table of the objective functions before and after optimization.

[0358]

[0359] It can be seen that each objective function has been improved to varying degrees before and after optimization. Among them, the lift-to-drag ratio has increased by 10.29%, the mass has been reduced by 8.87%, and the electromagnetic shielding rate has decreased by 21.32%; the aerodynamic objective, the structural objective, and the electromagnetic objective have all been improved significantly, indicating that the design effect of this embodiment is good; in terms of design efficiency, in the multidisciplinary optimization design, the optimization time and the number of design variables usually have a power function or exponential function relationship. Here, is taken. The secondary optimization method of the present invention decomposes a multidisciplinary optimization problem with 55-dimensional design variables into optimization problems with 25-dimensional, 22-dimensional, and two 4-dimensional design variables. Theoretically, the efficiency is increased by 2.65 times; in addition, in the structural discipline analysis, it takes 4.17 minutes to calculate once using the solid element finite element method, and it takes 7 seconds to calculate once using the structural rapid evaluation method of the present invention, and the efficiency is increased by more than 30 times; in terms of calculation accuracy, the calculation result error between the structural rapid evaluation method of the present invention and the solid element finite element method is within 10%. Among them, the deformation error in the main concerned wingtip direction is only -3.21%, which meets the accuracy requirements in engineering.

[0360] The following table is a comparison of the calculation results of the initial configuration node 9 (wingtip section):

[0361]

[0362] In summary, while taking into account both efficiency and accuracy, the present invention can achieve the aerodynamic / structural / electromagnetic multidisciplinary coupling optimization design of the integrated joined-wing layout.

[0363] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Without departing from the principles and spirit of the present invention, those of ordinary skill in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A multidisciplinary coupling optimization design method for integrated joint wing layout, characterized by: The steps include: Step 1: Establishing the initial configuration of the integrated coupled wing layout aircraft and parameterizing the configuration; determining the design objectives, design conditions and design variables of the multidisciplinary optimization of the integrated coupled wing layout aircraft, as well as the design space range of each design variable; Said design goals include high lift-to-drag ratio, low structural mass, and low electromagnetic obstruction; The design variables include aerodynamic design variables, structural design variables, and electromagnetic design variables; Step 2: Divide the design variables determined in step 1 into global design variables and local design variables; Step 3: Perform system-level optimization; The system-level design objectives are high lift-to-drag ratio, low structural mass, and low electromagnetic shielding rate. The global design variables are used as system-level optimization design variables, and the local design variables are fixed as constants. The CFD method is used for aerodynamic target analysis, the equivalent beam method based on engineering beam theory is used for structural target analysis, and the improved electromagnetic shielding rapid evaluation method is used for electromagnetic target analysis. The genetic algorithm is used for optimization to obtain the global optimal design result. Step 4: Perform subsystem-level optimization separately, including aerodynamic optimization, structural optimization, and electromagnetic optimization; The aerodynamic optimization takes high lift-to-drag ratio as the aerodynamic optimization design goal, takes the local design variables in the aerodynamic design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the structural design variables, and the local design variables in the electromagnetic design variables as constants; adopts the CFD method to analyze the aerodynamic target, and uses the genetic algorithm to search for the optimal solution to obtain the local optimal design result of the aerodynamic discipline; The structural optimization takes low structural quality as the structural optimization design goal, takes the local design variables in the structural design variables as the design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic design variables after aerodynamic optimization, and the local design variables in the electromagnetic design variables as constants; combined with the aerodynamic load obtained by the aerodynamic target analysis using the CFD method in the aerodynamic optimization process, the equivalent beam method based on the engineering beam theory is used for structural target analysis, and the genetic algorithm is used for optimization to obtain the local optimal design result of the structural discipline; The electromagnetic discipline optimization takes low electromagnetic shielding rate as the electromagnetic discipline design goal, takes local design variables in the electromagnetic discipline design variables as design variables, and fixes the optimized global design variables obtained in step 3, the local design variables in the aerodynamic discipline design variables after aerodynamic discipline optimization, and the local design variables in the structural discipline design variables after structural discipline optimization as constants; An improved electromagnetic shielding rapid assessment method is used to analyze electromagnetic targets, and a genetic algorithm is used to find the optimal solution to obtain the local optimal design result of the electromagnetic discipline. The improved electromagnetic shielding rapid assessment method is used to perform electromagnetic target analysis, which specifically includes the following steps: Step a2: In the finite element software, for the integrated connected-wing aircraft with the antenna array layout to be calculated, the thickness of the patch antenna is ignored, the aircraft surface is divided into the patch antenna area and other areas, triangular mesh units are generated according to the area, the leading and trailing edges of the wing are meshed, and the vertex coordinates of the mesh units are stored; Step b2: importing the integrated joint wing aircraft in step a2 into electromagnetic analysis software, in the electromagnetic analysis software, setting the field source on the ground, and setting a receiving end on the integrated joint wing aircraft; generating a plane wave as an incident wave in the direction from the field source to the receiving end; the incident wave is composed of a number of rays carrying a set electric field, and the number of rays is equal to the number of grid units in step a2; performing a judgment on the intersection of the ray and the grid unit, if the ray intersects with the grid unit, it is determined that the ray is directed to the bright area; Step c2: Tracing the field intensity of the ray directed to the bright area, and using the improved bouncing ray method to calculate the field intensity distribution of each patch antenna on the integrated linked-wing aircraft; specifically, the following sub-steps are included: Step c2.1: For the ray directed toward the bright area, assuming that the ray propagates in a straight line in a uniform medium, the electric field strength at the intersection point where the ray propagates in a straight line is calculated using the following formula: : in, is the grid cell differential at the source, Differentiate the grid cell at the intersection point; the complex expression Used to describe the propagation characteristics of plane electromagnetic waves in space. is an imaginary unit, is the wave number, is the distance from the intersection point to the field source, is the wavelength of the ray; is the electric field strength of the source; Step c2.2: Determine the intersection type between the ray and the integrated linked-wing aircraft model, which includes intersection in the plane and intersection at the edge; if the ray intersects with the grid unit in the plane, use the geometric optics method to simulate the ray reflection and calculate the ray direction vector after reflection. and the electric field strength of the reflected rays ; If the intersection is at the edge, the improved consistent diffraction theory is used to simulate the ray diffraction and calculate the ray direction vector after diffraction And the electric field strength of the rays after diffraction , the specific process is: First, solve the ray direction vector after diffraction : Get the vertex coordinates of the edge grid unit and its adjacent grid units, and get the normal vector of the edge grid unit according to the vertex coordinates of the edge grid unit ; Obtain the edge vector based on the common edges of the edge grid cell and its adjacent grid cells ; Assume the incident direction vector of the ray is , the incident direction vector and the normal vector The angle between ; The incident direction vector of the ray With edge vector The angle between ; Treat the diffraction phenomenon as the incident ray around the edge vector Deflect an angle , let the direction vector of the ray after diffraction be , the direction vector of the ray after diffraction With normal vector The angle between the two is the diffraction angle ; According to the incident direction vector of the ray and the normal vector Calculate the angle of incidence : According to the formula Calculate the deflection angle ,in is the diffraction coefficient; Assumptions , is the edge vector The three-axis components of the edge vector are calculated by the following expression Rotation The transformation matrix : Will Decomposed into parallel Vector and perpendicular to Vector , the specific expression is: Then the direction vector of the ray after diffraction is for: Then solve the electric field intensity of the ray after diffraction : According to the incident direction vector of the ray and the edge vector calculate : According to the diffraction direction vector of the ray and the normal vector Calculate the diffraction angle : The electric field intensity of the ray after diffraction is calculated by the following formula : In the formula, is the electric field strength at the point where the ray propagates in a straight line. At this time, the intersection point is the diffraction point. is the diffraction coefficient tensor, is the defocus distance, is the distance from the observation point to the diffraction point; The diffraction coefficient tensor is calculated by : The function in the above expression , , and intermediate variables The explanation is as follows: when : when : Step c2.3: trace each ray after reflection and diffraction, and determine the intersection between the ray and the grid unit. If the ray intersects with the grid unit, the ray is directed to the bright area; Step c2.4: Repeat steps c2.1 to c2.3 until the electric field carried by the traced ray is less than the set value or it is directed to a dark area; Step c2.5: For each patch antenna in the patch antenna area of ​​the integrated linked wing aircraft, the field strength distribution of the patch antenna is obtained by the following process: For each grid cell corresponding to the patch antenna, the sum of the electric fields carried by the rays directly, reflected, and diffracted on each grid cell is calculated to obtain the total electric field of each grid cell on the patch antenna: The patch antenna The total electric field of each grid unit is obtained to obtain the field strength distribution of the patch antenna; each patch antenna is traversed to obtain the field strength distribution of each patch antenna area on the integrated joint wing aircraft; Step d2: Based on the field intensity distribution of all patch antenna areas and the electric field strength of the field source, the electromagnetic shielding rate of the integrated linked wing aircraft is calculated. ; Step 5: Determine whether the convergence end condition is met. If not, return to step 3 to continue system-level optimization and subsystem-level optimization until the convergence end condition is met. End the optimization and obtain the multidisciplinary coupling optimization design result of the integrated connection wing layout.

2. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 1 is characterized by: In step 1, the configuration is parameterized and the design variables determined are as follows: the aerodynamic design variables include two-dimensional airfoil design parameters of the front wing and the rear wing and three-dimensional shape layout parameters of the connecting wing; The two-dimensional airfoil adopts the CST parameterization method, where the airfoil function is expressed by multiplying the shape function by the shape function, where the coefficient of the shape function is the design parameter of the two-dimensional airfoil; The 3D shape layout parameters include: Front wingspan , front wing root chord length , the front wing root is slightly larger than , front wing slight twist angle , front wing sweep angle , front wing anhedral angle ; Rear wingspan , rear wing root chord length , the rear wing root is slightly larger than , rear wing slight twist angle , rear wing sweep angle , rear wing anhedral ; Front wing leading edge position , the leading edge position of the rear wing ; End plate height ratio ; in, is the height difference between the front wing and the rear wing at the wing root, is the difference in height between the front wing and the rear wing at the end plate; The structural design variables include structural size parameters and structural layout parameters; Among them, the structural size parameters include: First front spar web thickness: , the thickness of the web of the second front spar: ; The flange area of ​​the first front spar: , the lower flange area of ​​the first front wing spar: ; Second front wing spar upper flange area: , the lower flange area of ​​the second front wing spar: ; Front rib thickness: ; Front wing skin thickness: ; First rear spar web thickness: , the second rear spar web thickness: ; The flange area of ​​the first rear spar: , the lower flange area of ​​the first rear wing spar: ; Second rear wing spar upper flange area: , the lower flange area of ​​the second rear wing spar: ; Rear wing rib thickness: ; Rear wing skin thickness: ; Structural layout parameters include: First front spar position: , second front spar position: ; Front wing rib spacing: ; First rear spar position: , second rear spar position: ; Rear wing rib spacing: ; The electromagnetic design variables include antenna spacing , assuming that m antennas are arranged on the side of the fuselage, then ;in, is the distance between the leading edge of the front wing and the first antenna, is the distance between the first antenna and the second antenna, and so on. is the distance between the m-1th antenna and the mth antenna.

3. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 2 is characterized by: The global design variables in step 2 include all the design parameters of the two-dimensional airfoil in the aerodynamic design variables and the following parameters in the three-dimensional shape layout parameters of the connecting wing: Front wingspan , front wing root chord length , the front wing root is slightly larger than , front wing sweep angle ; Rear wingspan , rear wing root chord length , the rear wing root is slightly larger than , rear wing sweep angle ; Front wing leading edge position , the leading edge position of the rear wing ; End plate height ratio ; The local design variables in step 2 include the aerodynamic local design variables , local design variables of structural disciplines , Electromagnetic discipline local design variables ; Among them, the local design variables of aerodynamics Including the front wing twist angle in the three-dimensional shape layout parameters of the connecting wing , front wing anhedral angle ; Rear wing slightly twisted , rear wing anhedral ; Local design variables of structural disciplines Including all structural size parameters and structural layout parameters; Local Design Variables in Electromagnetic Disciplines Including all antenna spacing .

4. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 1 is characterized by: The mathematical model of system-level optimization in step 3 is as follows: in, The target weight coefficient for lift-to-drag ratio optimization is: Optimize the target weight coefficient for structural quality, is the target weight coefficient of electromagnetic shielding rate; is the lift-to-drag ratio, For structural quality, is the electromagnetic shielding rate, is the lift-to-drag ratio of the initial configuration, is the structural mass of the initial configuration, is the electromagnetic shielding rate of the initial configuration; the constraints include aerodynamic constraints, structural constraints, and electromagnetic constraints, specifically: Aerodynamic constraints: cruise lift coefficient Not less than the lift coefficient setting value , pitch moment coefficient change Not greater than the torque change setting value ; Structural discipline constraints: cruise structure equivalent stress Not more than the allowable stress , maximum wing deformation Do not exceed the deformation setting value , the torsional deformation of the wing under yaw control Do not exceed the torsion angle setting value ; Electromagnetic discipline constraints: electromagnetic shielding rate The electromagnetic shielding rate is not greater than the initial configuration .

5. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 1 is characterized by: The specific mathematical model of aerodynamic optimization in step 4 is as follows: in, is the lift-to-drag ratio; the constraints are: cruise lift coefficient Not less than the lift coefficient setting value , pitch moment coefficient change Not greater than the torque change setting value ; The specific mathematical model of structural optimization is as follows: in, is the structural mass; the specific constraints are: cruise structure equivalent stress Not more than the allowable stress , maximum wing deformation Do not exceed the deformation setting value , the torsional deformation of the wing under yaw control Do not exceed the torsion angle setting value ; The specific mathematical model for electromagnetic optimization is as follows: in, is the electromagnetic shielding rate; the specific constraints are: electromagnetic shielding rate The electromagnetic shielding rate is not greater than the initial configuration .

6. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 1 is characterized by: In step b2, the process of determining whether the ray is directed toward the bright area is as follows: for each ray, traverse each grid unit on the surface of the integrated linked-wing aircraft model: For each grid cell, first perform a single occlusion check: Obtain the ray direction vector based on the source position and the integrated wing aircraft model position , get the normal vector of the mesh cell according to the vertex coordinates of the mesh cell ; like , it is determined to be non-intersecting, the grid unit is marked as a dark area, and the ray is directed to the dark area; if , it is determined to be an intersection, the grid unit is marked as a bright area, and the ray is directed toward the bright area; Then consider the mutual occlusion between grid cells and make multiple occlusion judgments: For a grid cell marked as a bright area, draw a line segment from its center point to the field source, and determine whether this line segment intersects with other grid cells. If the line segment intersects with other grid cells, the grid cell is blocked by other grid cells, and the grid cell is re-marked as a dark area; Finally, all the rays directed toward the bright area are obtained.

7. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 1 is characterized by: The method for determining the intersection type in step c2.2 is that if the angle between the normal of a bright area grid unit and the normal of an adjacent grid unit is less than 10 degrees, the intersection type is determined to be an intersection within the plane; if the angle between the normals is greater than or equal to 10 degrees, the intersection type is determined to be an intersection at the edge, and the bright area grid unit is an edge grid unit.

Images