CN116956782B - Nonlinear flutter analysis method - Google Patents

Nonlinear flutter analysis method Download PDF

Info

Publication number
CN116956782B
CN116956782B CN202311216406.9A CN202311216406A CN116956782B CN 116956782 B CN116956782 B CN 116956782B CN 202311216406 A CN202311216406 A CN 202311216406A CN 116956782 B CN116956782 B CN 116956782B
Authority
CN
China
Prior art keywords
wing
matrix
aerodynamic
object plane
flutter
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202311216406.9A
Other languages
Chinese (zh)
Other versions
CN116956782A (en
Inventor
马宇航
史亚云
王波
周佳星
贺龙
邹毅
陈耿
杨体浩
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202311216406.9A priority Critical patent/CN116956782B/en
Publication of CN116956782A publication Critical patent/CN116956782A/en
Application granted granted Critical
Publication of CN116956782B publication Critical patent/CN116956782B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a nonlinear flutter analysis method, which belongs to the technical field of wing structural design, and comprises the following steps of: s1: acquiring an initial value of a flight parameter; s2: according to the current flight parameters, analyzing the static pneumatic elasticity of the wing structure under large deformation to obtain a static pneumatic elasticity analysis result; s3: carrying out flutter analysis according to the static pneumatic elasticity analysis result to obtain flutter speed and frequency; s4: judging whether the flutter speed is matched with the incoming flow speed, if so, entering S5, otherwise, entering S6; s5: outputting the flutter speed and the frequency as nonlinear flutter analysis results; s6: and (2) adjusting the flight parameters according to the flutter speed and returning to the step (S2).

Description

Nonlinear flutter analysis method
Technical Field
The invention relates to the technical field of wing structural design, in particular to a nonlinear flutter analysis method.
Background
In the flight of an aircraft, aerodynamic loads borne by the wing can cause the wing to deform greatly, the geometrical nonlinear effect is obvious, and the structural rigidity of the wing can be changed.
In aircraft structural analysis, the nonlinear problem mainly includes: material, geometry and boundary conditions are nonlinear. For boundary nonlinearity, the boundary condition is the main source of nonlinear effects with the motion of the object; the nonlinear effect of material nonlinearity mainly comes from nonlinearity of stress-strain relation, and constitutive equation is different from linear problem; the geometrical nonlinearity is different from the geometrical equation, and the stress-strain relation is changed due to the occurrence of large strain, so that the corresponding balance, principal structure and virtual work equation are required to be redefined. The geometrical nonlinearity mainly includes: large displacement small strain, large displacement large strain, change in magnitude, direction, or boundary support conditions of external load due to structural deformation, and the like.
Disclosure of Invention
The invention aims to provide a nonlinear flutter analysis method for improving applicability and reliability in the wing design process.
The technical scheme for solving the technical problems is as follows:
the invention provides a nonlinear flutter analysis method, which comprises the following steps:
s1: acquiring an initial value of a flight parameter;
s2: according to the current flight parameters, analyzing the static pneumatic elasticity of the wing structure under large deformation to obtain a static pneumatic elasticity analysis result;
s3: carrying out flutter analysis according to the static pneumatic elasticity analysis result to obtain flutter speed and frequency;
s4: judging whether the flutter speed is matched with the incoming flow speed, if so, entering S5, otherwise, entering S6;
s5: outputting the flutter speed and the frequency as nonlinear flutter analysis results;
s6: and (2) adjusting the flight parameters according to the flutter speed and returning to the step (S2).
Optionally, the static aeroelastic analysis result includes a structural node nonlinear displacement, and the S2 includes:
s21: solving aerodynamic force by using a surface element method/RANS (random access network) consideration transition method according to the current flight parameters;
s22: interpolating the aerodynamic force into target structural points of the wing structure to obtain the force of each structural point;
s23: obtaining the displacement of each current structural point according to the relation between the stress condition of the force of each structural point and the rigidity matrix;
s24: and judging whether the deformation of the current wing structure is converged or not according to the displacement of each current structural point, if so, outputting the displacement of each current structural point as the nonlinear displacement of the structural node, otherwise, updating the pneumatic surface and returning to S21.
Optionally, in S21, according to the current flight parameter, solving aerodynamic force by using a bin method includes:
according toControl equation of time unsteady surface element method, determining action on the firstjPressure coefficients on the individual object plane elements;
the saidThe control equation of the time unsteady surface element method is as follows:
wherein the subscriptExpress object plane->Indicates wake up>Is->All panels at time->Surface of each objectMatrix of aerodynamic influence coefficients of object plane dipole intensity on element,/->Is->All panels are at the first timeA matrix of wake dipole intensity aerodynamic influence coefficients on individual object plane elements +.>Is->Moment object plane source strong aerodynamic force influence coefficient matrix; />Is->Dipole intensity of the object plane element at moment +.>Is thatDipole intensity of the moment trail surface element; />Is->The object plane source at the moment is strong，/>Is->A moment object plane surface element vector matrix;/>、/>and->Respectively->The moment is infinite free incoming flow speed, the movement speed of an object plane due to structural vibration and the disturbance speed of gusts; subscript->Represents free incoming flow at infinity,>representing the vibration of the object plane structure->Indicating gusts;
wherein:is->Pressure coefficients at the individual object plane element control points; />And->Acting at +.>Static pressure at the control point of each object plane element and static pressure of free incoming flow at infinity; />Is the atmospheric density; />Is->Local fluid velocities at the individual object plane bin control points; />Is->Velocity profile at the control point of the individual object plane element,/->Time is;
determining the pneumatic load of the surface element according to the pressure coefficient and the current flight parameter;
pneumatic loading of the doughThe method comprises the following steps:
wherein,for reference area->Is->The surface element areas of the individual object surfaces;
and obtaining aerodynamic force and aerodynamic moment according to the pressure coefficient and the aerodynamic load of the surface element.
Optionally, in S21, according to the current flight parameter, solving aerodynamic force by using the RANS to consider a transition method includes:
solving a flow field conservation variable according to a three-dimensional compressible unsteady N-S equation and an intermittent factor transportation equation;
the three-dimensional compressible unsteady N-S equation is:
wherein,is the coordinate of any point in the flow field, +.>Is a flow field conservation variable +.>、/>And->Respectively->Non-stick flux in three directions, +.>，/>And->Respectively->Three-directional adhesive flux, +.>For air density->Speed is +.>Components in three directions>Total internal energy per unit volume; />Is the sign of the partial derivative;
the intermittent factorThe transport equation in conservation form is:
wherein,for density (I)>For time (I)>To generate items->For dissipating items->Molecular viscosity>For vortex viscosity>Is constant and->=1.0；/>The values are 1, 2, 3, < >>Is->And->Coordinates representing three directions, +.>Is->And->Is->Speeds in three directions; />Is the sign of the partial derivative;
determining a pressure coefficient and a lift coefficient according to the flow field conservation variable;
the pressure coefficientAnd lift coefficient->The method comprises the following steps of:
wherein,is surface static pressure>Is free incoming flow static pressure, +.>Is the free-flowing dynamic pressure of the incoming flow,is the lower surface pressure coefficient,/->Is the upper surface pressure coefficient,/-, is->Is the normalized chordwise location;
coefficient of surface frictionIs surface shear stress->Ratio to free incoming flow pressure:
and obtaining aerodynamic force according to the pressure coefficient, the lift coefficient and the surface friction coefficient.
Aerodynamic forces are obtained by multiplying the coefficients by dynamic pressure and reference area.
Alternatively, the aerodynamic force is interpolated into target structural points of the wing structure using radial basis function interpolation to obtain the forces of the respective structural points.
Optionally, in S24, the aerodynamic force is recalculated by substituting the displacement of each current structural point into the geometric deformation to obtain a new displacement, and whether the difference between the new displacement and the displacement of each current structural point is within a preset threshold is determined, if yes, the deformation of the current wing structure converges, otherwise, the deformation does not converge.
Optionally, the static pneumatic elastic analysis result comprises a mass matrix and a stiffness matrix of nonlinear displacement of a structural node and nonlinear static deformation of the wing;
the step S3 comprises the following steps:
s31: according to the quality matrix and the rigidity matrix of the nonlinear static deformation of the wing, performing quasi-modal analysis in a aeroelastic deformation balance state to obtain the inherent vibration characteristic of the wing;
s32: calculating the unsteady aerodynamic force of the wing curved surface according to the nonlinear displacement of the structural node;
s33: and calculating the flutter speed and frequency by utilizing a Vg method according to the inherent vibration characteristics of the wing and the unsteady aerodynamic force of the curved surface of the wing.
Optionally, in the step S31, the inherent vibration characteristics of the wing include an inherent frequency when the wing can vibrate in only one degree of freedom of heave and an inherent frequency when the wing can vibrate in only one degree of freedom of pitch,
natural frequency of the wing only capable of performing sinking-floating single-degree-of-freedom vibrationThe method comprises the following steps:
natural frequency of the wing when only making pitching single-degree-of-freedom vibrationThe method comprises the following steps:
wherein,and->Stiffness matrix->And of elements in (1)，/>Wing mass in unit span, +.>Mass moment of inertia of unit extended wing to rigid center,>dimensionless turning radius of unit extended wing to rigid center, < >>For the chord length of the wing。
Optionally, the S32 includes:
s321: determining a substantial derivative matrix according to the nonlinear displacement of the structural node;
the substantial derivative matrixThe method comprises the following steps:
wherein,、/>respectively the real part and the imaginary part of the matrix of the substantial derivatives, < >>Is a unit wash down,/->Is the displacement of aerodynamic grid points, +.>For the wash-down matrix induced by initial angle of attack, camber and twist, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
s322: obtaining the unsteady aerodynamic force of the wing curved surface according to the substantial derivative matrix, the aerodynamic force influence coefficient matrix and the integral matrix;
the aerodynamic influence coefficient matrixThe method comprises the following steps:
wherein,is a unit wash down,/->Is->Pressure on individual cells->Is dynamic pressure (or->Indicate->A plurality of object plane surface element control points;
the integration matrixThe method comprises the following steps:
wherein,for the forces and moments acting on aerodynamic grid points +.>Is->Pressure on the individual cell, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
the wing curved surface is not always aerodynamicThe method comprises the following steps:
optionally, the inherent vibration characteristic of the wing includes a natural frequency of the wing, and the S33 includes:
s331: using a vibration equation under a modal coordinate as a flutter equation:
wherein,is the natural frequency of the wing>A mass matrix, a damping matrix and a stiffness matrix, respectively, < >>Is an unsteady aerodynamic influence coefficient matrix, which is the Mach number of the incoming flow>And reduced frequency->Is a function of (2) to reduce the frequencyFrom dimension parameters->、/>、/>Determining the following relation: />；/>For the reference chord length of the wing model, +.>For a dimensional flight speed>For flowing dynamic pressure->Is a modal coordinate magnitude vector;
s332: aerodynamic forces are incorporated into the mass term, resulting in:
s333: to obtain the vibration velocity, the frequency is converted into a dimensional velocity and multiplied by a damping termThe method comprises the following steps of:
wherein,is the incoming flow density; when->When the wing happens to correspond to the critical flutter condition, the characteristic value of the equationThe method comprises the following steps:
damping deviceAnd a dimensional flight speed->Can use the real part of the eigenvalue +.>And imaginary part->Expressed as:
the definition of the reduced frequency can obtain the flutter frequency：
Wherein,is the circumference ratio.
The invention has the following beneficial effects:
the low-reliability static gas bomb analysis method is based on a bin method, the high-reliability static gas bomb analysis method is based on an RANS equation, and a structure solver selects Nastran large-scale finite element analysis software; the pneumatic and structural solving modules are loosely coupled, information transmission is realized through the data exchange module, the dynamic grid technology adopts an RBF/TFI mixed deformation grid method, and the static pneumatic elasticity analysis method and the Nastran structural finite element solver are coupled to form nonlinear flutter analysis content of the wing, so that the applicability and reliability in the wing design process can be enhanced.
Drawings
FIG. 1 is a flow chart of a nonlinear flutter analysis method of the present invention;
FIG. 2 is a schematic diagram of the structure of a mixed mesh deformation method calculation mesh and a background mesh according to the present invention;
fig. 3 is a schematic structural diagram of a deformation unit obtained by using a TFI difference method in the hybrid mesh deformation method of the present invention.
Detailed Description
The principles and features of the present invention are described below with reference to the drawings, the examples are illustrated for the purpose of illustrating the invention and are not to be construed as limiting the scope of the invention.
The invention provides a nonlinear flutter analysis method, which is shown by referring to fig. 1, and comprises the following steps:
s1: acquiring an initial value of a flight parameter;
wherein, the flight parameters include data such as flight altitude, air density, flight speed and the like.
S2: according to the current flight parameters, analyzing the static pneumatic elasticity of the wing structure under large deformation to obtain a static pneumatic elasticity analysis result;
optionally, the static aeroelastic analysis result includes a structural node nonlinear displacement, and the S2 includes:
s21: solving aerodynamic force by using a surface element method/RANS (random access network) consideration transition method according to the current flight parameters;
according to the current flight parameters, solving aerodynamic forces by using a face element method comprises:
at the position ofThe control equation of the unsteady bin method is:
wherein the subscriptExpress object plane->Indicates wake up>Is->All panels at time->Matrix of aerodynamic influence coefficients of object plane dipole intensities on individual object plane elements +.>Is->All panels at time->A matrix of wake dipole intensity aerodynamic influence coefficients on individual object plane elements +.>Is->Moment object plane source strong aerodynamic force influence coefficient matrix; />Is->Dipole intensity of the object plane element at moment +.>Is thatDipole intensity of the moment trail surface element; />Is->The object plane source at the moment is strong，/>Is->A moment object plane surface element vector matrix;、/>and->Respectively->The moment is infinite free incoming flow speed, the movement speed of an object plane due to structural vibration and the disturbance speed of gusts; subscript->Represents free incoming flow at infinity,>representing the vibration of the object plane structure->Indicating gusts.
The intensity distribution of the object plane surface element can be obtained through a control equation, and then the pressure coefficient acting on the j object plane surface element is as follows:
wherein:is->Pressure coefficients at the individual object plane element control points; />And->Acting at +.>Static pressure at the control point of each object plane element and static pressure of free incoming flow at infinity; />Is the atmospheric density; />Is->Local fluid velocities at the individual object plane bin control points; />Is->Velocity profile at the control point of the individual object plane element,/->For time (I)>Is the sign of the partial derivative.
According to the current flight parameters, the RANS is utilized to consider a transition method, and solving aerodynamic force comprises the following steps:
the three-dimensional compressible unsteady N-S equation under the coordinate system is calculated, and can be written as follows:
wherein,is the coordinate of any point in the flow field, +.>Is a flow field conservation variable +.>、/>And->Respectively->Non-stick flux in three directions, +.>，/>And->Respectively->Three-directional adhesive flux, +.>For air density->Speed is +.>Components in three directions>Is the total internal energy per unit volume.
To solve the control equation, time and space dispersion is required. The time dispersion adopts an implicit propulsion method, the space dispersion adopts a finite volume method, and the calculation of the non-viscous flux adopts a Roe format of three-order windward MUSCL interpolation. A spark-Allmaras one-equation turbulence model (SA turbulence model for short) is used to close the N-S equation set. In addition, multiple grids, residual value fairing and parallel computing technology are introduced to improve the computing solution efficiency.
Transition model is implemented by enabling effective intermittent factors +.>Superimposed on the generation term and dissipation term of the turbulent energy equation to simulate and control the transition process. The flow field variables are solved by the transport equation.
The intermittent factorThe transport equation in conservation form is:
wherein,for density (I)>For time (I)>To generate items->For dissipating items->Molecular viscosity>For vortex viscosity>Is constant and->=1.0；/>The values are 1, 2, 3, < >>Is->And->Coordinates representing three directions, +.>Is->And->Is->Speeds in three directions; />Is the sign of the partial derivative; for->I.e. +.>，And the same is true.
S22: interpolating the aerodynamic force into target structural points of the wing structure to obtain the force of each structural point;
the aerodynamic force is interpolated into target structural points of the wing structure by using a radial basis function interpolation method, so that the force of each structural point is obtained.
The radial basis function interpolation method specifically comprises the following steps:
when the radial basis function interpolation (RBF) method is used for carrying out data conversion on discrete points, the method is simple, accurate, rapid and time-saving, and can complete data exchange among various models by utilizing the coordinate relation of the discrete points. The interpolation function of RBF is:
wherein,is->Interpolation function value at coordinates, +_>For the total number of interpolation nodes>For interpolation node +.>Is->Coordinates of interpolation nodes, +_>Is->Weight coefficient corresponding to each interpolation node, < ->As a radial basis function, whereinRespectively indicated abscissa, ordinate and altitude,/->Representing a linear polynomial and being 0 in the mesh deformation problem;
if it isEqual to the displacement at the boundary grid point, then there is:
wherein the subscriptRepresenting the flow field junction>、/>And->Respectively indicate->The boundary nodes of the object plane areA displacement component of the direction, and->、/>、/>，Representing boundary grid points, ">Indicating that the 1 st object plane boundary node is +.>Displacement component of direction, ++>Represent the firstThe boundary point of the object plane is +.>Displacement component of direction, ++>Indicating that the 1 st object plane boundary node is +.>Displacement component of direction, ++>Indicate->The boundary point of the object plane is +.>Displacement component of direction, ++>Indicating that the 1 st object plane boundary node is atDisplacement component of direction, ++>Indicate->The boundary point of the object plane is +.>Displacement component of direction, ++>、/>And->Respectively indicate->A matrix of weight coefficients with undetermined direction, and +.>、/>、，/>An interpolation matrix representing object plane nodes, an，/>Radial basis function representing object plane node, +.>And->Representing two boundary grid points, and the radial basis function is equal to +.>And->The distance expression between them is:
；
if it isAs with the displacement of the spatial grid points, there are:
wherein,、/>and->Respectively indicate->The flow field of the individual space is->Displacement component of direction, subscript->Representing the spatial flow field,/->For the space grid point->Interpolation matrix for spatial flow field and，/>representing the radial basis functions of the spatial flow field,
according toAnd->Obtaining boundary grid points and space grid pointsThe displacement relation between the two is as follows:
。
it can be seen from this that,the matrix is the key to acquiring the boundary and spatial grid point displacements. The radial basis function interpolation method has a simple data structure. However, in practical application, the border grid quantity +.>And space grid quantity->Is quite bulky, making the process particularly inefficient.
Therefore, the invention utilizes the mixed grid deformation method to transmit the structural deformation displacement to the whole calculation space domain, and forms a set of sparse background grids.
The hybrid mesh morphing method combines a high-robustness RBF method with a high-efficiency TFI method. The specific implementation steps are as follows:
(1) According to the focus of the aeroelastic deformation research, selecting proper control points from space and boundary grid points (the control points can be increased in the expanding direction for bending problems, and chord direction control points are required to be increased for torsion problems so as to ensure the grid quality after deformation), forming a set of sparse background grid, wherein the total amount of the grid is as shown in figure 2，/>Taking one control point from four points at intervals in three directions to obtain a +.>Is a background grid of (c).
(2) An RBF interpolation method is adopted forDeforming each vertex of the background mesh;
(3) Based on the high-quality deformed background grid obtained in (2), a deformation unit (as shown in fig. 3, the deformation unit isA local tile in the background grid) displacement of other grid nodes. Thus, the dynamic mesh deformation is ended.
The infinite interpolation method (Transfinite Interpolation, TFI) is a multidimensional interpolation method that can transfer the displacement of the object plane boundary points into the whole computational spatial domain by simple algebraic interpolation. Since this part is the prior art, the present invention will not be described in detail.
S23: obtaining the displacement of each current structural point according to the relation between the stress condition of the force of each structural point and the rigidity matrix;
specifically, the invention calculates the displacement of each current structural point by using Nastran, namely, the relation between the stress condition of each structural point and the rigidity matrix is input into Nastran for nonlinear calculation, and the output result comprises the displacement of each current structural point.
S24: and judging whether the deformation of the current wing structure is converged or not according to the displacement of each current structural point, if so, outputting the displacement of each current structural point as the nonlinear displacement of the structural node, otherwise, updating the pneumatic surface and returning to S21.
According to the invention, the aerodynamic force is recalculated by substituting the displacement of each current structural point into the geometric deformation to obtain new displacement, whether the difference between the new displacement and the displacement of each current structural point is within a preset threshold value is judged, if so, the deformation of the current wing structure is converged, otherwise, the deformation is not converged.
S3: carrying out flutter analysis according to the static pneumatic elasticity analysis result to obtain flutter speed and frequency;
optionally, the static pneumatic elastic analysis result comprises a mass matrix and a stiffness matrix of nonlinear displacement of a structural node and nonlinear static deformation of the wing;
the Nastran nonlinear calculation result file '. DBALL' file and '. MASTER' file comprise nonlinear static deformation calculation quality matrix and rigidity matrix after static deformation convergence. And adding a DAMP sentence in the linear flutter calculation '. BDF' file, and carrying out flutter analysis by restarting, wherein the finally obtained flutter speed is the nonlinear flutter speed.
Based on this, the S3 includes:
s31: according to the quality matrix and the rigidity matrix of the nonlinear static deformation of the wing, performing quasi-modal analysis in a aeroelastic deformation balance state to obtain the inherent vibration characteristic of the wing;
the natural vibration characteristics of the wing include the natural frequency when the wing can only vibrate in a heave single degree of freedom and the natural frequency when the wing can only vibrate in a pitch single degree of freedom,
natural frequency of wing capable of sinking and floating single degree of freedom vibrationThe method comprises the following steps:
natural frequency of wing capable of only pitching single-degree-of-freedom vibrationThe method comprises the following steps:
wherein,and->Stiffness matrix->And of elements in (1)，/>Wing mass in unit span, +.>Mass moment of inertia of unit extended wing to rigid center,>dimensionless turning radius of unit extended wing to rigid center, < >>For the chord length of the wing。
S32: calculating the unsteady aerodynamic force of the wing curved surface according to the nonlinear displacement of the structural node;
the method specifically comprises the following steps:
s321: determining a substantial derivative matrix according to the nonlinear displacement of the structural node;
matrix of substantial derivativesThe method comprises the following steps:
wherein,、/>respectively the real part and the imaginary part of the matrix of the substantial derivatives, < >>Is a unit wash down,/->Is the displacement of aerodynamic grid points, +.>For the wash-down matrix induced by initial angle of attack, camber and twist, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
s322: obtaining the unsteady aerodynamic force of the wing curved surface according to the substantial derivative matrix, the aerodynamic force influence coefficient matrix and the integral matrix;
aerodynamic force influence coefficient matrixThe method comprises the following steps:
wherein,is a unit wash down,/->Is->Pressure on individual cells->Is dynamic pressure (or->Indicate->A plurality of object plane surface element control points;
the integration matrixThe method comprises the following steps:
wherein,for the forces and moments acting on aerodynamic grid points +.>Is->Pressure on the individual cell, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
the wing curved surface is not always aerodynamicThe method comprises the following steps:
s33: and calculating the flutter speed and frequency by utilizing a Vg method according to the inherent vibration characteristics of the wing and the unsteady aerodynamic force of the curved surface of the wing.
Based on the above technical solution, S33 includes:
s331: using a vibration equation under a modal coordinate as a flutter equation:
wherein,is the natural frequency of the wing>A mass matrix, a damping matrix and a stiffness matrix, respectively, < >>Is an unsteady aerodynamic influence coefficient matrix, which is the Mach number of the incoming flow>And reduced frequency->Is a function of (2) to reduce the frequencyFrom dimension parameters->、/>、/>Determining the following relation: />；/>For the reference chord length of the wing model, +.>For a dimensional flight speed>For flowing dynamic pressure->Is a modal coordinate magnitude vector;
s332: aerodynamic forces are incorporated into the mass term, resulting in:
s333: to obtain the vibration velocity, the frequency is converted into a dimensional velocity and multiplied by a damping termThe method comprises the following steps of:
wherein,is the incoming flow density; when->When the wing happens to correspond to the critical flutter condition, the characteristic value of the equationThe method comprises the following steps: />
Damping deviceAnd a dimensional flight speed->Can use the real part of the eigenvalue +.>And imaginary part->Expressed as:
the definition of the reduced frequency can obtain the flutter frequency：
Wherein,is the circumference ratio.
S4: judging whether the flutter speed is matched with the incoming flow speed, if so, entering S5, otherwise, entering S6;
if the obtained flutter velocity does not coincide with the incoming flow velocity when the static aerobullet analysis is carried out, matching iteration is needed, namely, after the flutter analysis is finished, the incoming flow velocity used in the static aerobullet analysis is updated. And then carrying out iterative computation until the flutter speed is matched with the incoming flow speed, and ending the iteration, wherein the flutter speed at the moment is the real nonlinear flutter speed under the attack angle.
Selecting proper according to flight environment、/>、/>Value, will be the formulaThe method is regarded as solving the complex eigenvalue problem to obtain a complex eigenvalue:
this complex eigenvalue may be determined byAnd->Is indicative of the incoming flow velocity +.>The method comprises the following steps:
wherein,to the incoming flow density.
S5: outputting the flutter speed and the frequency as nonlinear flutter analysis results;
s6: and (2) adjusting the flight parameters according to the flutter speed and returning to the step (S2).
The foregoing description of the preferred embodiments of the invention is not intended to limit the invention to the precise form disclosed, and any such modifications, equivalents, and alternatives falling within the spirit and scope of the invention are intended to be included within the scope of the invention.

Claims (7)

1. A nonlinear chatter analysis method, characterized in that the nonlinear chatter analysis method comprises:
s1: acquiring an initial value of a flight parameter;
s2: according to the current flight parameters, analyzing the static pneumatic elasticity of the wing structure under large deformation to obtain a static pneumatic elasticity analysis result;
s3: carrying out flutter analysis according to the static pneumatic elasticity analysis result to obtain flutter speed and frequency;
s4: judging whether the flutter speed is matched with the incoming flow speed, if so, entering S5, otherwise, entering S6;
s5: outputting the flutter speed and the frequency as nonlinear flutter analysis results;
s6: adjusting flight parameters according to the flutter speed and returning to S2;
the static pneumatic elasticity analysis result comprises nonlinear displacement of a structural node, and the S2 comprises:
s21: solving aerodynamic force by using a surface element method/RANS (random access network) consideration transition method according to the current flight parameters;
s22: interpolating the aerodynamic force into target structural points of the wing structure to obtain the force of each structural point;
s23: obtaining the displacement of each current structural point according to the relation between the stress condition of the force of each structural point and the rigidity matrix;
s24: judging whether the deformation of the current wing structure is converged or not according to the displacement of each current structural point, if so, outputting the displacement of each current structural point as the nonlinear displacement of the structural node, otherwise, updating the pneumatic surface and returning to S21;
the static pneumatic elasticity analysis result comprises a structural node nonlinear displacement, a quality matrix of nonlinear static deformation of the wing and a rigidity matrix;
the step S3 comprises the following steps:
s31: according to the quality matrix and the rigidity matrix of the nonlinear static deformation of the wing, performing quasi-modal analysis in a aeroelastic deformation balance state to obtain the inherent vibration characteristic of the wing;
s32: calculating the unsteady aerodynamic force of the wing curved surface according to the nonlinear displacement of the structural node;
s33: calculating the flutter speed and frequency by utilizing a Vg method according to the inherent vibration characteristics of the wing and the unsteady aerodynamic force of the curved surface of the wing;
the natural vibration characteristics of the wing include a natural frequency of the wing, and the S33 includes:
s331: using a vibration equation under a modal coordinate as a flutter equation:
wherein,is the natural frequency of the wing>A mass matrix, a damping matrix and a stiffness matrix, respectively, < >>Is an unsteady aerodynamic influence coefficient matrix, which is the Mach number of the incoming flow>And reduced frequency->Is reduced in frequency +.>From dimension parameters->、/>、/>Determining the following relation: />；/>Is the reference chord length of the wing model,for a dimensional flight speed>For flowing dynamic pressure->Is a modal coordinate magnitude vector;
s332: aerodynamic forces are incorporated into the mass term, resulting in:
s333: to obtain the vibration velocity, the frequency is converted into a dimensional velocity and multiplied by a damping termThe method comprises the following steps of:
wherein,is the incoming flow density; when->The wing corresponds exactly to the critical flutter situation, the characteristic value of this equation +.>The method comprises the following steps:
damping deviceAnd a dimensional flight speed->Can use the real part of the eigenvalue +.>And imaginary part->Expressed as:
the definition of the reduced frequency can obtain the flutter frequency：
Wherein,is the circumference ratio.
2. The nonlinear flutter analysis method according to claim 1, wherein in S21, solving aerodynamic forces using a face element method according to current flight parameters comprises:
according toControl equation of time unsteady surface element method, determining action on the firstjPressure coefficients on the individual object plane elements;
the saidThe control equation of the time unsteady surface element method is as follows:
wherein the subscriptExpress object plane->Indicates wake up>Is->All panels at time->Matrix of aerodynamic influence coefficients of object plane dipole intensities on individual object plane elements +.>Is->All panels at time->A matrix of wake dipole intensity aerodynamic influence coefficients on individual object plane elements +.>Is->Moment object plane source strong aerodynamic force influence coefficient matrix; />Is->Dipole intensity of the object plane element at moment +.>Is->Dipole intensity of the moment trail surface element; />Is->The object plane source at the moment is strong，/>Is->A moment object plane surface element vector matrix;、/>and->Respectively->The moment is infinite free incoming flow speed, the movement speed of an object plane due to structural vibration and the disturbance speed of gusts; subscript->Represents free incoming flow at infinity,>representing the vibration of the object plane structure->Indicating gusts;
the pressure coefficient is as follows:
wherein:is->Pressure coefficients at the individual object plane element control points; />And->Acting at +.>Static pressure at the control point of each object plane element and static pressure of free incoming flow at infinity; />Is the atmospheric density; />Is->Local fluid velocities at the individual object plane bin control points; />Is->Velocity profile at the control point of the individual object plane element,/->Time is; />Is the sign of the partial derivative;
determining the pneumatic load of the surface element according to the pressure coefficient and the current flight parameter;
pneumatic loading of the doughThe method comprises the following steps:
wherein,for reference area->Is->The surface element areas of the individual object surfaces;
and obtaining aerodynamic force and aerodynamic moment according to the pressure coefficient and the aerodynamic load of the surface element.
3. The nonlinear flutter analysis method according to claim 1, wherein in S21, according to the current flight parameter, solving aerodynamic forces by using the RANS consideration transition method comprises:
solving a flow field conservation variable according to a three-dimensional compressible unsteady N-S equation and an intermittent factor transportation equation;
the three-dimensional compressible unsteady N-S equation is:
wherein,is the coordinate of any point in the flow field, +.>Is a flow field conservation variable +.>、/>And->Respectively->Non-stick flux in three directions, +.>，/>And->Respectively->Three-directional adhesive flux, +.>For air density->Speed is +.>Components in three directions>Total internal energy per unit volume;
the intermittent factorThe transport equation in conservation form is:
wherein,for density (I)>For time (I)>To generate items->For dissipating items->Molecular viscosity>For vortex viscosity>Is constant and->=1.0；/>The values are 1, 2, 3, < >>Is->And->Coordinates representing three directions, +.>Is->And->Is->Speeds in three directions;
determining a pressure coefficient and a lift coefficient according to the flow field conservation variable;is the sign of the partial derivative;
the pressure coefficientAnd lift coefficient->The method comprises the following steps of:
wherein,is surface static pressure>Is free incoming flow static pressure, +.>Is the free-flowing dynamic pressure of the incoming flow,is the lower surface pressure coefficient,/->Is the upper surface pressure coefficient,/-, is->Is the normalized chordwise location;
coefficient of surface frictionIs surface shear stress->And come and go freelyRatio of flow pressure:
and obtaining aerodynamic force according to the pressure coefficient, the lift coefficient and the surface friction coefficient.
4. The nonlinear flutter analysis method according to claim 1, wherein in S22, the aerodynamic force is interpolated into target structural points of the wing structure by using a radial basis function interpolation method to obtain forces of the respective structural points.
5. The nonlinear flutter analysis method according to any one of claims 1 to 4, wherein in S24, the aerodynamic force is recalculated by substituting the displacement of each current structural point into the geometric deformation to obtain a new displacement, and whether the difference between the new displacement and the displacement of each current structural point is within a preset threshold value is determined, if yes, the deformation of the current wing structure converges, otherwise, the deformation does not converge.
6. The nonlinear flutter analysis method according to claim 1, wherein in S31, the natural vibration characteristics of the wing include a natural frequency when the wing can vibrate in only a heave single degree of freedom and a natural frequency when the wing can vibrate in only a pitch single degree of freedom,
natural frequency of the wing only capable of performing sinking-floating single-degree-of-freedom vibrationThe method comprises the following steps:
natural frequency of the wing when only making pitching single-degree-of-freedom vibrationThe method comprises the following steps:
wherein,and->Stiffness matrix->Element of->，Wing mass in unit span, +.>Mass moment of inertia of unit extended wing to rigid center,>dimensionless turning radius of unit extended wing to rigid center, < >>Is +.>。
7. The non-linear chatter analysis method of claim 1, wherein said S32 comprises:
s321: determining a substantial derivative matrix according to the nonlinear displacement of the structural node;
the essence is thatDerivative matrixThe method comprises the following steps:
wherein,、/>respectively the real part and the imaginary part of the matrix of the substantial derivatives, < >>Is a unit wash down,/->Is the displacement of aerodynamic grid points, +.>For the wash-down matrix induced by initial angle of attack, camber and twist, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
s322: obtaining the unsteady aerodynamic force of the wing curved surface according to the substantial derivative matrix, the aerodynamic force influence coefficient matrix and the integral matrix;
the aerodynamic influence coefficient matrixThe method comprises the following steps:
wherein,is a unit wash down,/->Is->Pressure on individual cells->Is dynamic pressure (or->Indicate->A plurality of object plane surface element control points;
the integration matrixThe method comprises the following steps:
wherein,for the forces and moments acting on aerodynamic grid points +.>Is->Pressure on the individual cell, subscript +.>Indicate->Aerodynamic grid points>Indicate->A plurality of object plane surface element control points;
the wing curved surface is not always aerodynamicThe method comprises the following steps:
。
CN202311216406.9A 2023-09-20 2023-09-20 Nonlinear flutter analysis method Active CN116956782B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311216406.9A CN116956782B (en) 2023-09-20 2023-09-20 Nonlinear flutter analysis method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311216406.9A CN116956782B (en) 2023-09-20 2023-09-20 Nonlinear flutter analysis method

Publications (2)

Publication Number Publication Date
CN116956782A CN116956782A (en) 2023-10-27
CN116956782B true CN116956782B (en) 2023-12-01

Family

ID=88455046

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311216406.9A Active CN116956782B (en) 2023-09-20 2023-09-20 Nonlinear flutter analysis method

Country Status (1)

Country Link
CN (1) CN116956782B (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103310060A (en) * 2013-06-19 2013-09-18 西北工业大学 Transonic limit cycle flutter analysis method
CN106096088A (en) * 2016-05-31 2016-11-09 中国航空工业集团公司西安飞机设计研究所 A kind of propeller aeroplane WHIRL FLUTTER ANALYSIS method
CN110287505A (en) * 2019-03-20 2019-09-27 北京机电工程研究所 Stability of aircraft analysis method
CN113111430A (en) * 2021-03-06 2021-07-13 北京航空航天大学 Elastic aircraft flight dynamics modeling method based on nonlinear aerodynamic order reduction
CN115422654A (en) * 2022-08-21 2022-12-02 西北工业大学 CFD/CSD technology-based efficient high-precision flutter time domain analysis method for cross/supersonic aircraft

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2995426B1 (en) * 2012-09-11 2014-09-05 Airbus Operations Sas METHOD FOR SIMULATION OF INSTATIONARIAN AERODYNAMIC LOADS ON AN EXTERNAL AIRCRAFT STRUCTURE
US20220300672A1 (en) * 2021-03-19 2022-09-22 Government Of The United States, As Represented By The Secretary Of The Air Force Methods and Systems for Analyzing and Predicting Aeroelastic Flutter on Configurable Aircraft

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103310060A (en) * 2013-06-19 2013-09-18 西北工业大学 Transonic limit cycle flutter analysis method
CN106096088A (en) * 2016-05-31 2016-11-09 中国航空工业集团公司西安飞机设计研究所 A kind of propeller aeroplane WHIRL FLUTTER ANALYSIS method
CN110287505A (en) * 2019-03-20 2019-09-27 北京机电工程研究所 Stability of aircraft analysis method
CN113111430A (en) * 2021-03-06 2021-07-13 北京航空航天大学 Elastic aircraft flight dynamics modeling method based on nonlinear aerodynamic order reduction
CN115422654A (en) * 2022-08-21 2022-12-02 西北工业大学 CFD/CSD technology-based efficient high-precision flutter time domain analysis method for cross/supersonic aircraft

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"机翼跨声速非线性颤振及高效分析方法研究";刘南;《中国博士学位论文全文数据库 工程科技II辑》;第C031-4页 *
大展弦比机翼的几何非线性气动弹性分析;万仲;张军红;韩景龙;;江苏航空(S1);全文 *

Also Published As

Publication number Publication date
CN116956782A (en) 2023-10-27

Similar Documents

Publication Publication Date Title
CN113111430B (en) Elastic aircraft flight dynamics modeling method based on nonlinear aerodynamic order reduction
CN109933876B (en) Unsteady aerodynamic order reduction method based on generalized aerodynamic force
CN104182560B (en) aircraft flutter prediction analysis method and device
Sayed et al. Aeroelastic analysis of 10 MW wind turbine using CFD–CSD explicit FSI-coupling approach
Raveh CFD-based models of aerodynamic gust response
Murua 'Flexible aircraft dynamics with a geometrically-nonlinear description of the unsteady aerodynamics'
CN112580241B (en) Nonlinear aeroelastic dynamic response analysis method based on structure reduced order model
Hang et al. Analytical sensitivity analysis of flexible aircraft with the unsteady vortex-lattice aerodynamic theory
Wang et al. High-order simulation of aeronautical separated flows with a Reynold stress model
CN113094837A (en) Wind resistance design method of horizontal axis wind turbine blade under strong wind action
Artola et al. Modal-based nonlinear estimation and control for highly flexible aeroelastic systems
Imiela et al. Towards multidisciplinary wind turbine design using high-fidelity methods
CN116956782B (en) Nonlinear flutter analysis method
Luo et al. A 3D computational study of the flow-structure interaction in flapping flight
Xiong et al. Aeroelastic Modeling and CFD Simulation of Wind-Tunnel Scale Aspect Ratio 13.5 Common Research Model
Huang et al. Numerical method of static aeroelastic correction and jig-shape design for large airliners
Nguyen et al. Multi-point jig twist optimization of mach 0.745 transonic truss-braced wing aircraft and high-fidelity cfd validation
CN116956781B (en) RANS-based transition-considered static pneumatic elasticity analysis method
Wang Transonic static aeroelastic and longitudinal aerodynamic characteristics of a low-aspect-ratio swept wing
Luo et al. Strongly coupled fluid–structure interaction analysis of aquatic flapping wings based on flexible multibody dynamics and the modified unsteady vortex lattice method
Xiao et al. Wing Flutter Simulations Using an Aeroelastic Solver Based on the Predictor—Corrector Scheme
Zuo et al. Efficient aeroelastic design optimization based on the discrete adjoint method
Gov et al. Geometrically nonlinear model for gust response of very flexible wings using segmental modes
Liu et al. Experimental and numerical studies on a passively deformed coupled-pitching hydrofoil under the semi-activated mode
Chaoyuan et al. Aerodynamic characteristics of a pitching airfoil with leading-edge morphing

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant