Method for efficiently obtaining gust load of unmanned aerial vehicle

Abstract

The invention relates to a method for efficiently obtaining the gust load of an unmanned aerial vehicle, and belongs to the technical field of aircrafts and aeroelasticity. The method comprises the steps that firstly, a nonlinear wing beam model with strain as a variable is determined, and a section constitutive relation matrix of each beam section is obtained based on the nonlinear wing beam model; secondly, a nonlinear static aeroelastic model is constructed based on an unsteady vortex grid method, gust disturbance is introduced into the model, aerodynamic force corresponding to each beam section is obtained, and numerical solution is performed based on the nonlinear beam model to obtain corresponding update dependent variables; and finally, determining the wing section of each beam section, and obtaining a corresponding wing section load, namely a gust load, based on the updated dependent variable of each beam section and the corresponding section constitutive relation matrix.

Description

Technical Field

[0001] The present invention relates to the field of aircraft and aeroelastic technology, and in particular to a method for efficiently obtaining gust loads of unmanned aerial vehicles. Background Art

[0002] Flexible aircraft with large aspect ratio is a research object that has attracted much attention in the aeroelastic community at home and abroad. Due to its large aspect ratio, light weight, high flexibility and low frequency, it will show significant rigid body and elastic motion coupling, large structural deformation geometric nonlinearity and other effects. These characteristics will make the aeroelastic characteristics of this type of aircraft more complex and make it more sensitive to gust disturbances. Therefore, gust load calculation is indispensable for the aeroservoelastic analysis of large aspect ratio aircraft.

[0003] For flexible aircraft, gust loads usually refer to the bending moment and torque on the wing section under gust disturbance, both of which are internal forces. In general, people can obtain the vibration response of the elastic wing by solving the dynamic equation of the wing under airflow disturbance, as shown in the following formula: Among them, M, K, C and F are usually the mass matrix, stiffness matrix, damping matrix and external load column vector obtained by the structural finite element model. x contains the translational and rotational degrees of freedom of the structural finite element nodes. Since the wing section is usually irregular, it is usually difficult to integrate the load within its surface. Even if the elastic wing vibration response X at each moment is obtained, it is still difficult to directly obtain the section internal force and internal moment.

[0004] In the prior art, based on the idea of ​​force balance, the aerodynamic force is linearly integrated outward along the cross section of interest, and then the aerodynamic force is linearly integrated outward along the cross section of interest. The sum of the inertial forces of the outer wing section is obtained by summing the gravity distribution, ignoring the structural damping force. In the case of , the internal force is balanced with the aerodynamic force and inertial force, and the gust load can be calculated indirectly. However, this calculation method has complicated solution steps, low efficiency, and is difficult to compare with experiments. Summary of the invention

[0005] In view of the above problems, the present invention provides a method for efficiently obtaining the gust load of a UAV. The present invention first determines a nonlinear wing beam model with strain as a variable, and obtains the cross-sectional constitutive relationship matrix of each beam section based on the nonlinear wing beam model; secondly, a nonlinear static aeroelastic model is constructed based on the unsteady vortex lattice method, gust disturbances are introduced into the model, and the aerodynamic force corresponding to each beam section is obtained, and then the corresponding updated strain amount is obtained by numerical solution based on the nonlinear beam model; finally, the wing section of each beam section is determined, and the corresponding wing section load, i.e., the gust load, is obtained based on the updated strain amount of each beam section and the corresponding cross-sectional constitutive relationship matrix. The present invention is based on the strain beam model calculation, and the calculation method is faster. At the same time, the corresponding load can be obtained by obtaining the strain, which is efficient and direct.

[0006] The present invention provides a method for efficiently obtaining gust loads of an unmanned aerial vehicle, comprising:

[0007] Determine the finite element model of the wing structure,

[0008] Determine the beam reference line of the finite element model of the wing structure along the wing span direction; divide the beam reference line into multiple sections, determine the beam cross section corresponding to the chord direction of each beam section of the wing, and obtain multiple beam cross sections;

[0009] Applying unit loads to the finite element model of the wing structure to obtain the strain and internal force in each beam section; the unit loads include shear force and torque;

[0010] Based on the strain and internal force in each beam section, the wing equivalent stiffness of each beam section is obtained;

[0011] Based on the finite element model of the wing structure, the structural generalized mass matrix and generalized damping matrix of the wing are obtained;

[0012] A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix of the wing structure, the generalized damping matrix and the wing equivalent stiffness.

[0013] Preferably, the expression of the nonlinear beam model in step S1 is:

[0014]

[0015] Where ε is the strain; is the first-order derivative of strain with respect to time, is the second-order derivative of strain with respect to time; M SS is the structural generalized mass matrix; C SS is the generalized damping matrix; K SS is the generalized stiffness matrix; F SS It is an external force, including gravity load and unsteady aerodynamic load (including the force caused by gust disturbance).

[0016] Preferably, the strain of each beam cross section is expressed as:

[0017]

[0018] Among them, ε x，n is the tensile strain of the wing in the x-axis direction of the nth beam section, ε y，n is the shear strain of the wing in the y-axis direction at the nth beam section, ε z，n is the shear strain of the wing in the z-axis direction at the nth beam section, u x，j，n is the displacement of the wing at the end point j of the nth beam section in the x-axis direction, u x，i，n is the displacement of the end point i of the nth beam section in the x-axis direction, uz，u，n is the displacement of the nth beam section end point in the z-axis direction, u z，i，n is the displacement of the end point i of the nth beam section in the z-axis direction, u y，u，n is the displacement of the end point j of the nth beam section in the y-axis direction, u y，i，n is the displacement of the end point i of the nth beam section in the y-axis direction; κ x，n is the torsional strain of the wing in the x-axis direction at the nth beam section, κ y，n is the bending strain of the wing in the y-axis direction at the nth beam section, κ z，n is the bending strain of the wing in the z-axis direction at the nth beam section;

[0019] θ x，j，n is the rotation angle of the nth beam section end point j in the x-axis direction, θ x，i，n is the rotation angle of the nth beam section end point i in the x-axis direction, θ y，j，n is the rotation angle of the nth beam section end point in the y-axis direction, θ y，i，n is the rotation angle of the nth beam section end point i in the y-axis direction, θ z，j，n is the rotation angle of the nth beam section end point j in the z-axis direction, θ z，i，n is the rotation angle of the end point i of the nth beam section in the z-axis direction, l n is the length of the nth beam section, n=1, 2, 3…N, N represents the total number of beam sections.

[0020] Furthermore, the beam reference line is determined according to the reference coordinate system B, and the arc length parameter of the beam reference line is given by is a real number, and the configuration of the beam reference line at any time includes the beam reference line and the beam section;

[0021] The center position vector of the beam reference line is t is the time;

[0022] The basis vectors of the local coordinate system G of the beam section include: basis vector

[0023] G 1 (s, t), basis vector G 2 (s, t), basis vector three G 3 (s, t);

[0024] Basis vector 2G 2 (s, t), basis vector three G 3 (s, t) are two unit vectors orthogonal to each other along the main axis of the cross section at the reference line s of the beam. The basis vector G 1 (s, t) is defined by the right-hand rule, G 1 (s, t) =

[0025] G 2(s, t)×G 3 (s, t); when the wing is not deformed, the basis vector is G 1 (s, 0) is tangent to the center of the beam reference line, i.e., the basis vector G 1 (s, 0) = dR(s, 0) / ds.

[0026] Preferably, the internal force in each beam cross section includes strain set 1 and strain set 2;

[0027] The strain set 1 is a set related to shear stress and shear strain;

[0028] The strain set 2 is a set related to normal stress and normal strain;

[0029] Furthermore, the strain set 1 includes shear force in the y-axis direction, shear force in the z-axis direction, and torque in the x-axis direction; ε y，n is the shear strain of the wing in the y-axis direction at the nth beam section, ε z，n is the shear strain of the wing in the z-axis direction at the nth beam section; κ x，n is the torsional strain of the wing in the x-axis direction at the nth beam section;

[0030] The strain set 2 includes shear force in the x-axis direction, torque in the y-axis direction, and torque in the z-axis direction; ε x，n is the tensile strain of the wing in the x-axis direction of the nth beam section, κ y，n is the bending strain of the wing in the y-axis direction at the nth beam section, κ z，n is the bending strain of the wing in the z-axis direction at the nth beam section;

[0031] Step S2: determining a gust model, and obtaining a gust disturbance based on the gust model;

[0032] Preferably, the specific steps of obtaining the gust disturbance include:

[0033] The incoming flow velocity of the gust at each moment is obtained based on the gust model, and the incoming flow velocity of the gust at each moment is characterized as a gust disturbance.

[0034] Preferably, the gust model is a discrete gust model or a continuous turbulence model;

[0035] The discrete gust model is used for deterministic wind speed changes to characterize discrete wind shear, movement in the aircraft wake zone, or terrain-induced airflow;

[0036] The continuous turbulence model is used to conform to a practical gust disturbance model for characterizing a continuous random process.

[0037] Step S3, constructing a nonlinear static aeroelastic model based on an unsteady vortex lattice method;

[0038] Preferably, the specific steps of constructing the nonlinear static aeroelastic model include:

[0039] The coordinate system is determined based on the cambered surface of the finite element model of the wing; the x-axis of the coordinate system is along the incoming flow direction, i.e., the chord direction, the y-axis is horizontally to the right, i.e., the span direction, and the z-axis depends on the right-hand rule;

[0040] The arc surface of the finite element model of the wing is divided into several quadrilateral grids along the span direction and the chord direction, which are represented as multiple vortex grids.

[0041] Setting a plurality of vortex rings corresponding to the vortex grids;

[0042] Obtaining a vortex ring that is released along the chord direction from the trailing edge of the wing when the wing is in flight, which is defined as a tail vortex ring; and obtaining a gust disturbance caused by the tail vortex ring;

[0043] Determine the control points and aerodynamic action points of each vortex grid; establish the satisfying conditions of each vortex grid control point;

[0044] Based on the satisfying conditions of each vortex grid control point and Biot-Savart theorem, the aerodynamic strength of the vortex ring corresponding to each vortex grid is obtained.

[0045] Get the speed set related to the wing motion;

[0046] The aerodynamic pressure at the center of the vortex ring corresponding to each vortex grid is obtained by using the wing motion related velocity set and the aerodynamic strength of the vortex ring corresponding to each vortex grid;

[0047] Calculate the area of ​​each vortex grid and the corresponding normal vector;

[0048] A nonlinear static aeroelastic model is established based on the area, normal vector and aerodynamic pressure of each vortex cell at the center of the corresponding vortex ring.

[0049] Preferably, the conditions satisfied by each vortex lattice control point are the boundary conditions of the vortex lattice method not penetrating;

[0050] Exemplarily, the midpoint of the spanwise line formed by connecting the 1 / 4 positions on both sides of each vortex lattice along the chord direction is selected as the aerodynamic force action point; the midpoint of the spanwise line formed by connecting the 3 / 4 positions on both sides of each vortex lattice along the chord direction is selected as the control point;

[0051] Characterizing the control point of the vortex grid as the center point of the corresponding vortex ring;

[0052] Preferably, the wing motion related speed set includes the wing motion speed, the gust disturbance speed, and the induced speed caused by the tail vortex ring on the wing;

[0053] Preferably, the tail vortex ring is a vortex ring generated at the chordwise exit time t of the trailing edge of the wing when the wing flies along the flight trajectory, and the ring volume of the tail vortex ring at the exit time t is consistent with that of the tail vortex ring at the exit time t-1.

[0054] Furthermore, the specific steps of obtaining the vortex ring in step S3 include:

[0055] Set the vortex ring corresponding to each vortex grid;

[0056] Select the quarter of the two sides of the chord direction of the vortex grid of the j-th spanwise column and the i-th chordwise column as the vortex ring point 1 and the vortex ring point 2;

[0057] The quarter of the chord-wise sides of the vortex grid of the j-th spanwise column and the i+1-th chord-wise column are taken as the vortex ring point 3 and the vortex ring point 4;

[0058] Connect vortex ring point one, vortex ring point two, vortex ring point three and vortex ring point four in sequence to obtain the corresponding vortex ring.

[0059] Furthermore, when the wing is flying along the flight trajectory, only the vortex rings generated by each vortex grid at the initial moment are located on the wing surface, and the vortex rings generated at other moments are all separated by the trailing edge of the wing as the tail vortex rings;

[0060] Furthermore, the expression of the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid is:

[0061]

[0062] Among them, Δp ij，t is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid of the i-th chordal column and the j-th spanwise column at time t; U(t) is the component of the relative incoming flow velocity of the wing along the x-axis at time t, V(t) is the component of the relative incoming flow velocity of the wing along the y-axis at time t, W(t) is the component of the relative incoming flow velocity of the wing along the z-axis at time t, Γ ij is the vortex ring strength of the vortex grid attached to the i-th chordal column and the j-th spanwise column at time t, u w is the x-axis component of the induced velocity caused by the tail vortex on the wing, v w is the induced velocity component along the y-axis caused by the tail vortex on the wing, w w is the component of the induced velocity along the z-axis caused by the tail vortex on the wing, Δc ij is the chord length of the vortex grid of the i-th chordal column and the j-th spanwise column, Δb ij is the length of the vortex grid of the i-th chordal column and the j-th spanwise column, ρ is the atmospheric density; τ i is the tangent vector of the vortex ring of the i-th chordal column along the chord direction, τ j is the tangent vector of the vortex ring of the jth chordal column along the span direction.

[0063] The expression of the nonlinear static aeroelastic model is:

[0064] ΔF ij = = (ΔpΔs) ij n ij

[0065] Where, ΔF ij is the aerodynamic force of the vortex grid of the i-th chordal column and the j-th spanwise column; Δs is the vortex grid area; n ij is the normal vector of the vortex lattice from the i-th chord to the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex lattice.

[0066] Step S4, updating the strain output by the nonlinear beam model based on the nonlinear static aeroelastic model and the gust disturbance to obtain the updated strain of each beam section;

[0067] Preferably, the specific steps of obtaining the updated strain of each beam cross section in step S4 include:

[0068] S41. When m=1, it indicates the initial moment;

[0069] S42. Input the wing structural state quantity and working condition parameters of each beam section into the wing nonlinear static aeroelastic model, and introduce gust disturbance at the same time; obtain the time t of each beam section m The aerodynamic force;

[0070] S43. The time t of each beam section m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalent condition and the nonlinear time-domain numerical solution method, the updated strain of each beam section is obtained;

[0071] S44. Determine whether m is greater than or equal to M, where M represents the total number of moments. If not, set m+1, t m+1 =t m +Δt, Δt represents the time interval, and returns to step S42. If yes, the calculation is stopped to obtain the updated strain of each beam section;

[0072] Furthermore, the condition for the force equivalence is that the aerodynamic force, mass matrix, damping matrix and stiffness matrix are the same at each moment.

[0073] Step S5, characterizing each beam cross section of the wing as a corresponding wing section;

[0074] The wing equivalent stiffness of each beam section in step S1 is represented as a constitutive relationship matrix of each beam section;

[0075] The corresponding wing section load is obtained based on the updated strain of each beam section and the corresponding constitutive relationship matrix;

[0076] characterizing the corresponding airfoil section load as a gust load;

[0077] Preferably, the wing section load includes section bending moment, shear force and torque.

[0078] Preferably, the gust load expression is:

[0079] F n =K n ε′ n

[0080] Among them, F n represents the gust load of the nth beam section, K n is the constitutive relation matrix of the nth beam section, ε′ n Updated stress in the nth beam section segment.

[0081] Compared with the prior art, the present invention has at least the following beneficial effects:

[0082] The method for obtaining gust loads of the present invention is highly efficient, and the load distribution of the cross section can be directly obtained from the strain, and the load constructed based on the strain is easier; the present invention is used in gust tests to more conveniently measure the strain of the structure to be tested, and has an obvious conversion relationship with the generalized strain in the simulation, that is, the generalized strain can be calculated through the structural strain, and the simulation results are more easily verified and compared. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] The drawings are only for the purpose of illustrating particular embodiments and are not to be construed as limiting the invention.

[0084] Figure 1 Schematic diagram of the wing section load of the beam section in an embodiment of the present invention;

[0085] Figure 2 A flowchart of obtaining an update strain in an embodiment of the present invention;

[0086] Figure 3 A schematic diagram of aerodynamic grid division in an embodiment of the present invention;

[0087] Figure 4 Schematic diagram of the coordinates of the beam reference line in an embodiment of the present invention DETAILED DESCRIPTION

[0088] In order to more clearly understand the above-mentioned purpose, features and advantages of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments can be combined with each other. In addition, the present invention can also be implemented in other ways different from those described herein, and therefore, the protection scope of the present invention is not limited by the specific embodiments disclosed below.

[0089] A specific embodiment of the present invention, as Figure 1-4The present invention discloses a method for efficiently obtaining the gust load of a UAV. In order to illustrate the effectiveness of the method proposed by the present invention, the above technical solution of the present invention is described in detail through a specific embodiment below. The specific implementation steps are as follows:

[0090] The present invention provides a method for efficiently obtaining gust loads of an unmanned aerial vehicle, comprising:

[0091] Determine the finite element model of the wing structure,

[0092] Determine the beam reference line of the finite element model of the wing structure along the wing span direction; divide the beam reference line into multiple sections, determine the beam cross section corresponding to the chord direction of each beam section of the wing, and obtain multiple beam cross sections;

[0093] Applying unit loads to the finite element model of the wing structure to obtain the strain and internal force in each beam section; the unit loads include shear force and torque;

[0094] Based on the strain and internal force in each beam section, the wing equivalent stiffness of each beam section is obtained;

[0095] Based on the finite element model of the wing structure, the structural generalized mass matrix and generalized damping matrix of the wing are obtained;

[0096] A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix of the wing structure, the generalized damping matrix and the wing equivalent stiffness.

[0097] Preferably, the expression of the nonlinear beam model in step S1 is:

[0098]

[0099] Where ε is the strain; is the first-order derivative of strain with respect to time, is the second-order derivative of strain with respect to time; M SS is the structural generalized mass matrix; C SS is the generalized damping matrix; K SS is the generalized stiffness matrix; F SS It is an external force, including gravity load and unsteady aerodynamic load (including the force caused by gust disturbance).

[0100] Preferably, the strain of each beam cross section is expressed as:

[0101]

[0102] Among them, ε x，n is the tensile strain of the wing in the x-axis direction of the nth beam section, ε y，n is the shear strain of the wing in the y-axis direction at the nth beam section, ε z，nis the shear strain of the wing in the z-axis direction at the nth beam section, u x，j，n is the displacement of the wing at the end point j of the nth beam section in the x-axis direction, u x，i，n is the displacement of the end point i of the nth beam section in the x-axis direction, u z，u，n is the displacement of the end point j of the nth beam section in the z-axis direction, u z，i，n is the displacement of the end point i of the nth beam section in the z-axis direction, u y，u，n is the displacement of the end point j of the nth beam section in the y-axis direction, u y，i，n is the displacement of the end point i of the nth beam section in the y-axis direction; κ x，n is the torsional strain of the wing in the x-axis direction at the nth beam section, κ y，n is the bending strain of the wing in the y-axis direction at the nth beam section, κ z，n is the bending strain of the wing in the z-axis direction at the nth beam section;

[0103] θ x，j，n is the rotation angle of the nth beam section end point j in the x-axis direction, θ x，i，n is the rotation angle of the nth beam section end point i in the x-axis direction, θ y，j，n is the rotation angle of the nth beam section end point j in the y-axis direction, θ y,i，n is the rotation angle of the nth beam section end point i in the y-axis direction, θ z,j，n is the rotation angle of the nth beam section end point j in the z-axis direction, θ z，i，n is the rotation angle of the end point i of the nth beam section in the z-axis direction, l n is the length of the nth beam section, n=1, 2, 3…N, N represents the total number of beam sections.

[0104] Furthermore, the beam reference line is determined according to the reference coordinate system B, and the arc length parameter of the beam reference line is given by is a real number, and the configuration of the beam reference line at any time includes the beam reference line and the beam section;

[0105] The center position vector of the beam reference line is t is the time;

[0106] The basis vectors of the local coordinate system G of the beam section include: basis vector

[0107] G 1 (s, t), basis vector G 2 (s, t), basis vector three G 3 (s, t);

[0108] Basis vector 2G 2 (s, t), basis vector three G 3(s, t) are two unit vectors orthogonal to each other along the main axis of the cross section at the reference line s of the beam. The basis vector G 1 (s, t) is defined by the right-hand rule, G 1 (s, t) =

[0109] G 2 (s, t)×G 3 (s, t); when the wing is not deformed, the basis vector is G 1 (s, 0) is tangent to the center of the beam reference line, i.e., the basis vector G 1 (s, 0) = dR(s, 0) / ds.

[0110] Preferably, the internal force in each beam cross section includes strain set 1 and strain set 2;

[0111] The strain set 1 is a set related to shear stress and shear strain;

[0112] The strain set 2 is a set related to normal stress and normal strain;

[0113] Furthermore, the strain set 1 includes shear force in the y-axis direction, shear force in the z-axis direction, and torque in the x-axis direction; ε y，n is the shear strain of the wing in the y-axis direction at the nth beam section, ε z，n is the shear strain of the wing in the z-axis direction at the nth beam section; κ x，n is the torsional strain of the wing in the x-axis direction at the nth beam section;

[0114] The strain set 2 includes shear force in the x-axis direction, torque in the y-axis direction, and torque in the z-axis direction; ε x，n is the tensile strain of the wing in the x-axis direction of the nth beam section, K y，n is the bending strain of the wing in the y-axis direction at the nth beam section, κ z，n is the bending strain of the wing in the z-axis direction at the nth beam section;

[0115] Step S2: determining a gust model, and obtaining a gust disturbance based on the gust model;

[0116] Preferably, the specific steps of obtaining the gust disturbance include:

[0117] The incoming flow velocity of the gust at each moment is obtained based on the gust model, and the incoming flow velocity of the gust at each moment is characterized as a gust disturbance.

[0118] Preferably, the gust model is a discrete gust model or a continuous turbulence model;

[0119] The discrete gust model is used for deterministic wind speed changes to characterize discrete wind shear, movement in the aircraft wake zone, or terrain-induced airflow;

[0120] The continuous turbulence model is used to conform to a practical gust disturbance model for characterizing a continuous random process.

[0121] Step S3, constructing a nonlinear static aeroelastic model based on an unsteady vortex lattice method;

[0122] Preferably, the specific steps of constructing the nonlinear static aeroelastic model include:

[0123] The coordinate system is determined based on the cambered surface of the finite element model of the wing; the x-axis of the coordinate system is along the incoming flow direction, i.e., the chord direction, the y-axis is horizontally to the right, i.e., the span direction, and the z-axis depends on the right-hand rule;

[0124] The arc surface of the finite element model of the wing is divided into several quadrilateral grids along the span direction and the chord direction, which are represented as multiple vortex grids.

[0125] Setting a plurality of vortex rings corresponding to the vortex grids;

[0126] Obtaining a vortex ring that is released along the chord direction from the trailing edge of the wing when the wing is in flight, which is defined as a tail vortex ring; and obtaining a gust disturbance caused by the tail vortex ring;

[0127] Determine the control points and aerodynamic action points of each vortex grid; establish the satisfying conditions of each vortex grid control point;

[0128] Based on the satisfying conditions of each vortex grid control point and Biot-Savart theorem, the aerodynamic strength of the vortex ring corresponding to each vortex grid is obtained.

[0129] Get the speed set related to the wing motion;

[0130] The aerodynamic pressure at the center of the vortex ring corresponding to each vortex grid is obtained by using the wing motion related velocity set and the aerodynamic strength of the vortex ring corresponding to each vortex grid;

[0131] Calculate the area of ​​each vortex grid and the corresponding normal vector;

[0132] A nonlinear static aeroelastic model is established based on the area, normal vector and aerodynamic pressure of each vortex cell at the center of the corresponding vortex ring.

[0133] Preferably, the conditions satisfied by each vortex lattice control point are the boundary conditions of the vortex lattice method not penetrating;

[0134] Exemplarily, the midpoint of the spanwise line formed by connecting the 1 / 4 positions on both sides of each vortex lattice along the chord direction is selected as the aerodynamic force action point; the midpoint of the spanwise line formed by connecting the 3 / 4 positions on both sides of each vortex lattice along the chord direction is selected as the control point;

[0135] Characterizing the control point of the vortex grid as the center point of the corresponding vortex ring;

[0136] Preferably, the wing motion related speed set includes the wing motion speed, the gust disturbance speed, and the induced speed caused by the tail vortex ring on the wing;

[0137] Preferably, the tail vortex ring is a vortex ring generated at the chordwise exit time t of the trailing edge of the wing when the wing flies along the flight trajectory, and the ring volume of the tail vortex ring at the exit time t is consistent with that of the tail vortex ring at the exit time t-1.

[0138] Furthermore, the specific steps of obtaining the vortex ring in step S3 include:

[0139] Set the vortex ring corresponding to each vortex grid;

[0140] Select the quarter of the two sides of the chord direction of the vortex grid of the j-th spanwise column and the i-th chordwise column as the vortex ring point 1 and the vortex ring point 2;

[0141] The quarter of the chord-wise sides of the vortex grid of the j-th spanwise column and the i+1-th chord-wise column are taken as the vortex ring point 3 and the vortex ring point 4;

[0142] Connect vortex ring point one, vortex ring point two, vortex ring point three and vortex ring point four in sequence to obtain the corresponding vortex ring.

[0143] Furthermore, when the wing is flying along the flight trajectory, only the vortex rings generated by each vortex grid at the initial moment are located on the wing surface, and the vortex rings generated at other moments are all separated by the trailing edge of the wing as the tail vortex rings;

[0144] Furthermore, the expression of the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid is:

[0145]

[0146] Among them, Δp ij，t is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid of the i-th chordal column and the j-th spanwise column at time t; U(t) is the component of the relative incoming flow velocity of the wing along the x-axis at time t, V(t) is the component of the relative incoming flow velocity of the wing along the y-axis at time t, W(t) is the component of the relative incoming flow velocity of the wing along the z-axis at time t, Γ ij is the vortex ring strength of the vortex grid attached to the i-th chordal column and the j-th spanwise column at time t, u w is the x-axis component of the induced velocity caused by the tail vortex on the wing, v w is the induced velocity component along the y-axis caused by the tail vortex on the wing, w w is the component of the induced velocity along the z-axis caused by the tail vortex on the wing, Δc ij is the chord length of the vortex grid of the i-th chordal column and the j-th spanwise column, Δb ij is the length of the vortex grid of the i-th chordal column and the j-th spanwise column, ρ is the atmospheric density; τ i is the tangent vector of the vortex ring of the i-th chordal column along the chord direction, τ jis the tangent vector of the vortex ring of the jth chordal column along the span direction.

[0147] The expression of the nonlinear static aeroelastic model is:

[0148] ΔF ij =-(ΔpΔs) ij n ij

[0149] Where, ΔF ij is the aerodynamic force of the vortex grid of the i-th chordal column and the j-th spanwise column; Δs is the vortex grid area; n ij is the normal vector of the vortex lattice from the i-th chord to the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex lattice.

[0150] Step S4, updating the strain output by the nonlinear beam model based on the nonlinear static aeroelastic model and the gust disturbance to obtain the updated strain of each beam section;

[0151] Preferably, the specific steps of obtaining the updated strain of each beam cross section in step S4 include:

[0152] S41. When m=1, it indicates the initial moment;

[0153] S42. Input the wing structural state quantity and working condition parameters of each beam section into the wing nonlinear static aeroelastic model, and introduce gust disturbance at the same time; obtain the time t of each beam section m The aerodynamic force;

[0154] S43. The time t of each beam section m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalent condition and the nonlinear time-domain numerical solution method, the updated strain of each beam section is obtained;

[0155] S44. Determine whether m is greater than or equal to M, where M represents the total number of moments. If not, set m+1, t m+1 =t m +Δt, Δt represents the time interval, and returns to step S42. If yes, the calculation is stopped to obtain the updated strain of each beam section;

[0156] Furthermore, the condition for the force equivalence is that the aerodynamic force, mass matrix, damping matrix and stiffness matrix are the same at each moment.

[0157] Step S5, characterizing each beam cross section of the wing as a corresponding wing section;

[0158] The wing equivalent stiffness of each beam section in step S1 is represented as a constitutive relationship matrix of each beam section;

[0159] The corresponding wing section load is obtained based on the updated strain of each beam section and the corresponding constitutive relationship matrix;

[0160] characterizing the corresponding airfoil section load as a gust load;

[0161] Preferably, the wing section load includes section bending moment, shear force and torque.

[0162] Preferably, the gust load expression is:

[0163] F n =K n ε′ n

[0164] Among them, F n represents the gust load of the nth beam section, K n is the constitutive relation matrix of the nth beam section, ε′ n Updated stress in the nth beam section segment.

[0165] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by any technician familiar with the technical field within the technical scope disclosed by the present invention should be covered within the protection scope of the present invention.

Claims

1. A method for efficiently obtaining gust loads of unmanned aerial vehicles, characterized in that: include: Step S1, determining a wing structure finite element model and a beam reference line; Obtaining a plurality of beam cross sections based on the beam reference line; Applying a unit load to the finite element model of the wing structure to obtain the strain in each beam section; obtaining the wing equivalent stiffness of each beam section based on the strain in each beam section; Establishing a nonlinear beam model based on the strain in each beam section and the wing equivalent stiffness; Step S2: determining a gust model, and obtaining a gust disturbance based on the gust model; Step S3, constructing a nonlinear static aeroelastic model based on an unsteady vortex lattice method; Step S4, updating the strain output by the nonlinear beam model based on the nonlinear static aeroelastic model and the gust disturbance to obtain the updated strain of each beam section; Step S5, characterizing each beam cross section of the wing as a corresponding wing section; The wing equivalent stiffness of each beam section in step S1 is represented as a constitutive relationship matrix of each beam section; Based on the updated strain of each beam section and the corresponding constitutive relationship matrix, a corresponding wing section load is obtained; and the corresponding wing section load is characterized as a gust load.

2. The method for efficiently obtaining gust load of a UAV according to claim 1 is characterized in that: The specific steps of establishing the nonlinear beam model in step S1 include: Determine the finite element model of the wing structure, Determine the beam reference line of the finite element model of the wing structure along the wing span direction; divide the beam reference line into multiple sections, determine the beam cross section corresponding to the chord direction of each beam section of the wing, and obtain multiple beam cross sections; Apply unit load to the finite element model of the wing structure to obtain the strain and internal force in each beam section; Based on the strain and internal force in each beam section, the wing equivalent stiffness of each beam section is obtained; Based on the finite element model of the wing structure, the structural generalized mass matrix and generalized damping matrix of the wing are obtained; A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix of the wing structure, the generalized damping matrix and the wing equivalent stiffness.

3. The method for efficiently obtaining gust load of a UAV according to claim 2 is characterized in that: The expression of the nonlinear beam model in step S1 is: Where ε is the strain; is the first-order derivative of strain with respect to time, is the second-order derivative of strain with respect to time; M SS is the structural generalized mass matrix; C SS is the generalized damping matrix; K SS is the generalized stiffness matrix; F SS For external force.

4. The method for efficiently obtaining gust loads of a UAV according to claim 1 is characterized in that: The gust model is a discrete gust model or a continuous turbulence model.

5. The method for efficiently obtaining gust loads of a UAV according to claim 1 is characterized in that: The specific steps of constructing the nonlinear static aeroelastic model in step S3 include: Determine the coordinate system based on the arc surface of the wing finite element model; The arc surface of the finite element model of the wing is divided into several quadrilateral grids along the span direction and the chord direction, which are represented as multiple vortex grids. Setting a plurality of vortex rings corresponding to the vortex grids; Obtaining a vortex ring that is released along the chord direction from the trailing edge of the wing when the wing is in flight, which is defined as a tail vortex ring; and obtaining a gust disturbance caused by the tail vortex ring; Determine the control points and aerodynamic action points of each vortex grid; establish the satisfying conditions of each vortex grid control point; Based on the satisfied conditions of each vortex grid control point, the aerodynamic strength of the vortex ring corresponding to each vortex grid is obtained; Acquire a set of wing motion-related velocities; and obtain the aerodynamic pressure at the center of the vortex ring corresponding to each vortex grid using the set of wing motion-related velocities and the aerodynamic strength of the vortex ring corresponding to each vortex grid; Calculate the area of ​​each vortex grid and the corresponding normal vector; A nonlinear static aeroelastic model is established based on the area, normal vector of each vortex grid and the aerodynamic pressure at the center of the corresponding vortex ring.

6. The method for efficiently obtaining gust loads of a UAV according to claim 5 is characterized in that: The specific steps of obtaining the vortex ring in step S3 include: Set the vortex ring corresponding to each vortex grid; Select the quarter of the two sides of the chord direction of the vortex grid of the j-th spanwise column and the i-th chordwise column as the vortex ring point 1 and the vortex ring point 2; The quarter of the chord-wise sides of the vortex grid of the j-th spanwise column and the i+1-th chord-wise column are taken as the vortex ring point 3 and the vortex ring point 4; Connect vortex ring point one, vortex ring point two, vortex ring point three and vortex ring point four in sequence to obtain the corresponding vortex ring.

7. The method for efficiently obtaining gust loads of a UAV according to claim 5, characterized in that: The expression of the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid is: Among them, Δp ij,t is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid of the i-th chordal column and the j-th spanwise column at time t; U(t) is the component of the relative incoming flow velocity of the wing along the x-axis at time t, V(t) is the component of the relative incoming flow velocity of the wing along the y-axis at time t, W(t) is the component of the relative incoming flow velocity of the wing along the z-axis at time t, Γ ij is the vortex ring strength of the vortex grid attached to the i-th chordal column and the j-th spanwise column at time t, u w is the x-axis component of the induced velocity caused by the tail vortex on the wing, v w is the induced velocity component along the y-axis caused by the tail vortex on the wing, w w is the component of the induced velocity along the z-axis caused by the tail vortex on the wing, Δc ij is the chord length of the vortex grid of the i-th chordal column and the j-th spanwise column, Δb ij is the length of the vortex grid of the i-th chordal column and the j-th spanwise column, ρ is the atmospheric density; τ i is the tangent vector of the vortex ring of the i-th chordal column along the chord direction, τ j is the tangent vector of the vortex ring of the jth chordal column along the span direction.

8. The method for efficiently obtaining gust load of a UAV according to claim 7, characterized in that: The expression of the nonlinear static aeroelastic model is: ΔF ij =-(ΔpΔs) ij n ij Where, ΔF ij is the aerodynamic force of the vortex grid of the i-th chordal column and the j-th spanwise column; Δs is the vortex grid area; n ij is the normal vector of the vortex lattice from the i-th chord to the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex lattice.

9. The method for efficiently obtaining gust loads of a UAV according to claim 8, characterized in that: The specific steps of obtaining the updated strain of each beam section in step S4 include: S41. When m=1, it indicates the initial moment; S42. Input the wing structure state quantity and working condition parameters of each beam section into the nonlinear static aeroelastic model, and introduce gust disturbance at the same time; obtain the time t of each beam section m The aerodynamic force; S43. The time t of each beam section m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalent condition and the nonlinear time-domain numerical solution method, the updated strain of each beam section is obtained; S44. Determine whether m is greater than or equal to M, where M represents the total number of moments. If not, set m+1, t m+1 =t m +Δt, Δt represents the time interval, and the process returns to step S42. If yes, the calculation is stopped to obtain the updated strain of each beam section.

10. The method for efficiently obtaining gust loads of a UAV according to claim 1, characterized in that: The gust load expression is: F n =K n e' n Among them, F n represents the gust load of the nth beam section, K n is the constitutive relation matrix of the nth beam section, ε′ n Updated stress in the nth beam section segment.

Images