A method for efficiently obtaining gust loads of UAVs

Abstract

The present invention relates to a method for efficiently obtaining gust loads of an unmanned aerial vehicle (UAV), belonging to the technical field of aircraft and aeroelasticity. The present invention first determines a nonlinear wing beam model with strain as a variable, and obtains a cross-sectional constitutive relationship matrix of each beam cross section based on the nonlinear wing beam model; secondly, a nonlinear static aeroelastic model is constructed based on an unsteady vortex lattice method, gust disturbances are introduced into the model, and aerodynamic forces corresponding to each beam section are obtained, and then a numerical solution is performed based on the nonlinear beam model to obtain corresponding updated strain quantities; 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 quantities of each beam section and the corresponding cross-sectional constitutive relationship matrix.

Description

Technical Field

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

[0002] High-aspect-ratio flexible aircraft are a highly sought-after research topic in the aeroelastic community both at home and abroad. Due to their large aspect ratio, light weight, high flexibility, and low frequency, they exhibit significant effects such as coupling between rigid-body and elastic motion, large structural deformation, and geometric nonlinearity. These characteristics make the aeroelastic properties of such aircraft more complex and more sensitive to gust disturbances. Therefore, gust load calculation is indispensable for the aeroservoelastic analysis of high-aspect-ratio aircraft.

[0003] For flexible aircraft, gust loads usually refer to the bending moment and torque on the wing section of interest under gust disturbances, both of which are internal forces. Generally, the vibration response of the elastic wing can be obtained by solving the dynamic equation of the wing under airflow disturbances, as shown in the following formula: Here, M, K, C, and F are typically the mass matrix, stiffness matrix, damping matrix, and external load column vector obtained from the structural finite element model. x contains the translational and rotational degrees of freedom of the structural finite element nodes. Because wing profiles are often irregular, integrating the loads within their surface is often difficult. Even if the elastic wing vibration response X is determined at each moment, it is still difficult to directly determine the cross-sectional internal forces and moments.

[0004] In the prior art, the aerodynamic force is linearly integrated along the cross section of interest based on the idea of ​​force balance, and then the aerodynamic force is linearly integrated according to the force balance. The sum of the inertial forces of the outer wing section is obtained by summing the gravity distribution, ignoring the structural damping force. Under these conditions, the internal forces are balanced with the aerodynamic and inertial forces, and the gust loads can be calculated indirectly. However, this calculation method is complex, inefficient, and 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 a numerical solution is performed based on the nonlinear beam model to obtain the corresponding updated strain; 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 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 a UAV, comprising:

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

[0008] Determining a beam reference line of the finite element model of the wing structure along the span direction of the wing; dividing the beam reference line into multiple segments, determining the beam cross-section corresponding to the chord direction of each beam segment of the wing, and obtaining 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] The wing equivalent stiffness of each beam section is obtained based on the strain and internal force in each beam section;

[0011] Obtain the structural generalized mass matrix and generalized damping matrix of the wing based on the finite element model of the wing structure;

[0012] A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix, the generalized damping matrix of the wing structure 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 at the nth beam section, ε y，n is the shear strain in the y-axis direction of the wing at the nth beam section, ε z，n is the shear strain in the z-axis direction of the nth beam section of the wing, 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 nth beam section endpoint i in the x-axis direction, uz，u，n is the displacement of the nth beam section endpoint in the z-axis direction, u z，i，n is the displacement of the nth beam section endpoint i in the z-axis direction, u y，u，n is the displacement of the nth beam section end point j in the y-axis direction, u y，i，n is the displacement of the endpoint 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 endpoint i in the x-axis direction, θ y，j，n is the rotation angle of the nth beam section endpoint in the y-axis direction, θ y，i，n is the rotation angle of the nth beam section endpoint 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 nth beam section endpoint i 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 based on the reference coordinate system B, and the arc length parameter of the beam reference line is given is a real number, and the configuration of the beam reference line at any time includes the beam reference line and the beam cross 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 1

[0023] G1(s, t), basis vector two G2(s, t), basis vector three G3(s, t);

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

[0025] G2(s, t)×G3(s, t); when the wing is not deformed, the basis vector G1(s, 0) is tangent to the center position of the beam reference line, that is, the basis vector G1(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 second strain set 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 in the y-axis direction of the wing 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 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;

[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 variations to characterize discrete wind shear, motion in the aircraft wake, 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 the unsteady vortex lattice method;

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

[0039] A 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 is determined by 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 vortex rings corresponding to a plurality of 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 conditions for satisfying the control points of each vortex grid;

[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 velocity 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 that do not penetrate;

[0050] For example, the midpoint of the spanwise line formed by connecting the 1 / 4 points on both sides of each vortex lattice along the chord direction is selected as the aerodynamic action point; the midpoint of the spanwise line formed by connecting the 3 / 4 points 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 when the wing flies along the flight trajectory and the trailing edge of the wing escapes along the chord direction at the time t, which serves as the tail vortex ring; the amount of the tail vortex ring at the time of escape t is consistent with the amount of the tail vortex ring at the time of escape 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 chord-wise side of the vortex grid of the j-th spanwise row and the i-th chord-wise row as the vortex ring point 1 and vortex ring point 2;

[0057] The quarter of the chord-wise sides of the vortex grid of the j-th spanwise row and the i+1-th chord-wise row are used as vortex ring points 3 and 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 at the initial moment of each vortex grid 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 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] Where Δ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 trailing 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 trailing vortex on the wing, Δc ij is the chord length of the vortex grid in the i-th chordal column and the j-th spanwise column, Δb ij is the span 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 grid in the i-th chord direction and the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid.

[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 updated strains of each beam section;

[0067] Preferably, the specific steps of obtaining the updated strain of each beam 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 aerodynamic force;

[0070] S43. The beam sections are divided into the following sections: m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalence 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 the process returns to step S42. If yes, the calculation is stopped to obtain the updated strain of each beam section;

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

[0073] Step S5: characterize 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 nrepresents 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 a clear conversion relationship with the generalized strain in the simulation, that is, the generalized strain can be calculated from the structural strain, making it easier to verify and compare the simulation results. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0084] Figure 1 Schematic diagram of the wing section load of the beam segment 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 objects, features and advantages of the present invention, the present invention is further described in detail below with reference to 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. Therefore, the scope of protection of the present invention is not limited by the specific embodiments disclosed below.

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

[0090] The present invention provides a method for efficiently obtaining gust loads of a UAV, comprising:

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

[0092] Determining a beam reference line of the finite element model of the wing structure along the span direction of the wing; dividing the beam reference line into multiple segments, determining the beam cross-section corresponding to the chord direction of each beam segment of the wing, and obtaining 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] The wing equivalent stiffness of each beam section is obtained based on the strain and internal force in each beam section;

[0095] Obtain the structural generalized mass matrix and generalized damping matrix of the wing based on the finite element model of the wing structure;

[0096] A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix, the generalized damping matrix of the wing structure 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 at the nth beam section, ε y，n is the shear strain in the y-axis direction of the wing at the nth beam section, ε z，n is the shear strain in the z-axis direction of the nth beam section of the wing, 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 nth beam section endpoint i in the x-axis direction, u z，u，n is the displacement of the nth beam section end point j in the z-axis direction, u z，i，n is the displacement of the nth beam section endpoint i in the z-axis direction, u y，u，nis the displacement of the nth beam section end point j in the y-axis direction, u y，i，n is the displacement of the endpoint 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 endpoint 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 endpoint 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 nth beam section endpoint i 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 based on the reference coordinate system B, and the arc length parameter of the beam reference line is given is a real number, and the configuration of the beam reference line at any time includes the beam reference line and the beam cross 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 1

[0107] G1(s, t), basis vector two G2(s, t), basis vector three G3(s, t);

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

[0109] G2(s, t)×G3(s, t); when the wing is not deformed, the basis vector G1(s, 0) is tangent to the center position of the beam reference line, that is, the basis vector G1(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 second strain set 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 in the y-axis direction of the wing 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 at 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 variations to characterize discrete wind shear, motion in the aircraft wake, 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 the unsteady vortex lattice method;

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

[0123] A 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 is determined by 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 vortex rings corresponding to a plurality of 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 conditions for satisfying the control points of each vortex grid;

[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 velocity 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 that do not penetrate;

[0134] For example, the midpoint of the spanwise line formed by connecting the 1 / 4 points on both sides of each vortex lattice along the chord direction is selected as the aerodynamic action point; the midpoint of the spanwise line formed by connecting the 3 / 4 points 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 when the wing flies along the flight trajectory and the trailing edge of the wing escapes along the chord direction at the time t, which serves as the tail vortex ring; the amount of the tail vortex ring at the time of escape t is consistent with the amount of the tail vortex ring at the time of escape 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 chord-wise side of the vortex grid of the j-th spanwise row and the i-th chord-wise row as the vortex ring point 1 and vortex ring point 2;

[0141] The quarter of the chord-wise sides of the vortex grid of the j-th spanwise row and the i+1-th chord-wise row are used as vortex ring points 3 and 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 at the initial moment of each vortex grid 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 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] Where Δ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 trailing 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 trailing vortex on the wing, Δc ij is the chord length of the vortex grid in the i-th chordal column and the j-th spanwise column, Δb ij is the span 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.

[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 grid in the i-th chord direction and the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid.

[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 updated strains of each beam section;

[0151] Preferably, the specific steps of obtaining the updated strain of each beam 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 aerodynamic force;

[0154] S43. The beam sections are divided into the following sections: m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalence 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 the process returns to step S42. If yes, the calculation is stopped to obtain the updated strain of each beam section;

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

[0157] Step S5: characterize 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 nis 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 embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with this technical field within the technical scope disclosed by the present invention should be covered by the scope of protection of the present invention.

Claims

1. A method for efficiently obtaining gust loads of a UAV, characterized in that: include: Step S1, determining the wing structure finite element model and 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 the 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 updated strains of each beam section; Step S5: characterize 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; Obtaining corresponding wing section loads based on the updated strains of each beam section and the corresponding constitutive relationship matrix; characterizing the corresponding wing section loads as gust loads; 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 is the updated stress of the nth beam section.

2. The method for efficiently obtaining gust loads of a UAV according to claim 1, 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, Determining a beam reference line of the finite element model of the wing structure along the span direction of the wing; dividing the beam reference line into multiple segments, determining the beam cross-section corresponding to the chord direction of each beam segment of the wing, and obtaining 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; The wing equivalent stiffness of each beam section is obtained based on the strain and internal force in each beam section; Obtain the structural generalized mass matrix and generalized damping matrix of the wing based on the finite element model of the wing structure; A nonlinear beam model is established based on the strain in each beam section, the generalized mass matrix, the generalized damping matrix of the wing structure and the wing equivalent stiffness.

3. The method for efficiently obtaining gust loads of a UAV according to claim 2, 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, 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, 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 vortex rings corresponding to a plurality of 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 conditions for satisfying the control points of each vortex grid; Based on the conditions satisfied by each vortex grid control point, the aerodynamic strength of the vortex ring corresponding to each vortex grid is obtained; Obtaining a set of wing motion-related velocities; and obtaining 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 cell 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, 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 chord-wise side of the vortex grid of the j-th spanwise row and the i-th chord-wise row as the vortex ring point 1 and vortex ring point 2; The quarter of the chord-wise sides of the vortex grid of the j-th spanwise row and the i+1-th chord-wise row are used as vortex ring points 3 and 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: Where Δ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 trailing 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 trailing vortex on the wing, Δc ij is the chord length of the vortex grid in the i-th chordal column and the j-th spanwise column, Δb ij is the span 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 loads 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 grid in the i-th chord direction and the j-th span direction; Δp is the aerodynamic pressure at the center of the vortex ring corresponding to the vortex grid.

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 and working condition parameters of each beam section into the nonlinear static aeroelastic model, and introduce gust disturbances; obtain the time t of each beam section m aerodynamic force; S43. The beam sections are divided into the sections at time t m The aerodynamic force is input into the nonlinear beam model, and based on the force equivalence 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.

Images