An ultrasonic speed flutter analysis method based on a viscous local flow piston theory

Abstract

The present application belongs to the technical field of flutter analysis of supersonic vehicles, and proposes a supersonic flutter analysis method based on the viscous local flow piston theory. The method obtains the interference flow field of large attack angle and complex shape by numerically solving the N-S equation through CFD, and then obtains the effective aerodynamic shape of the vehicle according to the vorticity criterion; the surface of the vehicle is discretized to obtain a high-fidelity three-dimensional aerodynamic grid suitable for the piston theory; the flow field parameters and structural modal shapes at the effective shape of CFD are interpolated onto the piston grid by using the thin plate spline interpolation method, the generalized aerodynamic influence coefficient matrix is calculated, and the state space form of the structural dynamics equation is assembled; the flutter analysis is carried out by analyzing the eigenvalues of the state space equation. The analysis method can analyze any complex aerodynamic shape and consider the viscous effect of real gas, and has high calculation accuracy and efficiency.

Description

Technical Field

[0001] This invention belongs to the field of flutter analysis technology for supersonic high-speed vehicles, and specifically relates to a supersonic flutter analysis method based on the viscous local flow piston theory. Background Technology

[0002] Supersonic vehicles, due to their significant advantage of high speed, can effectively achieve rapid long-range deployment, efficient penetration, and high survivability, making them of extremely high strategic value in the field of military aviation. The next generation of supersonic and hypersonic vehicles under development will extensively utilize lightweight materials and employ large, thin-walled structural designs, resulting in lower natural vibration frequencies and more pronounced fluid-structure interaction problems. Therefore, flutter analysis is of great importance in the design of supersonic and hypersonic vehicles.

[0003] One of the key technologies in supersonic flutter analysis is determining unsteady aerodynamic forces. Previous algorithms can be divided into two categories: high-precision numerical simulation methods based on CFD / CSD coupling and engineering methods based on approximate theories, including piston theory and unified lifting surface theory. Numerical simulation methods suffer from high computational cost and low efficiency, making them unsuitable for preliminary aircraft design; while engineering methods struggle to analyze complex shapes such as wing-body assemblies and large angles of attack. Chinese invention patent CN104133933B discloses a hypersonic aeroelastic analysis method based on piston theory, but this method requires simplifying the wing-body assembly into a two-dimensional aerodynamic mesh, failing to consider the thickness effects of the wing and fuselage and the three-dimensional effects of low aspect ratio surfaces.

[0004] Based on the problems of the above two methods, the local flow piston theory, which combines high-precision CFD and piston theory, was developed in the 1990s, achieving a balance between computational efficiency and accuracy. In June 1995, Yang Bingyuan and Song Weili, in their paper "Calculation of Supersonic Flutter of High Angle-of-Attack Airfoils Using Local Flow Piston Theory" published in Volume 14, Issue 2 of "Vibration and Shock," first introduced CFD flow field parameters into the local flow piston theory and conducted flutter analysis. In September 2005, Zhang Weiwei, Ye Zhengyin, et al., in their paper "Research on Aeroelastic Calculation Method Based on Local Flow Piston Theory" published in Volume 37, Issue 5 of "Acta Mechanica Sinica," re-derived the local flow piston theory using the momentum theorem, laying the theoretical foundation for extending the local flow piston theory from two-dimensional airfoils to three-dimensional complex shapes. Chinese invention patent CN106508028B describes a supersonic and hypersonic flutter analysis method based on local flow piston theory and applicable to complex shapes. However, it does not consider the influence of viscous boundary layer on effective aerodynamic shape and only derives the generalized aerodynamic calculation method for symmetrical airfoil three-dimensional airfoil and rotating fuselage.

[0005] To address the shortcomings of existing technologies, this invention achieves aerodynamic calculations for arbitrarily complex aerodynamic shapes using a high-fidelity three-dimensional aerodynamic mesh, introduces effective shape correction to account for the viscous effects of real gases, and employs Kriging interpolation of the near-wall flow field to improve the efficiency of effective shape correction. Based on this, a supersonic flutter analysis method based on viscous local flow piston theory is proposed. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention proposes a supersonic flutter analysis method based on viscous local flow piston theory. This method obtains the perturbation flow field under large angles of attack and complex shapes by numerically solving the Navier-Stokes equations using CFD, and then obtains the effective aerodynamic shape of the aircraft based on the vorticity criterion. The aircraft surface is discretized to obtain a high-fidelity three-dimensional aerodynamic mesh suitable for piston theory. Thin-plate spline interpolation is used to interpolate the flow field parameters and structural mode shapes at the effective CFD shape onto the piston mesh, calculates the generalized aerodynamic influence coefficient matrix, and assembles the state-space structural dynamic equations. Flutter analysis is performed by analyzing the eigenvalues ​​of the state-space equations. This analysis method can analyze arbitrarily complex aerodynamic shapes and consider the viscous effects of real gases, balancing computational accuracy and efficiency.

[0007] The specific technical solution of the present invention is as follows:

[0008] A method for analyzing supersonic flutter based on the viscous in-flow piston theory includes the following steps:

[0009] Step S1: For the target aircraft, establish a structural finite element model, analyze and extract the translational modes, generalized mass M and generalized stiffness K of each node on the surface of the aircraft.

[0010] Step S2: For the target aircraft, establish a CFD model of the external flow field, solve the NS equations, and derive the near-wall flow field parameters;

[0011] Step S3: Interpolate the flow field near the wall to obtain the flow field parameters of each grid point on the wall, and calculate the effective shape of the aircraft under viscous effect according to the vorticity criterion.

[0012] Step S4: Establish an interpolation function based on the thin plate spline method to achieve parameter transfer between different grids;

[0013] Step S5: Calculate the generalized aerodynamic influence coefficient matrix based on the local flow piston theory;

[0014] Step S6: Assemble the state-space equations and analyze the flutter characteristics of the aircraft using the eigenvalues ​​of the state transition matrix.

[0015] Furthermore, the translational mode expression in step S1 is as follows:

[0016]

[0017] Where i represents the node number on the aircraft surface, and j represents the modal order. U represents the translational components of the i and j-th order mode shapes at surface nodes. xij u yij and u zij These represent the translational components of the j-th modal array of the i-th aircraft surface node in the global coordinate system of the structural finite element model in the x-axis, y-axis, and z-axis directions, respectively.

[0018] Transform to piston mesh coordinate system:

[0019]

[0020] Among them, L ps This represents the transformation matrix from the coordinate system of the structural finite element model to the coordinate system of the piston mesh.

[0021] Furthermore, the near-wall flow field parameters in step S2 are expressed as follows:

[0022] P o =[p o ρ o v xo v yo v zo ξ o a o ]

[0023] Among them, [v xo ,v yo ,v zo [p] represents the airflow velocity at node o in the flow field. o ρ o ξ o a o These represent the pressure, gas density, flow field vorticity, and local sound speed at node o in the flow field.

[0024] Furthermore, step S3 specifically includes:

[0025] Step S3-1: For the wall mesh point m in the CFD model, let its coordinates be x. m The wall normal vector is n m Extract all flow field nodes located above the wall grid point m in the near-wall flow field and use them as the data point set for flow field interpolation.

[0026] The data point set is extracted using the following method:

[0027]

[0028] Among them, U m Let d represent the data point set, and let x represent the flow field node in the data point set. d I3 represents the coordinates of the flow field node d, I3 represents the 3×3 identity matrix, and ε represents the preset threshold.

[0029] Step S3-2: Obtain the flow field parameters at any distance s along the normal direction of the wall grid point m using the ordinary Kriking method, expressed as:

[0030]

[0031] Where, x s Let λ represent the coordinates at point s. d P represents the weight. d This represents the flow field parameter vector at node d near the wall in the CFD, and its components are defined similarly to the near-wall flow field parameters P in step S2. o same;

[0032] Weight vector λ = [λ d We obtain the following formula:

[0033]

[0034] Where e =

[111] is a column vector consisting entirely of 1s, and the matrix... is the covariance matrix between known nodes, and vector r is the covariance between known nodes and nodes to be determined;

[0035] Step S3-3: Obtain the effective shape and the distance s between the wall surface according to the vorticity criterion. c , is represented as:

[0036]

[0037]

[0038]

[0039] Ma ∞ c ∞ ρ ∞ T ∞ μ ∞ These are the Mach number, sound velocity, density, temperature, and dynamic viscosity of the incoming flow from a distance, T. w Here, T is the wall temperature, T' is the Anderson reference temperature, and x is the wall temperature. m L represents the distance from grid point m on the wall to the leading edge of the wing. mξ represents the reference length, i.e., the wing chord length passing through grid point m on the wall. a1 and a2 are effective shape fitting coefficients; for a circular wing, a1 = 12.05, a2 = -0.8; for a rhomboid wing, a1 = 9.05, a2 = 0.19. c Let ξ(x) be the vorticity criterion at point m on the wall grid. s Let be the vorticity at any distance s along the normal direction from point m on the wall grid.

[0040] Let s gradually increase from 0 until it satisfies ξ(x) s )=ξ c At this time, s c That is, the distance s between the effective shape and the wall. c ;

[0041] Step S3-4: Assign the flow field parameters at the effective shape to the corresponding wall grid points;

[0042] Step S3-5: Perform steps 3-1 to 3-4 on the wall mesh points of all CFD models;

[0043] Step S3-6: Transform the flow field parameters of the wall grid points in the updated CFD model to the piston grid coordinate system;

[0044] Represented as:

[0045]

[0046] The converted flow field velocity is

[0047]

[0048] The transformed coordinates are

[0049]

[0050] Among them, [v xm ,v ym ,v zm [p] represents the airflow velocity at grid point m on the wall. m ρ m ξ m a m The values ​​for pressure, gas density, flow field vorticity, and local sound speed at grid point m on the wall are L, respectively. pa This represents the transformation matrix from the CFD model coordinate system to the piston mesh coordinate system. It is the coordinate of the origin of the CFD model in the piston mesh coordinate system.

[0051] Furthermore, step S4 specifically includes:

[0052] The flow field parameters at the wall mesh points of the updated CFD model are interpolated using the thin-plate spline method. Interpolate to the piston mesh and obtain the displacement interpolation matrix G from the structural finite element model to the piston mesh.

[0053] Furthermore, the expression for the generalized aerodynamic influence coefficient matrix in step S5 is as follows:

[0054]

[0055]

[0056] Where F represents the aerodynamic vector, and Let represent the generalized aerodynamic damping and generalized aerodynamic stiffness matrices, respectively; q represent the structural modal coordinate vector; Φ is the modal translational vibration mode matrix at the wall mesh points of the structural finite element model; S = diag(s) is the area weighted matrix of the piston mesh elements; N = diag(n) is a diagonal matrix composed of the normals n of each piston mesh element; A0 = diag(ρa) is a diagonal matrix composed of the product of the local air density ρ and the sound velocity a of each piston mesh element; A k =A0V k k = 1, 2, 3, where For each piston grid cell, the local flow velocity is at x k Components of direction The resulting diagonal matrix has x1, x2, and x3 directions corresponding to the x-axis, y-axis, and z-axis directions, respectively.

[0057] Furthermore, step S6 specifically includes:

[0058] Step S6-1: The dynamic equations in state space, with the structural modal coordinates q as the generalized displacements, are as follows:

[0059]

[0060] Where M and K are the generalized mass matrix and the generalized stiffness matrix, respectively;

[0061] Step S6-2: According to linear system theory, the necessary and sufficient condition for a linear time-invariant system not to experience chattering is: the system's state matrix... All eigenvalues ​​λ r (A) are all located in the left half of the complex plane, that is

[0062] Re[λ r [A]<0, r=1,2,...,n

[0063] In the formula, n is the number of eigenvalues ​​of matrix A;

[0064] Step S6-3: Calculate the speed state corresponding to each working condition under multiple operating conditions. and And its corresponding state matrix and eigenvalues, thereby obtaining the damping coefficient and vibration frequency corresponding to each velocity state;

[0065] Step S6-4: Plot the real and imaginary parts of the eigenvalues ​​corresponding to the velocities from low to high as root locus diagrams; plot the damping coefficients corresponding to the eigenvalues ​​as VG diagrams; plot the vibration frequencies corresponding to the eigenvalues ​​as VF diagrams; and find the first unstable velocity point as the flutter velocity boundary.

[0066] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0067] 1. Compared with existing CFD / CSD coupled algorithms, the supersonic flutter analysis method based on viscous local flow piston theory proposed in this invention can effectively improve the problems of large computational load and low computational efficiency.

[0068] 2. Compared with existing supersonic flutter analysis methods based on piston theory, the supersonic flutter analysis method based on viscous local flow piston theory proposed in this invention considers the viscosity effect of real gas, is suitable for high-fidelity three-dimensional aerodynamic meshes, and introduces effective automatic shape correction based on flow field interpolation, which improves accuracy and is more convenient for practical engineering applications.

[0069] 3. The supersonic flutter analysis method proposed in this invention, based on the viscous local flow piston theory, is not limited by complex aerodynamic shapes or high angle-of-attack flight conditions, and can analyze thermo-aeroelastic stability problems by introducing thermal modes. It is of great significance in the design of supersonic and hypersonic aircraft. Attached Figure Description

[0070] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the embodiments will be briefly introduced below. The features and advantages of the present invention can be more clearly understood by referring to the accompanying drawings. The accompanying drawings are schematic and should not be construed as limiting the present invention in any way. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0071] Figure 1 A schematic diagram of the piston mesh established for this invention;

[0072] Figure 2 A schematic diagram of interpolating the first-order modal translational matrix of the structure onto the piston mesh;

[0073] Figure 3 This is a schematic diagram illustrating effective shape correction based on the flow field parameters of the CFD model.

[0074] Figure 4 A VG diagram reflecting the variation of the damping coefficient of the eigenvalues ​​of the state matrix with the incoming flow velocity;

[0075] Figure 5 A VF diagram reflecting the variation of vibration frequency of the eigenvalues ​​of the state matrix with the incoming flow velocity. Detailed Implementation

[0076] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments of the present invention and the features thereof can be combined with each other.

[0077] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0078] Example 1

[0079] Taking a certain all-moving control surface as an example, the supersonic flutter analysis method based on the viscous local flow piston theory proposed in this invention is used for analysis.

[0080] Step 1: Establish a structural finite element model of a certain all-moving control surface and perform modal analysis.

[0081] The natural frequencies of the first 20 modes of the model are as follows:

[0082] Modal Frequency (Hz) Modal Frequency (Hz) 1 49.17 11 2674.12 2 53.90 12 2763.49 3 78.13 13 3519.24 4 472.65 14 3635.11 5 493.93 15 3905.66 6 713.78 16 4168.34 7 1223.24 17 4469.27 8 1580.67 18 4692.98 9 1750.23 19 5269.25 10 2085.50 20 5866.80

[0083] To facilitate subsequent calculations, the modal array normalizes the generalized mass.

[0084] Step 2: Establish a CFD model of a certain all-moving control surface and solve the NS equations, and make effective shape corrections based on the flow field parameters.

[0085] The flow pressure q ∞ The CFD calculation conditions corresponding to 110 kPa are shown in the table below:

[0086]

[0087] Maintaining the incoming Mach number M in flutter analysis ∞ Incoming flow density ρ ∞ and wall temperature T w The velocity of the incoming flow remains unchanged, but the velocity of the incoming flow, a, is increased. ∞ The incoming flow pressure is controlled, while the other parameters are determined according to the gas law.

[0088] After the calculation is completed, the flow field parameters are substituted into the following formula to obtain the boundary layer thickness s at each node. c

[0089]

[0090]

[0091] In the formula, the gas constant is taken as R = 287.0 J / (kg*K), and the reference length is taken as L. m =0.284m.

[0092] Step 3: Establish the piston mesh, transform the flow field parameters of the wall mesh of the CFD model after the structural finite element modal analysis results and effective shape correction to the coordinate system of the piston mesh, and interpolate them to the centroid of the piston mesh.

[0093] Although both are aerodynamic meshes, the difference between CFD model meshes and piston meshes lies in the fact that the flow field parameters in CFD models are defined at the nodes, while the local velocity, pressure, and normal vector required by piston theory are defined on the mesh. Therefore, the interpolation of modal displacements is from the wall nodes of the structural finite element model mesh to the piston mesh nodes; while the interpolation of flow field parameters is from the effective shape nodes of the CFD model to the centroid of the piston mesh.

[0094] The established piston grid is as follows Figure 1 As shown, it is divided into upper and lower parts, and each grid is only subject to aerodynamic forces on the outer normal side.

[0095]

[0096]

[0097] In the formula, n is the number of grid points with known parameters; w p Let be the parameter components to be interpolated; ε is a given constant, called a parameter. For a general flat function, ε = 10. -2 ~1, for functions with singularity, ε can be taken as 10. -5 ~10 -6 ; Let r be the d-th dimension coordinate of the i-th node; i 2 =||xx i || 2 Let be the distance between the point to be determined and the known point. The coefficients to be determined are obtained by solving the following matrix:

[0098]

[0099] In the formula in h is the distance between two known nodes.j The weighting coefficients corresponding to the j-th node are pre-defined by the calculator. When all h... j When h = 0, the fitted function passes through all node function values ​​exactly; when all h = 0, the fitted function passes through all node function values ​​exactly. j As the value approaches infinity, the approximation function tends to fit the least squares method.

[0100] The interpolation results for the first-order modal array are as follows: Figure 2 As shown, since the nodes of the structural finite element model mesh and the piston mesh are not in one-to-one correspondence, thin plate interpolation theory is needed to interpolate the array of the structural finite element model onto the piston mesh nodes. Since the applicable condition of thin plate interpolation theory is that the known point set and the point set to be interpolated are located on a smooth surface that is approximately overlapping, this invention adopts a partitioned interpolation method, that is, the upper and lower walls of the structural finite element model and the upper and lower walls of the piston mesh are interpolated accordingly.

[0101] Effective shape corrections such as Figure 3 As shown, the effective shape of the CFD model is first calculated based on the flow field parameters of the CFD model. Then, the distance between the wall mesh and the effective shape of the CFD model (i.e., the boundary layer thickness) is interpolated to the piston mesh in sections, and then the corrected three-dimensional piston mesh is constructed.

[0102] Step 4: Calculate the generalized aerodynamic influence coefficient matrix and state matrix based on piston theory.

[0103] The formula for calculating the generalized aerodynamic influence coefficient is as follows:

[0104]

[0105] in

[0106]

[0107] In the formula, n is the number of piston meshes, ρ l1 , ρ l2 …ρ ln Let a be the local air density of each piston grid cell. l1 a l2 …a ln V represents the sound velocity of each piston grid unit. lxn V represents the component of the local flow velocity in the x-direction for each piston grid cell. lyn V represents the component of the local flow velocity in the y-direction for each piston grid cell. lzn This represents the component of the local flow velocity in the z-direction for each piston grid cell.

[0108] Because the modal matrix normalizes the generalized mass, the form of the state matrix is ​​simplified:

[0109]

[0110] In the formula I is a diagonal matrix composed of the angular frequencies of the first 20 modal vibrations of the structure (in rad / s). 20 It is a 20×20 identity matrix.

[0111] Step 5: Perform eigenvalue analysis on the state matrix under different incoming flow pressures and plot the VG and VF diagrams:

[0112] Let matrix A(q) ∞ One eigenvalue of ) is λ j Then the VG graph is Follow q ∞ The curves change accordingly, with the curves corresponding to the first 10 modes as follows: Figure 4 As shown in the figure, the damping ratio of the second mode increases with increasing dynamic pressure, the damping ratio of the third mode decreases with increasing dynamic pressure, and the damping ratio of the other modes remains basically unchanged. When the incoming dynamic pressure q... ∞ When the pressure is 110 kPa, the second mode has already crossed from below the horizontal axis to above it, meaning that the object under analysis has experienced flutter.

[0113] VF diagram is Follow q ∞ A curve that changes with the changes, such as Figure 5 As shown. Since the variation amplitude of higher-order modes in the VG diagram is very small, the focus in the VF diagram is on the first three modes. According to flutter theory, when two modes are coupled, their vibration frequencies will approach each other as the incoming flow pressure increases. Based on this, it can be determined that the analyzed object is under the influence of incoming flow pressure q. ∞ The cause of flutter at 110 kPa is the coupling of the second and third modes.

[0114] Since the state-space equations are increasing in rank, each mode corresponds to two curves on the VG and VF diagrams. When q ∞ When the value is not too large, these two curves correspond to a pair of conjugate complex roots, and therefore completely overlap on the VG and VF plots. Therefore... Figure 4 , Figure 5 The number of curves and the number of modes are equal.

[0115] In summary, by observing whether the root locus plot contains eigenvalues ​​with real values ​​greater than 0, it can be determined whether the analyzed object is fluttering within the calculation range; by observing the modal branches crossing the horizontal axis in the VG plot, the critical flutter velocity and the mode that is crossing the axis can be determined; by observing the modal branches with similar frequencies in the VF plot, the order of the coupled modes causing aircraft flutter and the flutter frequency can be determined.

[0116] In this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral part; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication of two components or the interaction between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0117] In this invention, unless otherwise explicitly specified and limited, "above" or "below" the second feature can include direct contact between the first and second features, or contact between the first and second features through another feature between them. Furthermore, "above," "over," and "on top" of the second feature includes the first feature directly above or diagonally above the second feature, or simply indicates that the first feature is at a higher horizontal level than the second feature. "Below," "below," and "under" the second feature includes the first feature directly below or diagonally below the second feature, or simply indicates that the first feature is at a lower horizontal level than the second feature.

[0118] In this invention, the terms "first," "second," "third," and "fourth" are used for descriptive purposes only and should not be construed as indicating or implying relative importance. The term "multiple" refers to two or more unless otherwise expressly defined.

[0119] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for analyzing supersonic flutter based on the viscous local flow piston theory, characterized in that, Includes the following steps: Step S1: For the target aircraft, establish a structural finite element model, analyze and extract the translational modes, generalized mass M and generalized stiffness K of each node on the surface of the aircraft. Step S2: For the target aircraft, establish a CFD model of the external flow field, solve the NS equations, and derive the near-wall flow field parameters; Step S3: Interpolate the flow field near the wall to obtain the flow field parameters of each grid point on the wall, and calculate the effective shape of the aircraft under viscous effect according to the vorticity criterion. Step S4: Establish an interpolation function based on the thin plate spline method to achieve parameter transfer between different grids; Step S5: Calculate the generalized aerodynamic influence coefficient matrix based on the local flow piston theory; Step S6: Assemble the state-space equations and analyze the flutter characteristics of the aircraft using the eigenvalues ​​of the state transition matrix; Step S3 specifically includes: Step S3-1: For the wall mesh points of the CFD model m Let its coordinates be The wall normal vector is Extract grid points located near the wall in the near-wall flow field. m All the flow field nodes above are used as the data point set for flow field interpolation; The data point set is extracted using the following method: in, Represents a set of data points. d Represents the flow field nodes in the data point set. Represents flow field nodes d coordinates express identity matrix Indicates a preset threshold; Step S3-2: Obtain the wall mesh points using the ordinary kriging method. m any distance along the normal direction s The flow field parameters at that location are expressed as: in, express s The coordinates of the location Indicates weight, Represents flow field nodes d The flow field parameter vector at that location; Weight vector We obtain the following formula: in, It is a column vector consisting entirely of 1s, and a matrix. It is the known covariance matrix between nodes, a vector. Let the covariance be the difference between the known nodes and the nodes to be determined. Step S3-3: Obtain the effective shape and the distance between the wall based on the vorticity criterion. , is represented as: in These are the Mach number, sound velocity, density, temperature, and dynamic viscosity of the incoming flow from a distance. The wall temperature, For Anderson reference temperature, Represents wall grid points m Distance to the leading edge of the wing Indicates the reference length, i.e., the length passing through the wall grid points. m Wing chord length, , All are effective shape fitting coefficients. For wall grid points m The vorticity criterion at the location, For wall grid points m any distance along the normal direction s vorticity at that location; make s Start from 0 and gradually increase until the condition is met. At this time This is the distance between the effective shape and the wall surface. ; Step S3-4: Assign the flow field parameters at the effective shape to the corresponding wall grid points; Step S3-5: Perform steps S3-1 to S3-4 on the wall mesh points of all CFD models.

2. The supersonic flutter analysis method based on the viscous local piston theory according to claim 1, characterized in that, The translational mode expression in step S1 is as follows: in, i Indicates the number of the node on the aircraft surface. j Indicates the modal order. Represents surface nodes i , j Translational components of the first mode shape. , and The structural finite element model represents the first... i The first aircraft surface node j The first modal array in the global coordinate system of the structural finite element x axis, y shaft and z Translational component in the axial direction; Transform to piston mesh coordinate system: in, This represents the transformation matrix from the coordinate system of the structural finite element model to the coordinate system of the piston mesh.

3. The supersonic flutter analysis method based on the viscous local piston theory according to claim 1, characterized in that, The expression for the near-wall flow field parameters in step S2 is: in, flow field nodes o airflow velocity at that location These are the flow field nodes. o The pressure, gas density, flow field vorticity, and local sound speed at the location.

4. The supersonic flutter analysis method based on the viscous local piston theory according to claim 1, characterized in that, Step S3 further includes: Step S3-6: Transform the flow field parameters of the wall grid points in the updated CFD model to the piston grid coordinate system; Represented as: The converted flow field velocity is The transformed coordinates are in, Represents wall grid points m coordinates For wall grid points m airflow velocity at that location These are wall grid points. m The pressure, gas density, flow field vorticity, and local sound speed at that location. This represents the transformation matrix from the CFD model coordinate system to the piston mesh coordinate system. It is the coordinate of the origin of the CFD model in the piston mesh coordinate system.

5. The supersonic flutter analysis method based on the viscous local piston theory according to claim 4, characterized in that, Step S4 specifically includes: The flow field parameters at the wall mesh points of the updated CFD model are interpolated using the thin-plate spline method. Interpolate to the piston mesh and obtain the displacement interpolation matrix G from the structural finite element model to the piston mesh.

6. The supersonic flutter analysis method based on the viscous local piston theory according to claim 5, characterized in that, The expression for the generalized aerodynamic influence coefficient matrix in step S5 is as follows: Where F represents the aerodynamic vector, and Let represent the generalized aerodynamic damping matrix and the generalized aerodynamic stiffness matrix, respectively. Represents the structural modal coordinate vector. It is the modal translational vibration mode matrix at the mesh points on the wall of the structural finite element model. For the area weighting matrix of the piston mesh elements, It is composed of the normals of each piston mesh element A diagonal array; It is determined by the local air density of each piston grid unit. and speed of sound A diagonal matrix composed of products; ,in For the local flow velocity of each piston grid cell Components of direction The diagonal array formed The directions correspond to direction.

7. The supersonic flutter analysis method based on the viscous local piston theory according to claim 6, characterized in that, Step S6 specifically includes: Step S6-1: Using structural modal coordinates The dynamic equation for the generalized displacement in state space takes the following form: in, and These are the generalized mass matrix and the generalized stiffness matrix, respectively; Step S6-2: According to linear system theory, the necessary and sufficient condition for a linear time-invariant system not to experience chattering is: the system's state matrix... All eigenvalues All are located in the left half of the complex plane, that is In the formula n For matrix The number of eigenvalues; Step S6-3: Calculate the speed state corresponding to each working condition under multiple operating conditions. and , and its corresponding state matrix and eigenvalues, thereby obtaining the damping coefficient and vibration frequency corresponding to each velocity state; Step S6-4: Plot the real and imaginary parts of the eigenvalues ​​corresponding to the velocities from low to high as root locus diagrams, plot the damping coefficients corresponding to the eigenvalues ​​as VG diagrams, plot the vibration frequencies corresponding to the eigenvalues ​​as VF diagrams, and find the first unstable velocity point as the flutter velocity boundary.

Images