Multi-disciplinary coupling optimization design method for integrated joined wing layout

Abstract

The invention provides a multidisciplinary coupling optimization design method for integrated joined wing layout. The method comprises the following steps: firstly, determining a design target, a design working condition and a design variable of multi-disciplinary optimization of the joined wing layout; dividing design variables into global design variables and local design variables; taking a high lift-drag ratio, low structural quality and low electromagnetic shielding rate as a system-level design target, taking a global design variable as a design variable, fixing a local design variable, and carrying out system-level optimization; then, high lift-drag ratio, low structural quality and low electromagnetic shielding rate are taken as subsystem-level design targets, local design variables of corresponding subjects are taken as design variables, global design variables and local design variables of other subjects are fixed, and subsystem-level optimization of the subjects is carried out; and finally, judging whether a convergence condition is met or not, and if not, continuing system-level optimization and subsystem-level optimization until convergence. The pneumatic / structural / electromagnetic multi-disciplinary coupling optimization design of the integrated joined wing layout is realized.

Description

Technical Field

[0001] The present invention relates to the field of aircraft design, and in particular to a multidisciplinary coupling optimization design method for integrated linked wing layout. Background Art

[0002] The integrated design of an ultra-long-endurance unmanned aerial platform with a coupled-wing layout often requires a high lift-to-drag ratio, a small structural mass and a small electromagnetic obstruction rate. It is a typical aerodynamic / structure / electromagnetic multidisciplinary coupled optimization design problem. At the same time, the design objectives of multiple disciplines bring high-dimensional design variables. Therefore, the integrated design of an ultra-long-endurance unmanned aerial platform with a coupled-wing layout is a typical high-dimensional multi-objective problem.

[0003] At present, conventional multidisciplinary optimization methods are divided into two categories: holistic and distributed. The holistic approach regards the optimization problem as a holistic solution, such as the multidisciplinary feasible method and the single-discipline feasible method. These methods have high design variable dimensions and are prone to dimensional disasters, resulting in low design efficiency. Distributed approaches decompose the optimization problem into multiple levels of analysis, including collaborative optimization and parallel subspace optimization. The collaborative optimization method decouples the subsystem by introducing auxiliary design variables and consistency constraints, but it has the problem of insufficient convergence and easy premature maturity. The parallel subspace optimization method maintains multidisciplinary consistency through global sensitivity equations, but it needs to solve gradient information during optimization, which has the problem of difficult gradient solution and large amount of calculation.

[0004] In addition, when performing aerodynamic / structural / electromagnetic multidisciplinary optimization on the coupled wing layout, it is necessary to continuously calculate the structural characteristics of the current coupled wing layout scheme, and the structural characteristics analysis methods of traditional configuration aircraft, including the piecewise linear stiffness method and the solid unit finite element method, are not suitable for the integrated coupled wing configuration design. This is mainly because the piecewise linear stiffness method cannot adapt well to the special stiffness configuration of the coupled wing and has the problem of insufficient accuracy, while the solid unit finite element method requires the use of a large number of units and nodes, the analysis is time-consuming, the calculation cost is high, and the optimization design efficiency is low. Summary of the invention

[0005] To solve the above problems, the present invention proposes a multidisciplinary coupling optimization design method for an integrated coupling wing layout. By classifying the design variables and adopting a multidisciplinary system-level optimization and subsystem-level optimization hierarchical design method, the global design variables and local design variables can be reasonably optimized and the subsystem analysis time can be greatly reduced, thereby realizing an efficient multidisciplinary optimization design of the coupling wing layout. In addition, for the integrated coupling wing layout, the structural target calculation adopts a 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 coupling wing. In the electromagnetic discipline optimization, an improved electromagnetic shielding rapid evaluation method is adopted, which can quickly calculate the electromagnetic shielding rate of the antenna array of the integrated coupling wing layout aircraft, thereby providing an efficient and accurate evaluation method for the structural design and aerodynamic / structural / electromagnetic multidisciplinary design of the integrated coupling wing layout.

[0006] The technical solution of the present invention specifically comprises the following steps:

[0007] 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;

[0008] The design objectives include high lift-to-drag ratio, low structural mass and low electromagnetic shielding rate; the design variables include aerodynamic design variables, structural design variables and electromagnetic 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: Perform system-level optimization;

[0011] 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.

[0012] Step 4: Perform subsystem-level optimization separately, including aerodynamic optimization, structural optimization, and electromagnetic optimization;

[0013] 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;

[0014] 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;

[0015] The electromagnetic discipline optimization takes low electromagnetic shielding rate as the electromagnetic discipline design target, 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; adopts the improved electromagnetic shielding rapid evaluation method to perform electromagnetic target analysis, uses genetic algorithm to search for optimization, and obtains the local optimal design result of the electromagnetic discipline;

[0016] 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.

[0017] Furthermore, 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;

[0018] 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;

[0019] The 3D shape layout parameters include:

[0020] 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 ;

[0021] 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;

[0022] The structural design variables include structural size parameters and structural layout parameters;

[0023] Among them, the structural dimension parameters include: the thickness of the web of the first front wing beam: , the thickness of the web of the second front spar: ; Flange area on the first front wing beam: , the lower flange area of ​​the first front wing spar: ; Flange area on the second front wing beam: , the lower flange area of ​​the second front wing spar: ; Front wing rib thickness: ; Front wing skin thickness: ;Thickness of the web of the first rear wing spar: , the second rear spar web thickness: ; Flange area on the first rear wing spar: , the lower flange area of ​​the first rear wing spar: ; Flange area on the second rear wing spar: , the lower flange area of ​​the second rear wing spar: ; Rear wing rib thickness: ; Rear wing skin thickness: ;

[0024] 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: ;

[0025] 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.

[0026] Furthermore, the global design variables XS in step 2 include all 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:

[0027] 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 ;

[0028] The local design variables XL in step 2 include the aerodynamic local design variables , local design variables of structural disciplines , Electromagnetic discipline local design variables ;

[0029] 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 ;

[0030] Local design variables of structural disciplines Including all structural size parameters and structural layout parameters;

[0031] Local Design Variables in Electromagnetic Disciplines Including all antenna spacing .

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

[0033]

[0034]

[0035]

[0036]

[0037] 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:

[0038] 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 ;

[0039] 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 ;

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

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

[0042]

[0043]

[0044] 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 ;

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

[0046]

[0047]

[0048] 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 ;

[0049] The specific mathematical model for electromagnetic optimization is as follows:

[0050]

[0051]

[0052] in, is the electromagnetic shielding rate; the specific constraints are: electromagnetic shielding rate The electromagnetic shielding rate is not greater than the initial configuration .

[0053] Furthermore, in step 3 and step 4, an improved electromagnetic shielding rapid assessment method is used to perform electromagnetic target analysis, which specifically includes the following steps:

[0054] 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;

[0055] 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;

[0056] 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 joint-wing aircraft;

[0057] 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. .

[0058] Furthermore, 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 coupled-wing aircraft model:

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

[0060] 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 ;

[0061] 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;

[0062] Then consider the mutual occlusion between grid cells and make multiple occlusion judgments:

[0063] 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;

[0064] Finally, all the rays directed toward the bright area are obtained.

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

[0066] 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: :

[0067]

[0068] 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;

[0069] 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 ;

[0070] 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;

[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 it is directed to a dark area;

[0072] 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:

[0073] 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:

[0074] The patch antenna The total electric field of the grid units 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 linked 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 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.

[0076] Furthermore, in step c2.2, when the intersection type is intersection 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:

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

[0078] 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 ;

[0079] 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 ;

[0080] According to the incident direction vector of the ray and the normal vector Calculating the angle of incidence :

[0081] According to the formula , calculate the deflection angle ,in is the diffraction coefficient;

[0082] Assumptions , is the normal vector The three-axis components of the edge vector are calculated by the following expression Rotation The transformation matrix :

[0083]

[0084] Will Decomposed into parallel Vector and perpendicular to Vector , the specific expression is:

[0085]

[0086]

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

[0088] Then solve the electric field intensity of the ray 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 electric field intensity of the ray after diffraction is calculated by the following formula :

[0092] 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;

[0093] The diffraction coefficient tensor is calculated by :

[0094]

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

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

[0097] The function in the above expression , , 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 the integrated coupling wing layout, which mainly realizes the aerodynamic / structural / electromagnetic multidisciplinary coupling optimization of the coupling wing. In view of the problems of low efficiency, insufficient convergence, difficult gradient solution, low precision, etc. in the prior art, the present invention uses a multidisciplinary two-level optimization method to divide the design variables into global design variables and local design variables of each sub-discipline, which greatly reduces the optimization design dimension and solves the problem of low efficiency; the design space at the system level and the subsystem level does not overlap, and there is no need to introduce coupling auxiliary design variables and state coupling design variables, and the convergence is good; when classifying the design variables, the division is based on the coupling relationship of the various disciplines of the coupling wing, and there is no need to perform gradient solution, so the calculation amount is small; in addition, for the integrated coupling wing layout, the structural target calculation adopts finite element based on engineering beam theory, etc. The effective beam method, on the one hand, combines the classical structural force transmission model of engineering beam theory, and on the other hand, adds the finite element beam model, which can be used to analyze the special stiffness configuration of the connecting wing; in addition, in the optimization of electromagnetic discipline, an improved electromagnetic shielding rapid evaluation method is adopted, which is improved on the basis of the bouncing ray method. On the one hand, the diffraction coefficient is directly defined to quickly solve the diffraction angle; on the other hand, the third-order diffraction tensor coefficient in the traditional bouncing ray method is reduced in dimension to avoid magnetic field calculation, and the electric field of the diffraction ray can be quickly solved, which can provide an efficient and accurate evaluation method for the structural design and aerodynamic / structural / electromagnetic multidisciplinary design of the integrated connecting wing layout. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0106] Figure 2 The initial configuration of the connecting wing of the embodiment of the present invention;

[0107] Figure 3 A simplified electromagnetic analysis three-dimensional model of the connecting wing according to an embodiment of the present invention;

[0108] Figure 4 A detailed flow chart of a method for rapid equivalent analysis of structural characteristics according to an embodiment of the present invention;

[0109] Figure 5 A simplified three-dimensional model for analyzing the connecting wing structure according to an embodiment of the present invention;

[0110] Figure 6 It is a schematic diagram of the division of the connecting wing section of an embodiment of the present invention;

[0111] Figure 7 is a detailed flow chart of the electromagnetic shielding rate analysis method according to an embodiment of the present invention;

[0112] Figure 8 A schematic diagram of a grid for analyzing the electromagnetic shielding rate according to an embodiment of the present invention;

[0113] Figure 9 A comparison diagram of the aerodynamic / structural three-dimensional models of the front and rear connecting wings optimized in an embodiment of the present invention;

[0114] Figure 10 A comparison diagram of the electromagnetic three-dimensional models of the front and rear connection wings optimized in an embodiment of the present invention;

[0115] Among them: 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 spar; 10-second front wing spar; 11-front wing rib; 12-front wing skin; 13-first rear wing spar; 14-second rear wing spar; 15-rear wing rib; 16-rear wing skin; 17-single-side 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; DETAILED DESCRIPTION

[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 solutions of the present invention, the present invention is further described and fully described in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0117] This embodiment takes a double-beam coupled wing ultra-long flight time hovering UAV as an example. The initial configuration of the coupled wing aircraft is as follows: Figure 2As shown in FIG. 1 , the configuration is mainly composed of a fuselage, a front wing, a rear wing, a connecting end plate and a tail wing. Propellers for vertical take-off and landing are arranged at both ends of the two connecting end plates. There are a total of 4 propellers. For the sake of simplicity in the diagram Figure 2 The following adopts the multidisciplinary coupling optimization design method for the integrated link wing layout proposed by the present invention to perform aerodynamic, structural, and electromagnetic multidisciplinary optimization design on the link wing aircraft; in order to simplify the problem, a half-mold is taken for the design of the double-beam link wing ultra-long flight time hovering unmanned platform. Taking the electromagnetic analysis model as an example, the simplified electromagnetic analysis model is as follows Figure 3 As shown, the simplified model consists of a fuselage, connecting end plates, front wings and rear wings. The antenna layout is as follows: four patch antennas are arranged at equal intervals on the side of the fuselage, and the thickness of the antenna is ignored.

[0118] like Figure 1 As shown, the implementation flow chart of the scheme of the embodiment of the present invention specifically includes the following steps:

[0119] 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;

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

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

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

[0123] Vertical take-off and landing conditions: including hovering state and yaw control state. In the hovering state, the four propellers have the same pulling force, that is, ,in, is the left front propeller thrust, is the right front propeller thrust, is the left rear propeller thrust, is the thrust of the right rear propeller; in the yaw control state, the diagonal rotor increases to 1.5 times the lift, that is, , .

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

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

[0126] 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;

[0127] In this embodiment, only the front wing airfoil is optimized; the front wing 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, For class functions, is the shape function;

[0129]

[0130]

[0131]

[0132] In the formula, , In order to control the type parameters of the geometric shape, in this embodiment, , ; Coefficients of shape functions is the weight coefficient, is the Bernstein polynomial, n is its order, and in this embodiment, n=6;

[0133] That is, the design parameters of the two-dimensional airfoil. In this embodiment, the two-dimensional airfoil is parameterized by dividing the upper airfoil and the lower airfoil, and the weight coefficient is used. - Represents the design parameters of the upper wing surface, weight coefficient - represents the design parameters of the lower wing surface;

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

[0135]

[0136] The 3D shape layout parameters include:

[0137] 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 ;

[0138] 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;

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

[0140]

[0141] The structural design variables include structural size parameters and structural layout parameters;

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

[0143] First front spar web thickness: , the thickness of the web of the second front spar: ; Flange area on the first front wing beam: , the lower flange area of ​​the first front wing spar: ; Flange area on the second front wing beam: , the lower flange area of ​​the second front wing spar: ; Front wing rib thickness: ; Front wing skin thickness: ;Thickness of the first rear spar web: , the second rear spar web thickness: ; Flange area on the first rear wing spar: , the lower flange area of ​​the first rear wing spar: ; Flange area on the second rear wing spar: , the lower flange area of ​​the second rear wing spar: ; Rear wing rib thickness: ; Rear wing skin thickness: ; Structural layout parameters include: Position of the first front spar: , second front spar position: ; Front wing rib spacing: ; First rear spar position: , second rear spar position: ; Rear wing rib spacing: ;

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

[0145]

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

[0147]

[0148] The electromagnetic design variables include 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 ;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 second antenna.

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

[0150]

[0151] In this embodiment, based on the above design variables, the specific method of parameterizing the configuration is:

[0152] The 3D coupled wing configuration is parameterized by secondary development of CATIA software, including parametric modeling of aerodynamic / structural / electromagnetic disciplines. Specifically, when the digital model of the 3D coupled wing configuration is established by CATIA, the leading edge position coordinates of the wing root and the wing tip, the ratio of the wing root and the wing tip, and the rotation angle of the wing tip section are used as control parameters of the aerodynamic design variables. 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 are used as control parameters of the structural design variables. The coordinates of the leading edge of the antenna are used as control parameters of the electromagnetic design variables. A script file is recorded, and the changes of the 3D coupled wing configuration are controlled by modifying the above control parameters in the script file.

[0153] Among them, in aerodynamics, the span of the front and rear wings, the front wing sweep angle, the front wing dihedral, the rear wing sweep angle, the rear wing dihedral, the leading edge position 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 the wing tip; the root chord length and root-tip ratio of the front and rear wings are controlled by changing the ratio of the wing root and the wing tip; the wing tip torsion angle of the front and rear wings is controlled by changing the rotation angle of the wing tip section.

[0154] In the structural discipline, the skin thickness of the front and rear wings is controlled by changing the thick surface size of the skin, the web thickness and upper and lower flange areas of the front and rear wings are controlled by changing the sketch size of the beam section, the position of the front and rear wing spars is controlled by changing the coordinates of the beam section, the thickness of the front and rear wing ribs is controlled by changing the thick surface size of the ribs, and the rib spacing of the front and rear wings is controlled by changing the translation matrix size of the ribs.

[0155] Electromagnetic science controls the antenna spacing 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 aerodynamic design variables and the following parameters in the three-dimensional shape layout parameters of the connecting wing:

[0158] 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 ;

[0159] Local Design Variables XL includes local design variables for aerodynamics , local design variables of structural disciplines , Electromagnetic discipline local design variables ;

[0160] 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 ;

[0161] Local design variables of structural disciplines Including all structural size parameters and structural layout parameters;

[0162] Local Design Variables in Electromagnetic Disciplines Including all antenna spacing .

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

[0164] Specifically: high lift-to-drag ratio, small structural mass, and small electromagnetic shielding rate are set as system-level design goals, global design variables are used as system-level optimization design variables, and local design variables are fixed as constants; CFD method is used for aerodynamic target analysis, equivalent beam method based on engineering beam theory is used for structural target analysis, and improved electromagnetic shielding rapid evaluation method is used for electromagnetic target analysis; genetic algorithm is used for optimization to obtain the global optimal design result;

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

[0166]

[0167]

[0168]

[0169]

[0170] 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, which is 1 in this embodiment; They are lift-to-drag ratio, structural mass, and electromagnetic shielding rate. They are the lift-to-drag ratio, structural mass, and electromagnetic shielding rate of the initial configuration; the constraints include aerodynamic constraints, structural constraints, and electromagnetic constraints, specifically:

[0171] Aerodynamic constraints: cruise lift coefficient Not less than the lift setting value In this embodiment, , pitch moment coefficient change Not greater than the torque change setting value In this embodiment, ;

[0172] Structural discipline constraints: cruise structure equivalent stress Not more than the allowable stress , maximum wing deformation Do not exceed the deformation setting value In this embodiment, , unit: meter, wing torsional deformation under yaw control state Do not exceed the torsion angle setting value In this embodiment, ;

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

[0174] Step 4: Perform subsystem-level optimization separately, including aerodynamic optimization, structural optimization, and electromagnetic optimization;

[0175] Specifically: The aerodynamic optimization takes high lift-to-drag ratio as the aerodynamic optimization design target, 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 of the connected wing layout;

[0176] In this embodiment, the pneumatic optimization mathematical model is as follows:

[0177]

[0178]

[0179] in, is the lift-to-drag ratio; the constraints are: cruise lift coefficient Not less than the lift setting value In this embodiment, , pitch moment coefficient change Not greater than the torque change setting value In this embodiment, ;

[0180] 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 structural local optimal design result of the connection wing layout;

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

[0182]

[0183]

[0184] 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 In this embodiment, , the torsional deformation of the wing under yaw control Do not exceed the torsion angle setting value In this embodiment, ;

[0185] The electromagnetic optimization takes low electromagnetic shielding rate as the electromagnetic design target, takes local design variables in the electromagnetic design variables as 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 structural design variables after structural optimization as constants; adopts the improved electromagnetic shielding rapid evaluation method to perform electromagnetic target analysis, uses genetic algorithm to search for optimization, and obtains the electromagnetic local optimal design result of the connection wing layout;

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

[0187]

[0188]

[0189] in, is the electromagnetic shielding rate; the specific constraints are: electromagnetic shielding rate The electromagnetic shielding rate is not greater than the initial configuration .

[0190] In this embodiment, if Figure 4 As shown in step 3 and step 4, the equivalent beam method based on engineering beam theory is used to perform structural target analysis to obtain the cruise structure equivalent stress of the coupled wing configuration. , Maximum wing deformation , Wing torsional deformation under yaw control The specific steps for determining the overall structural quality are as follows. The equivalent beam method based on engineering beam theory can be found in Publication No. The Chinese patent application "Rapid equivalent analysis method of structural characteristics for integrated joint wing configuration design":

[0191] First, establish the global coordinate system , the origin of the global coordinate system Based on the design focus of the integrated wing layout aircraft, The axis points to the rear of the aircraft along the longitudinal symmetric section chord of the aircraft. The axis lies in the section plane perpendicular to Axially, The axis vertical section faces inwards of the aircraft, satisfying the right-hand system;

[0192] Step a1: Obtain the design variables of the coupling wing configuration scheme to be analyzed, regard the coupling wing to be analyzed as a spatial rigid frame structure composed of a front wing, a rear wing and a 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, not a flexible structure;

[0193] In this embodiment, Figure 2 The connecting wing shown in the figure is regarded as a spatial rigid frame structure composed of the front wing, the rear wing and the connecting end plate, which is simplified into a spatial six-degree-of-freedom rigid frame structure model; the simplified structural analysis three-dimensional model is as follows Figure 5 shown.

[0194] Step b1: In the spatial six-degree-of-freedom rigid frame structure model, the front wing, the rear wing and the connecting end plate are divided into a number of sections at set intervals and equal spacing, and the leading edge point coordinates of each section in the global coordinate system are obtained; wherein the sections of the front wing and the rear wing are parallel to the longitudinal symmetry plane of the aircraft where the connecting wing is located, and the section of the connecting end plate is perpendicular to the longitudinal symmetry plane of the aircraft where the connecting wing is located;

[0195] In this embodiment, the spatial six-degree-of-freedom rigid frame structure model is divided into sections according to the following two principles. The division diagram is shown in FIG. Figure 6 As shown:

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

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

[0198] The method for obtaining the coordinates of the leading edge points of each section is to first obtain the coordinates of the leading edge points of the sections at both ends of the front wing, rear wing and connecting end plate and the sections at the connection 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 remaining front edge points of the cross section are calculated through geometric relationships.

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

[0202] Let's take the front wing root section as an example, that is, section #1, and calculate its centroid coordinates. And get its centroid coordinate system , which specifically includes the following sub-steps:

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

[0204] In the formula, is the reduction factor or conversion factor, which represents 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 factor is ; is the skin thickness of the corresponding area of ​​the cross section. For example, in this embodiment, the cross section is the root cross section of the front wing. is the front wing skin thickness , is the arc length of the skin; is the thickness of the beam web in the corresponding area of ​​the cross section. In this embodiment, is the front spar web thickness , is the arc length of the beam web; is the kth type of spar flange area in the section corresponding region. In this embodiment, there are 4 types of spar flanges in the section corresponding region of the front wing root, which are the upper flange area of ​​the first front spar , the lower flange area of ​​the first front spar , flange area on the second front spar and the lower flange area of ​​the second front spar .

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

[0206]

[0207]

[0208] In the formula, is the Z-coordinate of the arc length element of the skin in the corresponding area of ​​the cross section in the global coordinate system, is the Z-coordinate of the arc length of the beam web in the corresponding area of ​​the section in the global coordinate system, is the Z-coordinate of the k-th wing spar flange in the corresponding area of ​​the cross section in the global coordinate system; is the X-coordinate of the skin arc length element in the corresponding area of ​​the cross section in the global coordinate system, is the X-coordinate of the arc length of the beam web in the corresponding area of ​​the section in the global coordinate system, is the X-coordinate of the kth wing spar 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 , get the centroid coordinate system corresponding to the section;

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

[0213] The details are as follows: Since the angle between the centroid coordinate system and the section principal inertia axis coordinate system is small, the centroid coordinate system is regarded as the principal inertia axis coordinate system, so the moment of inertia of the section divided by step b1 to the principal inertia axis, that is, the centroid axis, is and It is obtained by the following calculation formula:

[0214]

[0215] ;

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

[0217] The specific description is as follows, taking section #1 as an example:

[0218] Using engineering beam theory, the wing is regarded as a beam supported on the fuselage. The moment of inertia of the cross-section obtained in step d1 in the centroid coordinate system is and , the coordinates of the rigid center of the section It is obtained by the following calculation formula:

[0219]

[0220]

[0221] in, 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 vertical distance from the shear force of the skin arc length element in the corresponding area of ​​the section to the reference point of the set moment, The vertical distance from the shear force of the beam web arc length microelement in the corresponding area of ​​the section to the set moment reference point It is twice the area contained by the closed curve of the single closed chamber corresponding to the skin arc length microelement of the section. It is twice the area contained by the arc length of the beam web in the corresponding area of ​​the section corresponding to the closed curve of the single closed chamber. is the shear modulus of the reference material;

[0222] Move the origin of the centroid coordinate system to the rigid center coordinate system , get the rigid center coordinate system corresponding to the section;

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

[0224] Specifically, the following is explained taking section #1 as an example:

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

[0226]

[0227]

[0228] According to the formula

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

[0230] Step g1: The X-axis moment of inertia of each section in its respective rigid center coordinate system obtained in step f1 , Z-axis moment of inertia and torsional stiffness , establish the unit overall stiffness matrix of each unit and assemble to obtain the structural overall stiffness matrix of the spatial six-degree-of-freedom rigid frame structure model ; Specifically includes the following sub-steps:

[0231] Step g1.1: Treat each section as a node. The location of the node is the rigid center of the section. and nodes The local coordinate system of the unit is based on the node is the origin, and the coordinate axes are axis, axis, Axis, where Axis as Node and nodes Connection, Axis and Node The coordinate system of the rigid center of the section Axis coincidence, Axis and Axis and The axes form a right-handed coordinate system;

[0232] The simplified configuration of the connecting wing is a spatial six-degree-of-freedom rigid frame, so each node has 6 degrees of freedom, each unit has 12 degrees of freedom, and each unit is regarded as a beam unit with equal cross-section. According to the inertia moment and torsional stiffness of each section in the respective rigid center coordinate system obtained in step f1, the node and nodes The element local stiffness matrix of the constituent elements :

[0233]

[0234] In the formula, is Young's modulus, is the shear modulus, is the total cross-sectional area, and is the moment of inertia of the section rigid center coordinate system, , is the unit length, i.e. the distance between two sections; columns 1-6 of the matrix correspond to nodes of Direction displacement ( )、 Direction displacement ( )、 Direction displacement ( )、 Direction Angle ( )、 Direction Angle ( )、 Direction Angle ( ); Similarly, columns 7-12 of the matrix correspond to nodes 6 degrees of freedom;

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

[0236] The local stiffness matrix of the element in the local coordinate system is The element global stiffness matrix converted to the global coordinate system ,in

[0237]

[0238]

[0239] in is the coordinate transformation matrix, is the Euler transformation matrix, is the unit global stiffness matrix, It is the unit local coordinate system to the global coordinate system Euler angles of three-axis rotation;

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

[0241] Specifically, the unit global stiffness matrix of each unit is assembled into the structural global stiffness matrix in sequence according to the node number of the unit. The node number determines its position in the structural global stiffness matrix.

[0242] In this embodiment, the structural material properties of the skin, ribs and spars are combined to obtain the overall structural stiffness matrix The matrix size is The level is no longer displayed;

[0243] Step h1: According to the overall structural stiffness matrix obtained in step g1, the structural equation is established and solved to obtain the cruise structural equivalent stress of the connected wing configuration. , Wing deformation , Wing torsional deformation under yaw control and the overall structural quality, which specifically includes the following sub-steps:

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

[0245] 1) Load handling:

[0246] According to the aerodynamic load obtained by aerodynamic target analysis using CFD method in the aerodynamic optimization process, the node Load , the node Load It is expressed as:

[0247] In the formula, For Node Force loads in three directions in the global coordinate system, For Node Moment loads in three directions in the global coordinate system;

[0248] Assembling load vectors , expressed as:

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

[0250] Assemble to obtain load vector After that, combined with the overall stiffness matrix obtained in step g1 , the structural equation is established ,in is the displacement vector matrix to be solved, which contains displacements and rotations in three directions, namely Direction displacement, Direction displacement, Directional displacement; Direction angle, Direction angle, Direction angle;

[0251] 2) Boundary condition processing:

[0252] The displacement vector The positions of the nodes corresponding to the wing roots of the front and rear wings are limited to 0; and the displacement vectors of the nodes corresponding to the connecting end plates are limited to 0;

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

[0254] 1) By structural equation Available , the deflection and rotation angle of each node are obtained; the deflection and rotation angle of each node include the following six parameters:

[0255] Direction displacement: , Direction displacement: , Direction displacement: ; Direction Angle: , Direction Angle: , Direction Angle: ;

[0256] 2) According to the deflection and rotation angle of each node, the shape function of the deflection and rotation angle is interpolated for each unit, and the deflection and rotation angle of any point in the unit are obtained, and then the wing deformation of the connected wing configuration is obtained. , Wing torsional deformation under yaw control , as follows:

[0257] X-direction displacement shape function:

[0258] In the formula Is the unit in position The displacement along the X direction, is the unknown coefficient, obtained by interpolation, the Y-direction displacement shape function:

[0259] In the formula Is the unit in position The displacement along the Y direction, is the unknown coefficient, obtained by interpolation, the Z-direction displacement shape function:

[0260] In the formula Is the unit in position The displacement along the Z direction, The unknown coefficients are obtained by interpolation, and the shape function of the rotation angle around the X-axis is:

[0261] In the formula Is the unit in position The rotation angle around the X axis, The coefficients to be determined are obtained by interpolation, and the shape function of the rotation angle around the Y axis is:

[0262] In the formula Is the unit in position The rotation angle around the Y axis, is the unknown coefficient, obtained by interpolation, the shape function of the rotation angle around the Z axis:

[0263] In the formula Is the unit in position The rotation angle around the Z axis, is the unknown coefficient, obtained by interpolation, and the wing deformation of the connected wing configuration is obtained ; Wing torsional deformation under yaw control state ;

[0264] 3) Based on the obtained displacement shape functions in the X, Y, and Z directions, the equivalent stress of the unit structure of the connecting wing configuration is calculated , as follows:

[0265] According to the displacement shape function, the linear strain along the X-axis, Y-axis, and Z-axis is obtained:

[0266]

[0267] In the formula is the linear strain along the X-axis, is the linear strain along the Y-axis, is the linear strain along the Z axis;

[0268] According to the displacement shape function, the shear strain along the XY plane, YZ plane, and ZX plane is obtained:

[0269] , ,

[0270] In the formula 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 strain, the normal stress is calculated as follows:

[0272]

[0273]

[0274]

[0275] In the formula, is the Poisson's ratio of the material;

[0276] Based on the obtained shear strain, the shear stress is calculated 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: ; Get the eigenvalue , that is, the three principal stresses of the element , , , is the identity matrix;

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

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

[0283] 4) Calculate the structural mass of the entire connecting wing configuration based on the parameters of each section , the specific calculation formula is as follows:

[0284] m = ∑ i=1 n [D x i * A i+1 + A i 2 ] r mean

[0285] In the formula, , Node The total structural area of ​​the section, For Node The distance, n is the number of nodes, is the average density of the cross section.

[0286] In this embodiment, if Figure 7 As shown in step 3 and step 4, the improved electromagnetic shielding rapid assessment method is used to analyze the electromagnetic target and obtain the electromagnetic shielding rate of the coupled wing configuration. The specific steps are as follows:

[0287] 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;

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

[0289] 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;

[0290] The process of judging whether the ray is directed to the bright area is specifically as follows: for each ray, traverse each grid unit on the surface of the integrated joint wing aircraft;

[0291] For each grid cell, first perform a single occlusion check:

[0292] Obtain the ray direction vector based on the source position and the integrated wing aircraft position , get the normal vector of the mesh cell according to the vertex coordinates of the mesh cell ;

[0293] In this embodiment, in global coordinates , assuming that the integrated joint wing aircraft is at an altitude of 1000m and the source is 2000m to the left of the aircraft, then the receiving end coordinates are , field source coordinates , the ray direction vector After normalization, , set the initial electric field strength of the ray (field source electric field strength) to .

[0294] 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;

[0295] Then consider the mutual occlusion between grid cells and make multiple occlusion judgments:

[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, the grid cell is blocked by other grid cells, and the grid cell is re-marked as a dark area;

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

[0298] 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:

[0299] 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: :

[0300]

[0301] 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;

[0302] Step c2.2: Determine the intersection type between the ray and the integrated linked wing aircraft, which includes intersection in 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 ;

[0303] In this embodiment, the method for determining the intersection type is as follows: 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, it is determined that the intersection is in a plane; if the angle between the normals is greater than or equal to 10 degrees, it is determined that the intersection is at an edge, and the bright area grid unit is an edge grid unit.

[0304] ① When the intersection is in a plane, the geometric optics method (GO) is used to simulate the reflection of the ray and calculate the direction vector of the ray after reflection. and the electric field strength of the reflected rays The specific process is:

[0305] First solve the direction vector of the ray after reflection :

[0306] Get the normal vector of the mesh cell according to the vertex coordinates of the mesh cell ;

[0307] Assume the incident direction vector of the ray is , the incident direction vector of the ray and the normal vector of the grid cell The angle between ; Let the direction vector of the reflected ray be , the direction vector of the ray after reflection and the normal vector of the grid cell The angle of reflection is ; According to the law of reflection, the angle of incidence Equal to the reflection angle ; In the incident direction vector and the normal vector In the plane composed of Decomposed into normal vector Parallel vectors Sum and normal vector Vertical Vector ;

[0308] in That is exist The projection vector on

[0309] ;

[0310]

[0311] because ,so and About Normal Vectors Symmetrical, that is

[0312] get The expression is:

[0313] Then the electric field intensity of the reflected ray is solved by the following formula :

[0314] in, is the electric field strength at the intersection point where the ray straight line propagates from step 31. At this time, the intersection point is the reflection point. represents the reflection coefficient, is the grid cell differential 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, the improved consistent diffraction theory is used to simulate the diffraction of the ray and calculate the direction vector of the ray after diffraction. And the electric field strength of the rays after diffraction The specific process is:

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

[0317] 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 ;

[0318] The main shielding parts of the integrated wing aircraft are the connecting end plates and the front and rear wings. The diffraction problem is simplified to the flat plate diffraction problem. The incident direction vector of the ray is , the normal vector of the plate is taken as , and take the edge vector of the plate as , the incident direction vector of the ray and the normal vector The angle between ; The incident direction vector of the ray and the 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 and the normal vector The angle between the two is the diffraction angle ;

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

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

[0321] Assumptions , is the normal vector The three-axis components of the edge vector are calculated by the following expression Rotation The transformation matrix

[0322]

[0323] Will Decomposed into parallel Vector and perpendicular to Vector , the specific expression is:

[0324]

[0325]

[0326] Then the direction vector of the ray after diffraction is for:

[0327] Then solve the electric field intensity of the ray after diffraction :

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

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

[0330] According to the Keller cone, the electric field intensity of the ray after diffraction is calculated by the following formula: :

[0331] 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; As mentioned above, the main shielding components of the integrated joint-wing aircraft are the connecting end plates and the front and rear wings, simplifying the diffraction problem to a flat plate diffraction problem. For the flat plate diffraction problem, we have , ;

[0332] The diffraction coefficient tensor in the conventional bouncing ray method To reduce the dimension, we do not consider the spatial dimension characteristics of the electromagnetic field, but only consider the electric field solution, and calculate the diffraction coefficient tensor by the following formula: :

[0333]

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

[0335] D h = - e -jπ / 4 2π 2πbsin β 0 [ F bLa i s - I will i cos i s - I will i 2 + F( bLa( i s + I will i ) ) cos (( i s + I will 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 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;

[0344] 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 shoot towards the dark area;

[0345] In this embodiment, the setting value , that is, when the electric field intensity of the reflected or diffracted ray is less than When , it is considered that the ray is emitted 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 linked wing aircraft, the field strength distribution of the patch antenna is obtained by the following process:

[0347] 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:

[0348] The patch antenna The total electric field of the grid units 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 on the integrated linked-wing aircraft.

[0349] Step d2: Based on the field intensity distribution of each patch antenna and the electric field strength of the field source, the electromagnetic shielding rate of the integrated linked wing aircraft is calculated. ; Specifically:

[0350] Consider the four patch antennas as one area, and perform weighted averaging based on the grid unit area size according to the field strength distribution of each patch antenna obtained in step c2 to obtain the average electric field of the entire area. , the specific calculation formula is:

[0351] in, Indicates the area The total electric field of each grid cell is Indicates the area The area of ​​a grid cell, Represents the sum of the areas of all grid cells in the region; and Obtained through calculation of vertex coordinates of mesh cells;

[0352] Then the electromagnetic shielding rate of the integrated wing aircraft is calculated according to the following formula: :

[0353] in, Represents the electric field strength of the field source.

[0354] 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. The convergence end condition is that the objective function value of the population changes by less than 1% for five consecutive generations, or the maximum number of iterations reaches 100.

[0355] In this embodiment, the convergence condition is reached after 57 iterations, the convergence is good, and the main trend of the optimization is:

[0356] like Figure 9 and 10 The figure shows the comparison of the aerodynamic / structural and electromagnetic models before and after the optimization of this embodiment, where the black is the initial configuration before optimization and the white wireframe is the configuration after optimization. For the aerodynamic subject, the relative thickness of the front wing increases, the relative curvature of the wing increases, and the chord length of the front wing increases, which is beneficial to improving the lift of the front wing; for the structural subject, the thickness of the front wing rib increases, and the thickness of the rear wing beam web, the thickness of the rear wing skin, and the thickness of the rear wing rib decrease, which is due to the excess initial strength of the rear wing; for the electromagnetic subject, the antenna spacing increases, which is because the front wing is the main shielding component and the front wing area should be avoided as much as possible.

[0357] The following table is a comparison table of the objective function 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 which the lift-to-drag ratio has been increased by 10.29%, the mass has been reduced by 8.87%, and the electromagnetic shielding rate has been reduced by 21.32%; the aerodynamic objectives, structural objectives, and electromagnetic objectives have all been well improved, indicating that the design effect of this embodiment is good; in terms of design efficiency, the optimization time in the multidisciplinary optimization design is and the number of design variables The relationship is usually a power function or exponential function. Here we take The secondary optimization method of the present invention decomposes a multidisciplinary optimization problem of 55-dimensional design variables into an optimization problem of 25-dimensional, 22-dimensional, and two 4-dimensional design variables, which theoretically improves the efficiency by 2.65 times. In addition, in the structural discipline analysis, the solid unit finite element method takes 4.17 minutes to calculate once, and the structure rapid evaluation method of the present invention takes 7 seconds to calculate once, which improves the efficiency by more than 30 times. In terms of calculation accuracy, the calculation results of the structure rapid evaluation method of the present invention and the solid unit finite element method have an error within 10%, among which the wing tip is the main concern. The directional deformation error is only -3.21%, which meets the accuracy requirements of the project.

[0360] The following table compares the calculation results of the initial configuration node 9 (wing tip section):

[0361]

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

[0363] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention without departing from the principles and intent 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. 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 thickness of the web of the second rear spar: ; 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 3 and step 4, an improved electromagnetic shielding rapid assessment method is used to analyze electromagnetic targets, 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 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. .

7. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 6 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.

8. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 6 is characterized by: Step c2 specifically includes the following sub-steps: 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 ; Step c2.3: trace each ray after reflection and diffraction, and determine if the ray intersects with the grid unit. If the ray intersects with the grid unit, the ray is directed toward 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 the grid units 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 linked wing aircraft.

9. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 6 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.

10. The multidisciplinary coupling optimization design method for integrated joint wing layout according to claim 6 is characterized by: In step c2.2, when the intersection type is intersection 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 normal 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 : 。

Images