Rapid static aeroelasticity analysis method for unmanned aerial vehicle with high aspect ratio

Abstract

The invention discloses a rapid static aeroelasticity analysis method for a high-aspect-ratio unmanned aerial vehicle, and the method comprises the steps: firstly, carrying out the measurement of a geometric model, so as to construct an aerodynamic model in the form of a standard data card, and building a simplified beam model of the high-aspect-ratio unmanned aerial vehicle; carrying out aerodynamic solution based on a curved surface vortex grid method and statics solution based on a classical beam theory; and loading an aerodynamic load obtained by solving aerodynamic force into the structural model according to a force translation principle, solving torsion and bending of a beam through a beam theory, mapping structural parameters into aerodynamic parameters, and repeating the operation until calculation is converged. And finally, outputting final displacement, aerodynamic distribution and other results. Compared with a calculation result based on a CFD / CSD high-precision numerical method, the static aeroelasticity calculation result has the advantages that the maximum error of elastic deformation does not exceed 0.3%, and the precision is high.

Description

Technical Field

[0001] The present invention belongs to the technical field of unmanned aerial vehicles, and particularly relates to a method for rapidly analyzing the static aeroelasticity of a high aspect ratio unmanned aerial vehicle. Background Art

[0002] High altitude long endurance flying wing layout unmanned aerial vehicles often have design characteristics such as high aspect ratio and lightweight structure. Under the action of aerodynamic loads, their wings will produce large bending and torsional deformations. The elastic deformation of this structure will cause the redistribution of aerodynamic loads on the wings, thereby having a greater impact on the aerodynamic characteristics of the whole aircraft. Therefore, analyzing the aeroelastic deformation and load redistribution of high aspect ratio unmanned aerial vehicles is of great significance for strength verification. Domestic and foreign scholars have carried out relatively more research on geometrically nonlinear aeroelastic problems with high aspect ratio flexible wings as the object. Smith et al. and Garcia et al. respectively studied the static aeroelastic characteristics of high aspect ratio flexible wings based on the geometrically exact intrinsic beam theory and the three-dimensional geometrically nonlinear beam theory, combined with the Euler solver. Garcia studied the static aeroelasticity of high aspect ratio straight wings and swept wings at transonic speeds, and explored the relationship between transonic drag and structural bending-torsion coupling. The Zhou Zhou team at Northwestern Polytechnical University in China has also carried out relevant research. By writing a coupled solver for computational structural mechanics and computational fluid dynamics, the static aeroelastic problems of unmanned aerial vehicles with a layout similar to "Helios" have been studied. In the 1970s and 1980s, Hodges, Dowell et al. started from the accurate description of the motion and deformation of a beam with initial bending and torsion deformation, and derived the motion equation of an accurately large deformation geometrically nonlinear beam, namely the Hodges-Dowell equation, and thus established a mixed variational formula and a finite element algorithm. Tran and Petot et al. from the French Aerospace Research Institute first proposed a semi-empirical, unsteady, nonlinear two-dimensional aerodynamic force model, which has been improved to form the currently recognized nonlinear aerodynamic force model, called the ONERA model. Since the beginning of this century, Chinese aviation scientific and technical personnel have also carried out follow-up research on the geometrically nonlinear aeroelastic calculation methods based on the geometrically nonlinear beam and the ONERA nonlinear aerodynamic force model.

[0003] At present, in order to perform the static aeroelastic analysis of high aspect ratio unmanned aerial vehicles, CFD / CSD (Computational Fluid Dynamics / Computational Structural Dynamics) numerical simulation technology is generally directly used for fluid-structure coupling calculation. However, its calculation complexity and high calculation cost are not conducive to engineering design. Therefore, it is particularly important to develop a static aeroelastic method based on analytical or semi-analytical methods to achieve the rapid static aeroelastic analysis of high aspect ratio unmanned aerial vehicles. Summary of the Invention

[0004] To overcome the deficiencies of the prior art, the present invention provides a method for rapidly analyzing the static aeroelasticity of a high aspect ratio unmanned aerial vehicle. First, the geometric model is measured to construct an aerodynamic model in the form of a standard data card and establish a simplified beam model of the high aspect ratio unmanned aerial vehicle; then, the aerodynamic force solution based on the surface panel method and the static mechanics solution based on the classical beam theory are carried out; the aerodynamic loads obtained from the aerodynamic force solution are loaded into the structural model according to the principle of force translation, the torsion and bending of the beam are solved through the beam theory, and the structural parameters are mapped into aerodynamic parameters, and the above operations are repeated until the calculation converges. Finally, the results such as the final displacement and aerodynamic force distribution are output. The static aeroelastic calculation results of the present invention are compared with the calculation results based on the high-precision numerical method of CFD / CSD, and the maximum error of the elastic deformation does not exceed 0.3%, and the accuracy is relatively high.

[0005] The technical solution adopted by the present invention to solve its technical problems is as follows:

[0006] Step 1: Rapid calculation of the aerodynamic force of a high aspect ratio unmanned aerial vehicle based on the surface panel method;

[0007] Step 2: Check with the numerical simulation method to verify the accuracy of the rapid aerodynamic force calculation method;

[0008] Step 3: Static mechanics calculation of the unmanned aerial vehicle structure based on the classical beam bending-torsion theory;

[0009] Step 4: By constructing an interpolation interface for aerodynamic load-structural displacement, realize the data transfer between the steady aerodynamic force and the structural static mechanics, so as to realize the rapid convergence of the aeroelastic equilibrium state, and apply it to the static aeroelastic rapid calculation of a certain typical unmanned aerial vehicle.

[0010] Preferably, the specific content of Step 1 is as follows:

[0011] Step 1-1: Potential flow theory;

[0012] Given the potential function at any point in the known flow field, the flow field velocity q at that point is obtained:

[0013]

[0014] Among them, Φ represents the potential function, u, v, and w respectively represent the velocity components of the fluid in the x, y, and z directions, and i, j, and k represent the unit vectors in the Cartesian coordinate system, respectively along the x, y, and z axes;

[0015] The continuity equation of the fluid is Among them, ρ is the fluid density;

[0016] When the fluid is incompressible, it is considered that its density remains unchanged and mass is conserved, then The continuity equation is finally transformed into: That is the Laplace equation, which is the continuous state equation of an irrotational and incompressible fluid. Considering the velocity of the object surface, the boundary condition of the object surface can be expressed as: V represents the velocity of the object surface, and n represents the normal vector of the object surface, pointing to the outside of the object surface, which is used to define the flow velocity component in the normal direction in the boundary condition.

[0017] Step 1-2: Three-dimensional grid division;

[0018] Establish a reference coordinate system. Let the oncoming flow direction be the x-axis, the horizontal right direction be the y-axis, and the z-axis be determined by the right-hand screw rule. Divide the mid-arc surface of the wing into a series of columns along the y-axis and into a series of rows along the x-axis. At this time, the lifting surface has been discretized into several grids, that is, vortex grids. Assume that a horseshoe vortex with a constant strength of Γ is placed at any grid.

[0019] The horseshoe vortex consists of three segments of vortex lines. The attached vortex is placed at the 1 / 4 chord line of the grid, and the other two free vortices extend downstream to infinity along the x-axis from the two endpoints of the 1 / 4 chord line along the airflow. Select the midpoint of the 1 / 4 chord line of each vortex grid as the acting point of the force, and the midpoint of the 1 / 4 chord line of the vortex grid as the control point. The boundary condition of the object surface needs to be satisfied at the control point.

[0020] Step 1-3: Biot-Savart law;

[0021] It is known that the distributed vorticity of the finite-length vortex line AB is Γ AB , let Then the induced velocity generated by the vortex segment AB at point C is:

[0022]

[0023] Among them, w AB，C is the aerodynamic influence coefficient of the vortex line AB on point C;

[0024] A constant horseshoe vortex is placed on all vortex grids, and each horseshoe vortex consists of several segments of vortex lines. Therefore, the induced velocity of any horseshoe vortex on any point can be calculated by linear superposition through the Biot-Savart law:

[0025] v i，j = Γ i w i，j (3)

[0026] Among them, w i,j , v i,j are respectively the aerodynamic influence coefficient and the induced velocity of the horseshoe vortex on the i-th vortex grid at the control point of the j-th vortex grid, and Γ i is the vorticity of the horseshoe vortex on the i-th grid; i = 1,..., n, j = 1,..., n;

[0027] Superpose and solve the induced velocity of the horseshoe vortices on the entire wing surface at any vortex lattice control point:

[0028]

[0029] where, v i is the total induced velocity at the i-th vortex lattice control point;

[0030] Steps 1 - 4: Boundary conditions;

[0031] According to the Neumann boundary condition, the induced velocity at the control point should satisfy the boundary condition that the vortex lattice is impenetrable in the normal direction, i.e.:

[0032] (r i，c ×ω + v0 - v i )n i = 0 (5)

[0033] where, r i，c is the vector pointing from the center of gravity to the i-th vortex lattice control point; ω = [ω P ω Q ω R is the aircraft's angular velocity vector, ω P , ω Q , ω R are the roll angular velocity, pitch angular velocity, and sideslip angular velocity of the wing respectively; v i , n i are the induced velocity and unit normal vector at the i-th vortex lattice control point respectively; v0 = V ∞ [cos(α)cos(β) - cos(α)sin(β)sin(α)] is the oncoming air vector, V ∞ is the oncoming air velocity, α is the angle of attack, and β is the sideslip angle;

[0034] The magnitude of the horseshoe vortex vorticity Γ i on all vortex lattices can be solved through the boundary conditions;

[0035] Steps 1 - 5: Aerodynamic force calculation:

[0036] Calculate the induced velocity v i，v at the midpoint of the 1 / 4 chord line segment of the vortex lattice by the Biot - Savart law. In the i-th vortex lattice, the local flow velocity at the attached vortex should be V i，Local = v0 - v i，ν + r i，c ×ω; Use the Kutta - Jukovski theory to solve the aerodynamic force F i acting on the midpoint of the attached vortex:

[0037] F i = ρV i，Local ×Γi L i (6)

[0038] Among them, ρ is the air density, and L i is the vector of the attached vortex segment.

[0039] Preferably, the specific content of step 3 is as follows:

[0040] Step 3-1: Principle of force line translation;

[0041] When translating the force F acting on a rigid body to a point O, in order not to change the effect of the force F on the rigid body, it is necessary to attach a couple. The moment vector of this attached couple is equal to the moment vector of the original force with respect to point O;

[0042] Step 3-2: Bending deformation theory of a hollow circular cross-section conical beam;

[0043] In the case of plane bending, taking the axis before deformation as the x-axis and the axis perpendicular to the axis as the w-axis, with the right direction of the x-axis being positive and the upward direction of the w-axis being positive; the ordinate of any point on the deflection curve is represented by w, and the deflection curve is written as: w = w(x); during the deformation process, the angle θ by which the cross-section rotates relative to its original position around the neutral axis is called the rotation angle of the cross-section; the cross-section rotation angle θ is the angle between the normal of the deflection curve and the w-axis, and the differential relationship between the deflection and the rotation angle is: When bending deformation occurs, the relationship between the curvature and the bending moment is: Among them, ρ(x) is the radius of curvature, M(x) is the bending moment, I(x) is the moment of inertia, and E is Young's modulus;

[0044] The approximate differential equation of the deflection curve of the beam is:

[0045]

[0046] For a hollow circular cross-section conical beam, assuming that the outer diameter and inner diameter of the cross-section change with the length, the outer diameter of the beam is D(x), the inner diameter is d(x), and its moment of inertia is: The stiffness of the hollow circular cross-section conical beam is:

[0047] Step 3-3: Using the numerical integration method to find the deflection curve;

[0048] At discrete points, the trapezoidal rule is used for numerical integration to solve the approximate differential equation of the deflection curve:

[0049]

[0050] where △x = x i+1 -x i ;

[0051] w(x i+1 ) ≈ w(xi ) + θ(x i )△x (9)

[0052] Step 3 - 4: Twisting deformation theory of the hollow circular cross - section conical beam;

[0053] For the hollow circular cross - section conical beam, the second polar moment of inertia J(x) of a certain point on it is: The angle of twist on each differential end dx is: Where T is the resultant torque acting on the beam, and G is the shear modulus of the material. In a homogeneous and isotropic material: μ is the Poisson's ratio; the total angle of twist is obtained by integrating each micro - segment, and the total angle of twist is:

[0054]

[0055] Step 3 - 5: Solving the angle of twist of the beam by the numerical integration method;

[0056] First, discretize the beam. Divide the beam length L into N small segments, and the length of each segment is Use the trapezoidal rule for numerical integration to calculate the angle of twist:

[0057]

[0058] where x i = i△x;

[0059]

[0060] Preferably, the specific content of step 4 is:

[0061] The aerodynamic force calculated by the aerodynamic force module is distributed on many vortex grids and is transformed into the force {F 0} n and moment {M 0} n loaded onto each station in the structural model; the twisting and bending deformations {θ 0} n of the beam are obtained by solving the beam theory, and through the structure - aerodynamic parameter mapping, the bending of the beam is mapped to the change of the dihedral angle of the wing, and the twist of the beam is mapped to the change of the angle of attack, that is, the aerodynamic parameters {α i} n and {Γ i} n of the deformed UAV are obtained. At this time, i = 1, and the iteration starts;

[0062] Update the standard data card of the aerodynamic module, conduct aerodynamic analysis, and solve for the new aerodynamic force Again, load it into the structural model according to the principle of force translation, repeat the above process, and perform iteration until the convergence criterion is met; if it converges, it can be regarded as the stable flight of the high aspect ratio UAV, and the aerodynamic force data after stabilization and the deformation data of the UAV can be obtained.

[0063] A computer program that causes a computer to execute the above-mentioned rapid aeroelastic analysis method.

[0064] An electronic device, including: a processor and a memory; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory so that the electronic device executes the above-mentioned rapid aeroelastic analysis method.

[0065] A computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the above-mentioned rapid aeroelastic analysis method is implemented.

[0066] A chip, including: a processor, which is used to call and run a computer program from a memory, so that a device equipped with the chip executes the above-mentioned rapid aeroelastic analysis method.

[0067] A computer program product, the computer program product includes a computer storage medium, the computer storage medium stores a computer program, the computer program includes instructions that can be executed by at least one processor, and when the instructions are executed by the at least one processor, the above-mentioned rapid aeroelastic analysis method is implemented.

[0068] The beneficial effects of the present invention are as follows:

[0069] (1) By simplifying the geometric model of the high aspect ratio UAV, the skin is discarded, leaving only the main beam and the fuselage structure. Among them, the problem of the force and deformation of the main wing can be regarded as the problem of the force and deformation of a hollow circular cross-section conical composite beam along the span direction. This can facilitate the bending-torsion analysis of the high aspect ratio UAV.

[0070] (2) In terms of calculation accuracy, compared with the calculation results based on the high-precision CFD / CSD numerical method, the maximum error of the elastic deformation of the aeroelastic calculation results of the present invention does not exceed 0.3%, and the accuracy is relatively high.

[0071] (3) In terms of computational efficiency, the CFD / CSD high-precision numerical simulation method requires 13 million grids, and the average time consumption for a single verification case is 50 minutes. However, the fast algorithm of the present invention only requires 500 grids, and the average time consumption for a single verification case is 7.56 seconds, which is only 0.252% of the high-precision numerical method. That is, the computational efficiency is improved by two to three orders of magnitude. Description of the Drawings

[0072] Figure 1 Schematic diagram of the technical route for the fast calculation method of UAV aeroelasticity;

[0073] Figure 2 Schematic diagram of curved surface vortex lattice division;

[0074] Figure 3 Schematic diagram of the horseshoe vortex form;

[0075] Figure 4 Schematic diagram of the vortex line segment;

[0076] Figure 5 Schematic diagram of the grid distribution of a 45° swept wing, (a) three-dimensional view (b) top view;

[0077] Figure 6 Schematic diagram of the relative lift error (CFD vs. curved surface vortex lattice method) for 40 sets of working conditions;

[0078] Figure 7 Schematic diagram of beam deformation;

[0079] Figure 8 Schematic diagram of the structural model data (half model) of a high aspect ratio UAV;

[0080] Figure 9 Deformation diagram of the combined beam under the action of a single force;

[0081] Figure 10 Schematic diagram of a high aspect ratio UAV under a uniformly distributed force;

[0082] Figure 11 Deformation diagram of the combined beam under the action of a uniform force;

[0083] Figure 12 Schematic diagram of the details of the iteration part; ①: Aerodynamic force action point and aerodynamic force data, ②: Transferred force and moment data, ③: Bending deformation rotation angle and torsional deformation twist angle data, ④: Updated file in the form of a standard data card;

[0084] Figure 13 Schematic diagram of the comparison of the displacement of the main beam in Group1 (fast algorithm v.s. FEM);

[0085] Figure 14Schematic diagram of the relative error of the results of the two methods for Group1;

[0086] Figure 15 Schematic diagram of the bending and torsion of a high aspect ratio wing along the rigid axis;

[0087] Figure 16 Schematic diagram of the theorem of translation of force lines. (a) shows the initial force application situation, and (b) shows the force application situation after translating the force F originally located at point A to point B through the theorem of translation of force lines;

[0088] Figure 17 Verification of the bending deformation of a beam - deflection curve under the action of multiple concentrated forces;

[0089] Figure 18 Verification of the torsional deformation of a beam - angle of twist under the action of multiple torques Curve graph;

[0090] Figure 19 Schematic diagram of the relative position between the aerodynamic force application point and the main beam;

[0091] Figure 20 Schematic diagram of the iterative data flow of the fast calculation method;

[0092] Figure 21 Schematic diagram of the change in lift during the iteration process;

[0093] Figure 22 Schematic diagram of the change in tip displacement during the iteration process;

[0094] Figure 23 Schematic diagrams of the bending stress distribution, bending rotation angle, torque, and angle of twist on the main beam. (a) Schematic diagram of the bending stress of the main beam, (b) Schematic diagram of the bending rotation angle of the main beam, (c) Schematic diagram of the torque on the main beam, (d) Schematic diagram of the angle of twist of the main beam;

[0095] Figure 24 Schematic diagram of the lift (coefficient) distribution along the span of a high aspect ratio unmanned aerial vehicle (before deformation vs. after deformation);

[0096] Figure 25 Schematic diagram of the aerodynamic model of a high aspect ratio unmanned aerial vehicle (before deformation vs. after deformation);

[0097] Figure 26 Schematic diagram of the deformation of the main beam. Detailed implementation manner

[0098] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0099] The present invention is a full - aircraft aeroservoelastic deformation and load calculation method for a lightweight flexible solar unmanned aerial vehicle.

[0100] High-altitude long-endurance unmanned aerial vehicles (UAVs) have design characteristics such as large aspect ratios and lightweight materials. Under the conditions of level flight and forced maneuvers, the thin and light wings will produce large bending and torsional deformations under the action of aerodynamic loads, which in turn affect the aerodynamic force distribution of the whole aircraft. The large structural deformations will generate significant aeroelastic coupling stresses, directly affecting the safety assessment of structural strength.

[0101] The present invention aims to develop a computational model and a simplification scheme suitable for rapid static aeroelastic analysis of large-aspect-ratio UAVs in order to accurately analyze aeroelastic deformations and coupling loads under multiple working conditions under the current situation of limited computing resources.

[0102] The following is the main content of the invention:

[0103] a) A rapid method for calculating the aerodynamic forces of large-aspect-ratio UAVs based on the surface panel method;

[0104] b) Checking with numerical simulation methods to verify the accuracy of the rapid method for calculating aerodynamic forces;

[0105] c) A method for calculating the static mechanics of UAV structures based on the classical beam bending-torsion theory;

[0106] d) By constructing an interpolation interface for aerodynamic loads-structural displacements, realizing the data transfer between steady aerodynamic forces and structural statics, thus achieving rapid convergence of the aeroelastic equilibrium state, and applying it to the rapid static aeroelastic calculation of a typical UAV.

[0107] The present invention has developed a rapid method for calculating the static aeroelasticity of large-aspect-ratio UAVs. Figure 1 The technical route of the rapid method for calculating the static aeroelasticity of large-aspect-ratio UAVs is shown. The rapid aeroelastic calculation method can be divided into three parts: preprocessing, rapid aeroelastic calculation method, and postprocessing. The preprocessing part includes measuring the geometric model to construct an aerodynamic model in the form of a standard data card and establishing a simplified beam model of a large-aspect-ratio UAV. As Figure 15 shown in the schematic diagram of the bending and torsion of a large-aspect-ratio wing along the rigid axis.

[0108] The rapid aeroelastic calculation method part includes an aerodynamic force solution module based on the surface panel method and a static mechanics solution module based on the classical beam theory. The aerodynamic loads obtained by solving the aerodynamic force module are loaded into the structural model according to the principle of force translation. The torsion and bending of the beam are solved through the beam theory, and the structural parameters are mapped to aerodynamic parameters. Repeat the above operations until the calculation converges. Finally, the output of the final displacements, aerodynamic force distributions, etc. is the postprocessing part.

[0109] Step 1: Establish a rapid method for calculating the aerodynamic forces of UAVs based on the surface panel method;

[0110] The Vortex Lattice Method is a mature and rapid aerodynamic calculation method, and its core theory was established in the late 1930s. In 1943, Faulkner officially named this method the Vortex Lattice Method (VLM). Since the classical Vortex Lattice Method assumes that the lifting surface is a plane, the simplest wing is a rectangular plane wing. Some aircraft use this plane form, but the layouts of most aircraft are more complex, and it has poor adaptability to curved lifting surfaces. Based on this, Tomas Melin developed the Tornado program, which developed a new curved surface Vortex Lattice Method on the basis of the classical Vortex Lattice Method theory. It has been verified that the aerodynamic calculation results are in good agreement with the experimental data. The calculation efficiency and accuracy of the Vortex Lattice Method make it play an important role in both the conceptual design stage and the iterative design stage of aircraft. In this project, the curved surface Vortex Lattice Method is used for the rapid calculation of the aerodynamics of UAVs.

[0111] Potential flow theory:

[0112] In a flow field, the integral of the velocity vector along any closed curve is called the circulation. If the fluid is irrotational, any two points are arbitrarily selected, and its line integral is not affected by the integration path. At this time, there is a potential function in the flow field. Given the potential function at any point in the flow field, the flow field velocity at that point can be obtained:

[0113]

[0114] The continuity equation of the fluid is where ρ is the fluid density.

[0115] When the fluid is incompressible, its density is considered unchanged and mass is conserved, then The continuity equation finally transforms into: That is the Laplace equation, which is the continuous state equation of an irrotational and incompressible fluid. Considering the object surface has a moving velocity, the object surface boundary condition can be expressed as:

[0116] Three-dimensional grid division:

[0117] Establish a reference coordinate system. Let the oncoming flow direction be the x-axis, the horizontal right direction be the y-axis, and the z-axis be determined by the right-hand screw rule. As Figure 2 shown, the mid-arc surface of the wing is divided into a series of numbers along the y-axis direction and into several rows along the x-axis direction. At this time, the lifting surface has been discretized into several grids (also called vortex lattices). It is assumed that a horseshoe vortex with a constant strength of Γ is placed at any grid.

[0118] As Figure 3It is shown in the form of a horseshoe vortex, which consists of three segments of vortex lines. The attached vortex is placed at the 1 / 4 chord line of the grid, and the other two free vortices extend downstream to infinity along the x-axis from the two endpoints of the 1 / 4 chord line following the airflow. The midpoint of the 1 / 4 chord line of each vortex panel is selected as the acting point of the force, and the midpoint of the 1 / 4 chord line of the vortex panel is the control point, where the boundary condition of the solid surface needs to be satisfied.

[0119] The division of the three-dimensional vortex panel can accurately simulate the lifting surface and provide a basis for calculating the aerodynamic force on the complex aerodynamic surface. By specifying parameters such as the dihedral angle, the twist angle of the wing tip relative to the wing root, the installation angle, the airfoil parameters, and the sweep angle, the complex and diverse wing shapes can be well described, and the three-dimensional vortex panel division is well applicable to such complex shapes.

[0120] Biot-Savart law:

[0121] As Figure 4 shown, given that the vorticity distribution of the finite-length vortex line AB is Γ AB , let then the induced velocity generated by the vortex segment AB at point C is:

[0122]

[0123] where, w AB，C is the aerodynamic influence coefficient of the vortex line AB on point C. A constant horseshoe vortex is placed on all vortex panels, and each horseshoe vortex consists of several segments of vortex lines. Therefore, the induced velocity of any horseshoe vortex on any point can be calculated by linear superposition through the Biot-Savart law:

[0124] v i，j = Γ i w i，j (3)

[0125] where, w i,j , v i,j are respectively the aerodynamic influence coefficient and the induced velocity of the horseshoe vortex on the i-th vortex panel at the control point of the j-th vortex panel, and Γ i is the vorticity of the horseshoe vortex on the i-th grid. The induced velocity of the horseshoe vortices on the entire wing surface at the control point of any vortex panel is solved by superposition:

[0126]

[0127] where, v i is the total induced velocity at the control point of the i-th vortex panel.

[0128] Boundary condition:

[0129] According to the Neumann boundary condition, the induced velocity at the control point should satisfy the boundary condition that the vortex panel is non-penetrable in the normal direction, that is:

[0130] (r i，c ×ω + v0 - v i )n i =0 (5)

[0131] Wherein, r i，c is the vector pointing from the center of gravity to the i-th vortex lattice control point; ω = [ω P ω Q ω R is the aircraft rotational angular velocity vector, ω P , ω Q , ω R are the wing roll angular velocity, pitch angular velocity, and sideslip angular velocity respectively; v i , n i are the induced velocity and unit normal vector at the i-th vortex lattice control point respectively; v0 = V ∞ [cos(α)cos(β) - cos(α)sin(β)sin(α)] is the oncoming air vector, V ∞ is the oncoming air velocity, α is the angle of attack, and β is the sideslip angle; all the horseshoe vortex vorticities Γ i on the vortex lattice can be solved through the boundary conditions.

[0132] Aerodynamic force calculation:

[0133] The induced velocity v i，v at the midpoint of the 1 / 4 chord line segment of the vortex lattice can be calculated by the Biot - Savart law. Therefore, in the i-th vortex lattice, the local flow velocity at the attached vortex should be V i，Local = v0 - v i，ν + r i，c ×ω. The Kutta - Jukovski theory is used to solve the aerodynamic force F i acting on the midpoint of the attached vortex:

[0134] F i = ρV i，Local ×Γ i L i (6)

[0135] Wherein, ρ is the air density, and L i is the vector of the attached vortex line segment. The above is the basic principle of the curved surface vortex lattice method. The flexible arrangement of horseshoe vortices can accurately establish the true aerodynamic force model of the large - deformation wing.

[0136] Step 2: Conduct accuracy verification of the rapid aerodynamic force calculation method based on the numerical simulation method;

[0137] The selected verification example is a 45° swept - wing. The grid division is shown in Figure 5 , and the lift - curve slopes obtained by using the curved surface vortex lattice method of the present invention and other methods respectively The comparison of the calculation results is shown in Table 1. From the data in the table, it can be seen that the calculation results of the surface vortex lattice method are consistent with those obtained by other lifting surface calculation methods. This shows that it is accurate and feasible to calculate the lift line slope of a high aspect ratio unmanned aerial vehicle using the surface vortex lattice method.

[0138] Table 1 Comparison of lift coefficient slopes

[0139]

[0140] In the method provided by the present invention, 40 sets of verification working conditions (divided into 8 groups in total, see Table 1) are set, covering different incoming flow velocities, altitudes and angle of attack parameters (-2°, -1°, 0°, 1°, 2°) to verify the accuracy of the method. The specific parameter configurations are shown in Table 2.

[0141] Table 2 Calculation altitude and incoming flow velocity of verification working conditions

[0142]

[0143]

[0144] In order to verify the accuracy and applicable range of the total lift calculation results of this surface vortex lattice method for a certain high aspect ratio unmanned aerial vehicle under various working conditions, the detailed data of the calculation altitude and incoming flow velocity of the 8 working conditions adopted in this verification are given in Table 2. Figure 6 The relative error diagram of the 40 sets of total lift results obtained by the surface vortex lattice method and CFD simulation calculation is given. From Figure 6 the relative error data in it, it can be seen that the calculation accuracy of this surface vortex lattice method is relatively high, the calculation results are highly consistent with the CFD simulation calculation results, and its average error is only 1.83%. Therefore, various analyses of aerodynamic forces and aerodynamic characteristics of a certain high aspect ratio unmanned aerial vehicle can be carried out through this surface vortex lattice method. Step 3: Establish a calculation method for the structural statics of an unmanned aerial vehicle based on the classical beam bending and torsion theory;

[0145] Theorem of parallel transfer the line of action of a force:

[0146] When translating the force F acting on a rigid body to a certain point O, in order not to change the effect of the force F on the rigid body, a couple must be added. The moment vector of this added couple is equal to the moment vector of the original force with respect to point O. The present invention transfers the aerodynamic force acting on the skin to the main beam based on this theorem.

[0147] Bending deformation theory of a hollow circular cross-section conical beam:

[0148] As Figure 7As shown in the figure, in the case of plane bending, taking the axis before deformation as the x-axis and the axis perpendicular to the axis as the w-axis, with the positive direction of the x-axis to the right and the positive direction of the w-axis upward. The ordinate of any point on the deflection curve is represented by w, so the deflection curve can be written as: w = w(x). During the deformation process, the angle θ through which the cross-section rotates relative to its original position around the neutral axis is called the rotation angle of the cross-section. The cross-section rotation angle θ is the angle between the normal of the deflection curve and the w-axis, and the differential relationship between the deflection and the rotation angle is: The relationship between the curvature and the bending moment during bending deformation is: where ρ(x) is the radius of curvature, M(x) is the bending moment, I(x) is the moment of inertia, and E is Young's modulus.

[0149] The approximate differential equation of the deflection curve of the beam is:

[0150]

[0151] For a tapered beam with a hollow circular cross-section, assuming that the outer diameter and inner diameter of the cross-section vary with the length, the outer diameter of the beam is D(x), the inner diameter is d(x), and its moment of inertia is: The stiffness of the tapered beam with a hollow circular cross-section is:

[0152] Numerical integration method for finding the deflection curve:

[0153] At discrete points, the trapezoidal rule is used in this paper for numerical integration to solve the approximate differential equation of the deflection curve:

[0154]

[0155] where △x = x i+1 -x i

[0156] w(x i+1 )≈w(x i )+θ(x i )△x (9)

[0157] The torsional deformation theory of a tapered beam with a hollow circular cross-section:

[0158] For a tapered beam with a hollow circular cross-section, the polar moment of inertia J(x) of a certain point on it is: The angle of twist on each differential end dx is: where T is the resultant torque acting on the beam, G is the shear modulus of the material, and in a homogeneous and isotropic material: μ is Poisson's ratio. The total angle of twist is obtained by integrating each micro-segment, and the total angle of twist is:

[0159]

[0160] Numerical integration method for solving the torsional angle of a beam:

[0161] First, discretize the beam. Divide the beam length L into N small segments, and the length of each segment is In this paper, the trapezoidal rule is used for numerical integration to calculate the torsional angle:

[0162]

[0163] where x i = iΔx.

[0164]

[0165] Step 4: Check and verify with the deformation of the beam calculated by the finite element method;

[0166] Based on the structural model of a high aspect ratio unmanned aerial vehicle and the theory of a hollow circular cross-section beam, a modeling analysis is carried out for a certain high aspect ratio unmanned aerial vehicle. The specific model and data are as Figure 8 , where for the semi-model analysis of the high aspect ratio unmanned aerial vehicle, the entire beam structure is divided into 3 parts, with lengths L1, L2, and L3 respectively. The first section of the beam is a straight beam with a hollow circular cross-section, with an outer diameter of D0 and an inner diameter of d0. The fuselage is located at the connection between the first section of the beam and the second section of the beam. The second section of the beam is a tapered beam with a hollow circular cross-section, with an outer diameter of D0 and an inner diameter of d0 at the starting point, and an outer diameter of D1 and an inner diameter of d1 at the length L2. The third section of the beam is also a tapered beam with a hollow circular cross-section, which is smoothly connected to the second section of the beam with an included angle of β, that is, the dihedral angle is β, with an outer diameter of D1 and an inner diameter of d1 at the starting point, and an outer diameter of D L , and an inner diameter of d L . As Figure 17 and Figure 18 shown.

[0167] (1) Apply a force F = 500 N to the tip of the beam to verify the accuracy comparison between the self-developed algorithm and the CFD method.

[0168] To verify the accuracy comparison of the bending deformation of the beam calculated by the self-developed algorithm and the CFD method when the beam is under a single acting force, Figure 9 the structural model data (semi-model) of a certain high aspect ratio unmanned aerial vehicle are given in Figure 9 and the deformation diagram of the combined beam under a single force is given in Figure 9 . The blue, red, and green solid lines in Figure 9The blue, red, and green circles in [figure] indicate the beam deformation curves of the three-section beam under ANSYS calculation. From each curve in the figure, it can be seen that the MATLAB data fits well with the ANSYS data, meeting the specific requirements of the problem. Therefore, the deformation results of the composite beam obtained by rapid MATLAB analysis can be applied to the analysis of the beam model of the high-aspect-ratio unmanned aerial vehicle.

[0169] (2) A force F = 50 N is uniformly applied every 5 m at each location of the composite beam, as Figure 10 , to verify the accuracy comparison between the self-developed algorithm and the CFD method.

[0170] To verify the accuracy comparison of the bending deformation of the beam calculated by the self-developed algorithm and the CFD method when the beam is under uniformly distributed force, Figure 10 a schematic diagram of a high-aspect-ratio unmanned aerial vehicle under uniformly distributed force is given in Figure 11 and a deformation diagram of the composite beam under uniformly distributed force is given in Figure 11 The blue, red, and green solid lines in [figure] indicate the beam deformation curves of the three-section beam under the self-developed MATLAB algorithm, Figure 11 The blue, red, and green circles in [figure] indicate the beam deformation curves of the three-section beam under ANSYS calculation. From Figure 11 each curve in [figure], it can be seen that the MATLAB data fits well with the ANSYS data, meeting the specific requirements of the problem. Therefore, the deformation results of the composite beam obtained by rapid MATLAB analysis can be applied to the analysis of the beam model of the high-aspect-ratio unmanned aerial vehicle in the present invention.

[0171] Step Five: Construct an interpolation interface for aerodynamic load - structural displacement to achieve rapid convergence of the aeroelastic equilibrium state;

[0172] As Figure 19 , the aerodynamic force calculated by the aerodynamic force module is distributed on many vortex lattices and is transformed into the force {F 0} n and moment {M 0} n applied to each station in the structural model, as shown in Figure 16 . The beam theory solves to obtain the torsional and bending deformations {θ 0} n and Through the structure - aerodynamic parameter mapping, the bending of the beam is mapped to the change in the dihedral angle of the wing, and the torsion of the beam is mapped to the change in the angle of attack, that is, the aerodynamic parameters {α i} n and {Γ i} n of the deformed unmanned aerial vehicle are obtained. At this time, i = 1, and the iteration begins.

[0173] Update the standard data card of the aerodynamic module, conduct aerodynamic analysis, and solve for the new aerodynamic force Again, according to the principle of force translation, load it into the structural model, repeat the above process, and perform iteration until the convergence criterion is met. If it converges, it can be regarded as the stable flight of the high aspect ratio unmanned aerial vehicle, and the aerodynamic force data after stabilization and various deformation data of the unmanned aerial vehicle can be obtained. The details of the iterative part are as Figure 12 , which shows the data calculation method between the aerodynamic module and the static analysis module

[0174] Step 6: Verification of the fast calculation method for the static aeroelasticity of the unmanned aerial vehicle;

[0175] Compare with the finite element results using the fast algorithm respectively to verify the fast calculation method for the static aeroelasticity of the high aspect ratio unmanned aerial vehicle. The example working conditions are selected from Group 1 in the combination of oncoming flow velocity and altitude. The specific values are shown in Table 3. The angles of attack in each group are taken as -2°, 0°, and 2° respectively. Calculate the static aeroelasticity of the unmanned aerial vehicle using numerical simulation and the fast algorithm respectively

[0176] Table 3 Calculation altitude and oncoming flow velocity corresponding to the verification example working conditions

[0177]

[0178] Working condition Figure 13 and Figure 14 show the comparison diagram of the main beam displacements calculated by numerical simulation and the fast algorithm at three angles of attack of -2°, 0°, and 2° under the height-velocity combination of Group 1 in Table 3, as well as the relative error of the results of the two methods along the spanwise distribution

[0179] Example:

[0180] Figure 20 shows the calculation flow chart of the fast calculation method for the static aeroelasticity of the unmanned aerial vehicle. The vector variables in the figure are marked with "{}"; the superscript is the number of iterations, and 0 means the variable obtained under the rigid shape; the subscript is the vector dimension. For example, for the aerodynamic force distribution calculated by the vortex lattice method, the dimension is the number of vortex lattices, that is, m dimensions

[0181] For the preprocessing part of constructing the aerodynamic model in the form of a standard data card and establishing the simplified beam model of the high aspect ratio unmanned aerial vehicle. Measure the geometric model of the unmanned aerial vehicle, construct the aerodynamic model input in the form of a standard data card, and establish the simplified beam model. The high aspect ratio unmanned aerial vehicle is divided into 100 wing segments. Therefore, as Figure 20 the dimension n of the relevant vector in is 100. Identify the dihedral angle {Γ 0} n and the angle of attack {α 0}n 、airfoil, chord length, etc. Similarly, the simplified beam model is also divided into 100 segments, and the stiffness data {EI} n and {GJ} n are configured on these 100 segments respectively. The load application positions are at 100 stations.

[0182] Static aeroelastic analysis of a certain high aspect ratio unmanned aerial vehicle under the working conditions of an altitude of 13 km, an incoming flow velocity of 17.411 m / s, and an angle of attack of -2°:

[0183] The variation diagram of the lift force during the iteration process is as shown in Figure 21 ; The variation diagram of the tip displacement during the iteration process is as shown in Figure 22 .

[0184] The schematic diagrams of the bending stress distribution, bending angle, torque, and torsional angle on the main beam after stabilization are respectively as shown in Figure 23 (a)-(d) in.

[0185] The diagrams of the lift force, lift coefficient, drag force, and drag coefficient along the span direction of the high aspect ratio unmanned aerial vehicle after stabilization are as shown in Figure 24 .

[0186] The aero model of the high aspect ratio unmanned aerial vehicle is as shown in Figure 25 , showing the comparison of the vortex lattice method aero grids before and after deformation.

[0187] The deformation diagram of the main beam after stabilization is as shown in Figure 26 , and the distribution of the bending stress in the main beam is shown in the form of a contour map.

Claims

1. A rapid analysis method for the static aeroelasticity of a high aspect ratio unmanned aerial vehicle, characterized in that It includes the following steps: Step 1: Rapid aerodynamic calculation of high aspect ratio UAVs based on the surface vortex lattice method; Step 2: Check with the numerical simulation method to verify the accuracy of the rapid aerodynamic calculation method; Step 3: Static structural calculation of UAVs based on the classical beam bending-torsion theory; Step 4: By constructing an interpolation interface for aerodynamic load-structural displacement, realize the data transfer between steady aerodynamic force and static structure, so as to achieve rapid convergence of the aeroelastic equilibrium state, and apply it to the rapid static aeroelastic calculation of a typical UAV.

2. The rapid analysis method of static aeroelasticity for a high aspect ratio unmanned aerial vehicle according to claim 1, wherein The specific content of Step 1 is as follows: Step 1-1: Potential flow theory; Given the potential function at any point in the flow field, obtain the flow field velocity q at that point: Among them, Φ represents the potential function, u, v, and w respectively represent the velocity components of the fluid in the x, y, and z directions, and i, j, and k represent the unit vectors in the Cartesian coordinate system, along the x, y, and z axes respectively; The continuity equation for a fluid is where ρ is the fluid density; When the fluid is incompressible, its density is considered constant, and mass is conserved. Then The continuity equation finally transforms into: That is, the Laplace equation, which is the continuous state equation for an incompressible and irrotational fluid. Considering that the surface has a moving velocity, the surface boundary condition can be expressed as: V represents the surface moving velocity, and n represents the surface normal vector, pointing to the outside of the surface, which is used to define the flow velocity component in the normal direction in the boundary condition. Step 1-2: Three-dimensional grid division; Establish a reference coordinate system, set the oncoming flow direction as the x-axis, the horizontal right as the y-axis, and determine the z-axis by the right-hand screw law; divide the mid-arc surface of the wing into several columns along the y-axis and several rows along the x-axis. At this time, the lifting surface has been discretized into several grids, that is, vortex lattices. Assume that a horseshoe vortex with a constant strength Γ is placed at any grid. The horseshoe vortex consists of three sections of vortex lines. The attached vortex is placed at the 1 / 4 chord line of the grid, and the other two free vortices extend downstream to infinity along the x-axis from the two endpoints of the 1 / 4 chord line along the airflow; select the midpoint of the 1 / 4 chord line of each vortex grid as the action point of the force, and the midpoint of the 1 / 4 chord line of the vortex grid as the control point. The control point needs to satisfy the physical boundary condition; Step 1-3: Biot-Savart law; The known finite-length vortex line AB has a vorticity distribution of Γ AB , let Then the induced velocity generated by the vortex segment AB at point C is:[[]]END]] where w AB，C is the aerodynamic influence coefficient of the vortex line AB on point C; A constant horseshoe vortex is placed on all vortex grids, and each horseshoe vortex consists of several sections of vortex lines. Therefore, the induced velocity of any horseshoe vortex on any point can be calculated by linear superposition through the Biot-Savart law: v i，j = Γ i w i，j (3) where w i,j , v i,j are the aerodynamic influence coefficient and the induced velocity of the horseshoe vortex on the i-th vortex lattice at the control point of the j-th vortex lattice, respectively, and Γ i is the vorticity of the horseshoe vortex on the i-th grid; i = 1, …, n, j = 1, …, n; Superpose and solve the induced velocity of the horseshoe vortices on the entire wing surface at any vortex grid control point: where v i is the total induced velocity at the i-th vortex lattice control point; Step 1-4: Boundary conditions; According to the Neumann boundary condition, the induced velocity at the control point should satisfy the boundary condition that the vortex grid is non-penetrable in the normal direction, that is: (r i，c × ω + v0 - v i )n i = 0 (5) where r i，c is the vector pointing from the center of gravity to the i-th vortex lattice control point; ω = [ω P ω Q ω R is the aircraft angular velocity vector, and ω P , ω Q , ω R are the rolling angular velocity, pitching angular velocity, and sideslip angular velocity of the wing, respectively; v i , n i are the induced velocity and unit normal vector at the i-th vortex lattice control point, respectively; v0 = V ∞ [cos(α)cos(β) - cos(α)sin(β)sin(α)] is the oncoming air vector, V ∞ is the oncoming air velocity, α is the angle of attack, and β is the sideslip angle; The magnitude of the horseshoe vortex vorticity Γ on all vortex cells can be solved through the boundary conditions i ; Step 1-5: Aerodynamic force calculation: Calculate the induced velocity v at the midpoint of the 1 / 4 chord segment of the vortex lattice by the Biot - Savart law i，v , in the i - th vortex lattice, the local flow velocity at the bound vortex should be V i，Local = v0 - v i，ν + r i，c × ω; Use the Kutta - Jukovski theory to solve for the aerodynamic force F acting on the midpoint of the bound vortex i : F i = ρV i，Local × Γ i L i (6) where ρ is the air density, and L i is the vector of the attached vortex segment.

3. A rapid aeroelastic analysis method for a high aspect ratio unmanned aerial vehicle according to claim 2, characterized in that The specific content of Step 3 is as follows: Step 3-1: Principle of force line translation; When translating the force F acting on a rigid body to a certain point O, in order not to change the action effect of the force F on the rigid body, a couple must be added. The moment vector of this added couple is equal to the moment vector of the original force with respect to point O; Step 3-2: Bending deformation theory of a hollow circular cross-section conical beam; In the case of plane bending, taking the axis before deformation as the x-axis and the axis perpendicular to the axis as the w-axis, with the positive direction of the x-axis to the right and the positive direction of the w-axis upward; the ordinate of any point on the deflection curve is represented by w, and the deflection curve is written as: w = w(x); during the deformation process, the angle θ through which the cross-section rotates relative to its original position about the neutral axis is called the rotation angle of the cross-section; the cross-section rotation angle θ is the angle between the normal of the deflection curve and the w-axis, and the differential relationship between the deflection and the rotation angle is: The relationship between the curvature and the bending moment during bending deformation is: where ρ(x) is the radius of curvature, M(x) is the bending moment, I(x) is the moment of inertia, and E is the Young's modulus; The approximate differential equation of the deflection curve of the beam is: For a tapered beam with a hollow circular cross-section, assuming that the outer and inner diameters of the cross-section vary with length, the outer diameter of the beam is D(x) and the inner diameter is d(x), and its moment of inertia is: The stiffness of the tapered beam with a hollow circular cross-section is: Step 3-3: Numerical integration method to find the deflection curve; At discrete points, use the trapezoidal rule for numerical integration to solve the approximate differential equation of the deflection curve: where △x = x i+1 - x i ; w(x i+1 )≈w(x i )+θ(x i )△x(9) Step 3-4: Torsional deformation theory of a hollow circular cross-section conical beam; For a hollow circular cross-section tapered beam, the second moment of area J(x) at a certain point is as follows: The angle of twist on each differential segment dx is: where T is the resultant torque acting on the beam, and G is the shear modulus of the material. In a homogeneous and isotropic material: μ is the Poisson's ratio; the total angle of twist is obtained by integrating over each differential segment, and the total angle of twist is: Step 3-5: Numerical integration method to solve the torsional angle of the beam; First, discretize the beam by dividing the beam length L into N small segments, each with a length of Use the trapezoidal rule for numerical integration to calculate the angle of twist: where x i = iΔx; 4. A rapid static aeroelastic analysis method for a high aspect ratio unmanned aerial vehicle according to claim 3, characterized in that, The specific content of Step 4 is as follows: The aerodynamic force calculated by the aerodynamic force module is distributed on many vortex cells and is transformed into the force {F 0} n and moment {M 0} n acting on the main beam according to the force translation principle, and loaded onto each station in the structural model; the torsional and bending deformations {θ 0} n of the beam are obtained by solving the beam theory, and through the structural-aerodynamic parameter mapping, the bending of the beam is mapped to the change in the dihedral angle of the wing, and the torsion of the beam is mapped to the change in the angle of attack, that is, the aerodynamic parameters {α i} n and {Γ i} n of the deformed UAV are obtained. At this time, i = 1, and the iteration begins; Update the standard data card of the aerodynamic module, perform aerodynamic analysis, and solve for the new aerodynamic force Again, load it into the structural model according to the principle of force translation, repeat the above process, and perform iteration until the convergence criterion is met; if it converges, it can be regarded as the stable flight of the high aspect ratio UAV, and the aerodynamic force data after stabilization and various deformation data of the UAV can be obtained.

5. A computer program, characterized in that, The computer program enables the computer to execute the method described in any one of claims 1 to 4.

6. An electronic device, characterized in that, It includes: A processor and a memory; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory, so that the electronic device executes the method according to any one of claims 1 to 4.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, the method according to any one of claims 1 to 4 is implemented.

8. A chip, characterized in that, Comprising: A processor, configured to call and run a computer program from a memory, so that a device installed with the chip executes the method according to any one of claims 1 to 4.

9. A computer program product, characterized in that, The computer program product includes a computer storage medium storing a computer program, the computer program including instructions executable by at least one processor, and when the instructions are executed by the at least one processor, the method according to any one of claims 1 to 4 is implemented.

Images