Static aeroelastic analysis method based on improved coupled cfd technology and local flow piston theory

Abstract

The application provides a static aeroelastic analysis method based on improved coupling CFD technology and local flow piston theory. First, modal analysis is performed on the aircraft to obtain a vibration mode; then, first steady CFD calculation is performed, local flow parameters are extracted, and an aerodynamic stiffness matrix in a modal coordinate system is established by using a piston theory; then, a static aeroelastic balance equation of the aircraft is established, and structural static deformation and balancing quantities of the aircraft are obtained by solving; finally, CFD grids are updated according to the deformation, second steady CFD calculation is performed, and the aerodynamic characteristics after the deformation are obtained. The application only needs two steady CFD calculations, the piston theory is only used for establishing the aerodynamic stiffness matrix to solve the deformation, and the final aerodynamic characteristics are directly obtained by the second CFD, so that the analysis efficiency is improved while the calculation precision is ensured; and the application is suitable for static aeroelastic analysis of complex aircraft such as a fuselage-rudder combination.

Description

Technical Field

[0001] This invention belongs to the field of aeroelasticity technology, specifically involving a static aeroelasticity analysis method based on improved coupled CFD technology and local flow piston theory. Background Technology

[0002] Due to its simplicity and efficiency, the piston theory is widely used in supersonic aeroelasticity problems. The classical piston theory is mainly applicable to thin airfoils. It assumes that the disturbance generated by the airfoil propagates approximately along the normal direction of the object's surface, and that the disturbance at a certain point on the airfoil has a very weak effect on other points. That is, the pressure at a certain point on the airfoil is only related to the downwash velocity at that point.

[0003] The local flow piston theory was first proposed by Chen Jingsong and Yang Bingyuan. This method relaxes the requirements for the object's shape and the angle of attack of the incoming flow. However, it uses the shock wave expansion wave method to calculate local flow parameters, which cannot accurately calculate the aerodynamic forces of complex shapes. Zhang Weiwei et al. proposed a local flow piston theory combined with steady CFD technology (hereinafter referred to as the LPT method). This method obtains local flow parameters by performing a single CFD calculation, and then uses the piston theory to calculate the aerodynamic stiffness matrix and aerodynamic damping matrix, which are coupled with the structural vibration equation for aeroelastic analysis. The existing LPT method directly calculates the aerodynamic characteristics after deformation using the local flow piston theory in static aeroelastic analysis. Under complex shape conditions, the calculation deviation of the pressure core position is large, which is difficult to meet the engineering accuracy requirements.

[0004] Currently, local flow piston theory combined with CFD technology is mainly used for the aeroelasticity analysis of supersonic panels or wing-like components, while the static aeroelasticity problem of complex-shaped aircraft is mainly analyzed by high-precision CFD methods. High-precision CFD methods obtain aerodynamic loads on the object surface by solving the flow field using CFD, then interpolate these loads to the structural mesh and solve for structural deformation. The deformation is then fed back to the CFD mesh for further CFD solving, iterating until the loads and displacements converge. Although high-precision CFD methods offer high accuracy, the iterative process requires multiple CFD solutions, resulting in large computational scale and long processing time. Summary of the Invention

[0005] The purpose of this invention is to solve the problems of large calculation scale, long time consumption and difficulty in meeting the engineering accuracy requirements under complex shape conditions in the existing static aeroelastic analysis technology for supersonic aircraft. Instead, it provides a static aeroelastic analysis method based on improved coupled CFD technology and local flow piston theory. This method can accurately and quickly obtain the static aeroelastic deformation, trim amount and aerodynamic characteristics of the aircraft after deformation through two steady CFD calculations.

[0006] To achieve the above objectives, the technical solution provided by this invention is:

[0007] The static aeroelastic analysis method based on improved coupled CFD technology and local flow piston theory includes the following steps:

[0008] Step 1: Establish the structural finite element model and aerodynamic surface element model of the target aircraft, perform modal analysis on the structural finite element model, and obtain the mode shapes of each order;

[0009] Step 2: Establish a CFD mesh for the undeformed state of the aircraft, perform the first steady flow field solution, and extract the local flow parameters on the aircraft surface, including local pressure, local density, local sound speed, and local flow velocity vector; interpolate these local flow parameters from the CFD mesh onto the aerodynamic surface element model, and establish the aerodynamic stiffness matrix in the modal coordinate system based on the local flow piston theory. ;

[0010] Step 3: Use the modal method to establish the static aeroelastic balance equation of the aircraft. Substitute the mode shape obtained in Step 1 and the aerodynamic stiffness matrix in the modal coordinate system established in Step 2 into the static aeroelastic balance equation to solve for the static aeroelastic deformation and balance amount of the aircraft.

[0011] Step 4: Based on the static aeroelastic deformation obtained in Step 3, update the CFD mesh of the aircraft to obtain the deformed CFD mesh; use CFD technology to perform a second steady flow field solution on the deformed CFD mesh to obtain the aerodynamic characteristics of the deformed aircraft.

[0012] In step 2 and step 4, only one steady flow field calculation based on CFD technology is performed each time, for a total of two steady flow field calculations based on CFD technology. In step 2, the local flow piston theory is only used to establish the aerodynamic stiffness matrix in the modal coordinate system to solve the static aeroelastic deformation of the aircraft in step 3. Furthermore, the local flow piston theory does not participate in the calculation of the aerodynamic characteristics of the deformed aircraft. The aerodynamic characteristics of the deformed aircraft are directly obtained from the second steady flow field solution in step 4.

[0013] Furthermore, in step 2, before interpolating the local flow parameters from the CFD mesh onto the aerodynamic surface element model, the expression for the pressure at each point on the surface of the aircraft when it is disturbed is first established based on the local flow piston theory:

[0014]

[0015] In the formula, Surface pressure; Local pressure; For local density; The washing speed is caused by the deformation of the object surface.

[0016] Furthermore, the characteristic is that the washing speed... Calculated using the following formula:

[0017]

[0018] In the formula, Let be the local flow velocity of the aerodynamic surface element on the object surface. This is the unit outward normal vector of the object surface in its undeformed state. It is the unit outward normal vector after the object surface has been deformed.

[0019] Furthermore, the aerodynamic stiffness matrix in the modal coordinate system established in step 2 satisfy:

[0020]

[0021] In the formula, For the overall modal aerodynamics of the aircraft; It is the modal aerodynamic force corresponding to the aerodynamic force in the undeformed state; It is the transpose of the centroidal mode shape matrix of the aerodynamic surface element; These are modal coordinates; It is the overall aerodynamic vector of the aircraft; It is the incoming flow density; It is the Mach number of the incoming flow; It is the incoming flow velocity;

[0022] Furthermore, the aerodynamic stiffness matrix in the modal coordinate system elements Calculated using the following formula:

[0023]

[0024] In the formula, This represents the total number of aerodynamic surface units; The number of the aerodynamic surface element indicates the number of the... One aerodynamic surface unit; , , and For the introduced local dimensionless parameters, where It is the ratio of the local pressure of the aerodynamic surface element on the object surface to the incoming flow pressure; It is the ratio of the local density of the aerodynamic surface element on the object surface to the incoming flow density; It is the ratio of the local sound velocity of the aerodynamic surface element on the object surface to the incoming flow sound velocity; It is the ratio of the local Mach number of the aerodynamic surface element on the object surface to the incoming Mach number; Indicates the first Unit outward normal vector of aerodynamic surface element under first-order mode shape; For the first Displacement of the first-order mode aerodynamic surface element; It is the area of ​​the aerodynamic surface unit.

[0025] Furthermore, step 2 also includes adjusting the aerodynamic stiffness matrix in the modal coordinate system based on the rigid body modes and elastic modes. The block-based step involves dividing the aerodynamic stiffness matrix in the modal coordinate system. The blocks are as follows:

[0026]

[0027] In the formula, Here is the aerodynamic stiffness matrix for rigid body modes. Here is the aerodynamic stiffness matrix under elastic modes. Let be the stiffness matrix of the rigid body modal aerodynamic forces induced by elastic modal displacements. is the stiffness matrix for the elastic aerodynamic forces induced by the rigid body modes.

[0028] Furthermore, in step 3, the static aeroelastic deformation of the aircraft is obtained by solving the following elastic deformation equation:

[0029]

[0030] In the formula, It is the structural stiffness matrix of the elastic modes; It is an elastic modal displacement; It is the transpose of the modal shape matrix; This is the partial derivative matrix of the distributed aerodynamic forces with respect to the balancing quantity at each structural node.

[0031] Furthermore, solving the static aeroelastic balancing equation in step 3 includes the following two scenarios:

[0032] Scenario 1: Given the angle of attack and rudder deflection, the elastic modal displacement can be obtained directly using the elastic deformation solution equation. Finally, the elastic deformation vectors of the structural nodes are used. With elastic modal displacement The relational expression yields the elastic deformation vector of the structural node;

[0033] Scenario 2: Given an acceleration trim amount, adjust the downwash velocity to cause a unit change in the aircraft's angle of attack or unit rudder deflection, thereby calculating the angle of attack and rudder deflection. Then, use the elastic deformation solution equation to obtain the elastic modal displacement. Finally, the elastic deformation vectors of the structural nodes are used. With elastic modal displacement The relational expression yields the elastic deformation vector of the structural node.

[0034] Furthermore, in step 4, the CFD mesh of the aircraft is updated using the radial basis function interpolation method, and a second steady flow field solution is performed using CFD technology to obtain the aerodynamic characteristics of the deformed aircraft. Specifically, this includes the following processes:

[0035] Step 4.1, Update the CFD mesh:

[0036] Step 1: Displacement boundary transfer;

[0037] displacement of the boundary nodes of the structured mesh Given the quantities, the RBF interpolation polynomial is constructed for the first time:

[0038]

[0039] In the formula, This represents the displacement of the boundary nodes of the structural mesh in the x, y, and z directions; Indicates the total number of boundary nodes of the structured mesh; These are the weighting coefficients. Indicates and The corresponding weighting coefficients; The function symbol for the radial basis functions; Represents the first in the structural boundary mesh The location of the radial basis function support point; Represents the first in the structural boundary mesh The location of the radial basis function support point; Represents the first in the structural boundary mesh The radial basis function support point to the first The distance between the support points of the radial basis functions;

[0040] The RBF interpolation polynomial constructed in the first step is rewritten as follows: In other words:

[0041] , , ;

[0042] in, This represents the displacement matrix of the boundary nodes of the structural mesh. Represents the radial basis function matrix;

[0043] The weight coefficient matrix is ​​obtained by solving the rewritten RBF interpolation polynomial. ;

[0044] Using the weight coefficient matrix And the coordinates of the CFD mesh boundary nodes, to calculate the CFD mesh boundary nodes. The displacement is obtained by measuring the displacement of the boundary nodes of the CFD mesh. :

[0045]

[0046] In the formula, Indicates the first The displacements of the boundary nodes of the CFD mesh in the x, y, and z directions; The total number of nodes at the CFD mesh boundary; Represents the boundary of the structured mesh. The weighting coefficients of the radial basis function support points in the x, y, and z directions; The function symbol for the radial basis functions; Represents the first in the CFD boundary mesh The location of the radial basis function support point; Represents the first in the structural boundary mesh The location of the radial basis function support point; Represents the first in the CFD mesh boundary The radial basis function support point is located in the structural boundary mesh. The distance between the support points of the radial basis functions;

[0047] Step 2: Update the internal grid;

[0048] The CFD spatial grid is updated based on the displacement information of the boundary nodes of the CFD grid, that is, the position of each boundary node in the CFD grid is updated.

[0049] Using the same transfer process as in the first step, where the displacement of the structural mesh boundary nodes is transferred to the CFD mesh boundary nodes, the displacement of the CFD mesh boundary nodes is transferred to all CFD mesh points in space via radial basis function interpolation, thus obtaining the displacement information of the nodes inside the CFD mesh. ;

[0050] Step 4.2: Using the displacement information of the boundary nodes of the CFD mesh and the displacement information of the internal nodes of the CFD mesh obtained in Step 4.1, update the node positions of the entire CFD mesh to obtain the deformed CFD mesh.

[0051] Step 4.3: Use CFD technology to perform a second steady flow field solution on the deformed CFD mesh to directly obtain the aerodynamic characteristics of the deformed aircraft, including lift coefficient, moment coefficient and pressure center position.

[0052] Furthermore, the radial basis function interpolation method uses Wendland's C2 function, which has the following form:

[0053]

[0054] In the formula, The function symbol representing the radial basis functions; , It is the position vector of any point in the computational domain. It is the first The position vector of the radial basis function support point. Represents any point in the computational domain and the first... The distance between the support points of the radial basis functions; Let be the radius of action of the radial basis functions; outside the radius of action, i.e. hour, .

[0055] Furthermore, the target aircraft mentioned in step 1 is a complex-shaped aircraft consisting of a rotating fuselage and control surfaces.

[0056] The advantages of this invention are:

[0057] 1. Compared with high-precision CFD methods, this invention only requires two CFD calculations, one for the undeformed state and one for the deformed state, saving the computational cost of multiple iterations and enabling rapid acquisition of static aeroelastic deformation and aerodynamic characteristics. Compared with the traditional local flow piston theory (LPT) method combined with CFD, this invention obtains more accurate aerodynamic characteristics through a second CFD calculation, significantly improving the calculation accuracy of the pressure core position and realizing the effective application of local flow piston theory in the static aeroelastic problems of complex-shaped aircraft.

[0058] 2. According to an embodiment of the present invention, under Mach 3 conditions, the computation time of the method of the present invention is reduced from 6 hours to 3.8 hours compared with the high-precision CFD method, a reduction of approximately 37%. For the aircraft provided in the embodiment of the present invention, under small angle of attack conditions of less than 10°, the pressure center movement calculated based on the method of the present invention maintains a good agreement with the high-precision CFD method, while the existing LPT method shows a significant deviation in pressure center movement under this condition. Therefore, the increased computation time of the method of the present invention compared with the existing LPT method results in a significant improvement in the accuracy of pressure center position calculation, achieving an optimized balance between accuracy and efficiency.

[0059] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0060] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0061] Figure 1 This is a flowchart of the static aeroelastic analysis method based on improved coupled CFD technology and local flow piston theory of the present invention;

[0062] Figure 2 This is a comparison diagram of fuselage deformation calculated by the method of the present invention and the high-precision CFD method at Mach number 3 and angle of attack 8° in an embodiment of the present invention.

[0063] Figure 3 This is a comparison chart of the changes in pressure core calculated using the LPT-C method of this invention, the existing LPT method, and the high-precision CFD method at different angles of attack. Detailed Implementation

[0064] The embodiments of the present invention are described in detail below. These embodiments are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0065] Reference Figure 1 The static aeroelastic analysis method based on improved coupled CFD technology and local flow piston theory provided in this embodiment includes the following steps:

[0066] Step 1: Model preparation and modal analysis.

[0067] First, for the target aircraft, a structural finite element model and an aerodynamic surface element model for piston theory calculations are established. The structural finite element model can be established using commercial finite element software such as Abaqus or Nastran; the aerodynamic surface element model is required to closely resemble the actual aerodynamic shape of the aircraft, therefore the positions of the aerodynamic surface elements in the model are often different from those of the structural elements.

[0068] Then, modal analysis was performed on the structural finite element model to obtain the mode shapes and frequencies of each mode of the aircraft. Let the mode shape matrix at the structural nodes be denoted as... It is divided into rigid body modal shape matrices as follows: And the elastic mode shape matrix is In this embodiment, all mode shape matrices are in mass-normalized form. Let the first... The frequency corresponding to the first mode shape is ;Pick If the model has 10 modalities, then the generalized mass matrix is: Generalized stiffness matrix is If only the elastic modes are considered, then the stiffness matrix of the elastic modes is denoted as... .

[0069] Because the structural unit and the aerodynamic surface unit are mismatched, it is necessary to obtain the mode shape at the centroid of the aerodynamic surface unit in order to achieve data transmission. In this embodiment, the mode shape at the centroid of the aerodynamic surface unit can be obtained from the mode shape matrix on the structural node. Derived through radial basis function interpolation, denoted as ,Will It is also divided into rigid body modes and elastic modes, that is ,in The rigid body mode shape matrix at the structural nodes. This represents the elastic mode shape matrix at the structural nodes.

[0070] Step 2: Establish the aerodynamic stiffness matrix in the modal coordinate system using the local flow piston theory combined with CFD technology.

[0071] Step 2.1: Create a CFD mesh for the undeformed state of the aircraft, considering an angle of attack of... For steady flight at Mach number Ma, the first steady flow field solution is obtained using CFD technology based on the Euler equations, yielding the flow field solution for the aircraft under steady flight conditions; local airflow parameters at various points on the extract surface, including local pressure, are also obtained. Local density Local speed of sound and local flow velocity vector Simultaneously, record the lift coefficient at this time. Torque coefficient .

[0072] Step 2.2: Interpolate the calculated local flow parameters onto the aerodynamic surface element model established in Step 1. First, when the aircraft is disturbed, establish the pressure at each point on the aircraft's surface according to the first-order piston theory. The expression for "can" is:

[0073] (1)

[0074] in, It is the washing speed, caused by the displacement of the object surface:

[0075] (2)

[0076] In the formula, Let be the local flow velocity of the aerodynamic surface element on the object surface. This is the unit outward normal vector of the object surface in its undeformed state. This is the unit outward normal vector after the object surface has undergone displacement. When the deformation is a small amount... Can be recorded as , which represents the angle between the external normal vectors before and after the deformation of the object surface.

[0077] Then, modal coordinates are introduced. Let the first modal coordinate be... The first modal coordinates are Before considering For the first mode, the total change in the outward normal direction of the object surface (i.e., the angle between the outward normal vectors before and after the deformation of the object surface) is:

[0078] (3)

[0079] In the formula, No. The unit outward normal vector of the object surface under first-order modal displacement.

[0080] The formula for calculating the surface pressure can be obtained from formulas (1) to (3):

[0081] (4)

[0082] Substitute the expression for the downward washing velocity into the pressure formula (4), and introduce the local dimensionless parameter. and After derivation, the modal coordinate expression of the pressure at various points on the surface of the aircraft can be obtained as follows:

[0083] (5)

[0084] In the formula, The incoming air pressure; The Mach number of the incoming flow; The speed of incoming sound; It is the incoming flow density.

[0085] Next, to establish the relationship between aerodynamic forces and modal coordinates, it is necessary to integrate the forces and moments. The object surface is divided into... There are discrete aerodynamic surface elements, each with an area of... The pressure on the aerodynamic surface element can be equivalent to a concentrated force acting on the centroid of the aerodynamic surface element. Therefore, based on the pressure on the aerodynamic surface element and the area of ​​the aerodynamic surface element, the force on each aerodynamic surface element in the physical coordinate system can be expressed as follows:

[0086] (6)

[0087] in, , , The components of the unit outward normal vector of the object surface in the x, y, and z directions of the global coordinate system; Let be the pressure of the nth aerodynamic surface element. Let n be the aerodynamic force of the nth aerodynamic surface element in the undeformed state. Let be the change in aerodynamic force on the nth aerodynamic surface element caused by deformation. By sequentially arranging and combining the forces acting on each aerodynamic surface element, the overall aerodynamic force vector of the aircraft can be obtained as follows:

[0088] (7)

[0089] In the formula, It is the overall aerodynamic vector of the aircraft; It refers to the overall aerodynamic force of the aircraft in its undeformed state; It is the change in the overall aerodynamic force of the aircraft caused by deformation.

[0090] Step 2.3: Transform the aerodynamic forces in the physical coordinate system to the modal coordinate system to obtain the aerodynamic stiffness matrix in the modal coordinate system. After the transformation, the overall modal aerodynamic forces of the aircraft can be expressed as:

[0091] (8)

[0092] in, It is the transpose of the centroidal mode shape matrix of the aerodynamic surface element; It is the incoming flow velocity; These are modal coordinates; For the overall modal aerodynamics of the aircraft; These are modal coordinates.

[0093] These are the modal aerodynamic forces corresponding to the aerodynamic forces in the undeformed state:

[0094] (9)

[0095] In the formula, The aerodynamic forces of the aircraft as a whole in its undeformed state. It is the rigid body mode shape of the centroid of the aerodynamic surface element; It is the elastic mode shape of the centroid of the aerodynamic surface element; It refers to the aerodynamic forces in the rigid body mode; It is the aerodynamic force in the elastic mode.

[0096] It is the aerodynamic stiffness matrix in modal coordinates, and its elements are... It can be calculated using the following formula:

[0097] (10)

[0098] In the formula, is the total number of aerodynamic surface elements; n is the number of the aerodynamic surface element, indicating the nth aerodynamic surface element; It is the ratio of the local pressure of each aerodynamic surface element to the incoming flow pressure. It is the ratio of the local density of the aerodynamic surface element on the object surface to the incoming flow density; It is the ratio of the local sound velocity of the aerodynamic surface element on the object surface to the incoming flow sound velocity; It is the ratio of the local Mach number of the aerodynamic surface element on the object surface to the incoming Mach number; Indicates the first The unit outward normal vector of an aerodynamic surface element under first-order mode shape; For the first Displacement of the centroid of the aerodynamic surface element in the first mode; It is the area of ​​the aerodynamic surface unit. This represents the unit outward normal vector of the object surface in its undeformed state.

[0099] aerodynamic stiffness matrix in modal coordinate system The system is divided into blocks based on rigid body modes and elastic modes, i.e. ,in Here is the aerodynamic stiffness matrix in rigid body modes. Here is the aerodynamic stiffness matrix under elastic modes. Let be the stiffness matrix of the rigid body modal aerodynamic forces induced by elastic modal displacements. is the stiffness matrix for the elastic aerodynamic forces induced by the rigid body modes.

[0100] Step 3: Solve the static aeroelastic balance equations of the aircraft based on the modal method.

[0101] The purpose of this step is to utilize the structural modes obtained in step 1 and the aerodynamic stiffness matrix in the modal coordinate system obtained in step 2. This study establishes and solves the static aeroelastic balance equations to obtain the rigid body motion, corresponding structural deformation, and aerodynamic forces required for the aircraft to maintain equilibrium considering structural elasticity; that is, the elastic deformation and rigid body balance of the aircraft. This embodiment uses a given angle of attack and rudder deflection angle, considering only the translational degree of freedom in the z-axis to maintain balance, as an example. Based on the principle of maintaining balance among the aerodynamic forces, inertial forces, and elastic forces at each structural node of an elastic aircraft, the following basic equations are established:

[0102] (11)

[0103] In the formula, For balancing angles, including angle of attack and rudder deflection; For each node, there is the elastic deformation vector; Here is the structural stiffness matrix; For the incoming flow pressure; This is the aerodynamic influence coefficient matrix; This is the partial derivative matrix of the distributed aerodynamic forces at each structural node with respect to the balancing angle. This invention only considers the translational degree of freedom in the z-axis, therefore the inertial force... It can be represented as:

[0104] (12)

[0105] In the formula, Let the total mass of the aircraft be the number of degrees of freedom considered. The distribution acceleration of each node, For the rigid body aircraft, the acceleration is given in the degrees of freedom. Let be the rigid body mode shape matrix for the considered degrees of freedom; It is the quality matrix.

[0106] Substituting equation (12) into equation (11), we get:

[0107] (13)

[0108] To transform the above equations into modal coordinates for solution, multiply them by the transpose of the elastic mode shape matrix on the left. The transpose of the rigid body modal shape matrix is ​​obtained, and the aerodynamic stiffness matrix in the elastic modal coordinates obtained in step 2 is used to divide the model into blocks. and Decoupling yields two key equations: one for solving elastic deformation and the other for solving rigid body balancing.

[0109] The derivation process of the elastic deformation solution equation is as follows:

[0110] Multiply equation (13) on the left by the transpose of the elastic modal shape matrix. According to the modal expansion theorem, the elastic deformation vector of the structural nodes can be determined. With elastic modal displacement The relationship is:

[0111] (14)

[0112] Furthermore, based on the orthogonality of the mode shape vectors with respect to the mass matrix and stiffness matrix... , And using the aerodynamic stiffness matrix in elastic modal coordinates Aerodynamic influence coefficient matrix Relationship:

[0113] (15)

[0114] Therefore, we can obtain:

[0115] (16)

[0116] In the formula, It is the structural stiffness matrix of the elastic modes.

[0117] Furthermore, elastic modal displacement It can be represented as:

[0118] (17)

[0119] The derivation process of the equations for balancing a rigid body is as follows:

[0120] Then multiply equation (13) on the left by the transpose of the rigid body mode shape matrix. And based on the orthogonality of the mode shape vectors with respect to the stiffness matrix and mass matrix... , And the rigid body aerodynamic forces caused by elastic modal displacements, determined by the aerodynamic stiffness matrix in the modal coordinate system. From the physical meaning, we can know:

[0121] (18)

[0122] Furthermore, it can be deduced that:

[0123] (19)

[0124] In the formula, It is the generalized mass matrix, which is the identity matrix.

[0125] Substituting formula (17) into formula (19) yields the equation for rigid body balancing. Based on this equation, the acceleration balancing amount can be calculated, which is the acceleration on the degrees of freedom considered for the rigid body aircraft. :

[0126] (20)

[0127] The solution to the above static aeroelastic balancing equations includes the following two scenarios:

[0128] Scenario 1: Solving for deformation given balancing amounts

[0129] Given trim values ​​such as angle of attack and rudder deflection, then and Given that the elastic modal displacement is obtained by directly using the elastic deformation solution equation shown in formula (17), the elastic modal displacement can be solved directly. Then, it is transformed into the elastic deformation of the structural nodes through equation (14). Meanwhile, the acceleration balance can be solved using the rigid body balance equation shown in formula (20), which is the acceleration on the degrees of freedom considered by the rigid body aircraft. .

[0130] Scenario 2: Finding the balancing quantity given acceleration

[0131] If the target acceleration is known, that is, the acceleration in the degrees of freedom considered by the rigid body aircraft. Then the washing speed in equation (1) can be changed. This causes the aircraft to exhibit a change in unit angle of attack or a change in unit rudder deflection, i.e.:

[0132] (twenty one)

[0133] Similarly:

[0134] (twenty two)

[0135] Therefore, the trim amount of the angle of attack or rudder deflection can be obtained by inverse solving formula (20). Then, based on the angle of attack and rudder deflection obtained from the inverse solution, the elastic modal displacement is solved using the elastic deformation solution equation shown in formula (17). Finally, it is transformed into structural node deformation through equation (14).

[0136] Step 4: Update the CFD mesh based on the static aeroelastic deformation, and then use CFD technology to analyze and obtain the aerodynamic characteristics of the aircraft under static aeroelastic deformation.

[0137] After obtaining the structural nodal deformations, the traditional method substitutes the deformations into the local flow piston theory expression to calculate the pressure changes on the surface, thereby calculating the total aerodynamic force. This is referred to below as the local flow piston theory method combined with CFD technology (LPT method). However, the LPT method has a large deviation in calculating the pressure center position and is inaccurate in calculating the aerodynamic characteristics after deformation. To improve the accuracy of aerodynamic characteristic analysis, the method of this invention uses the structural deformations obtained in step 3 to update the CFD mesh and uses CFD technology to perform a second flow field solution (LPT-C method) to obtain the accurate aerodynamic characteristics of the deformed aircraft, thus compensating for the insufficient accuracy of the traditional LPT method in calculating the aerodynamic forces after deformation.

[0138] First, based on the elastic deformation X of the structural nodes obtained in step 3, the CFD mesh of the aircraft is updated. Since the topology and node density of the CFD mesh boundary and the structural field mesh boundary are usually mismatched, the node displacement information required for mesh updating must be achieved through interpolation techniques. Radial Basis Function (RBF) interpolation, with its simple mathematical model, strong applicability, and independence from mesh topology, has been successfully applied in multi-field coupled numerical simulation problems. RBF interpolation is independent of the mesh and solver; the parameter values ​​at the target interpolation point depend only on the spatial location of that point, giving it advantages in complex meshes. Therefore, this invention uses the Radial Basis Function (RBF) interpolation method to achieve displacement information transfer at the interface between the structural mesh and the CFD mesh, and to handle the spatial mesh deformation of the CFD.

[0139] Regarding displacement information transfer at the interface between the structural mesh and the CFD mesh, the basic principle of the radial basis function interpolation method is to interpolate the parameters of the mesh nodes on one side of the coupling interface using radial basis functions, and then use the constructed radial basis function sequence to transfer the parameters to the mesh nodes on the other side of the coupling interface. Its basic form is as follows:

[0140] (twenty three)

[0141] In the formula, It is an interpolation function; This represents the total number of radial basis functions; The function symbol representing the radial basis functions; It is the position vector of any point in the computational domain; It is the first The position vector of the radial basis function support point; Represents any point in the computational domain and the first... The distance between the support points of the radial basis functions; Is with the first The weighting coefficients corresponding to the radial basis functions.

[0142] Preferably, the present invention uses Wendland's C2 function as the radial basis function, which has the following form:

[0143] (twenty four)

[0144] In the formula, , Let be the radius of action of the radial basis functions; outside the radius of action, i.e. hour, This radial basis function has advantages in mesh deformation interpolation, and its influence range is a circle in two-dimensional space and a sphere in three-dimensional space.

[0145] The CFD grid update process (i.e., the data transfer process) consists of two steps:

[0146] Step 1: Displacement Boundary Transfer

[0147] (1) First construction of RBF interpolation polynomial

[0148] displacement of the boundary nodes of the structured mesh Given the quantities, construct the RBF interpolation polynomial:

[0149] (25)

[0150] In the formula, This represents the displacement of the boundary nodes of the structural mesh in the x, y, and z directions; Indicates the total number of boundary nodes of the structured mesh; These are the weighting coefficients. Indicates and The corresponding weighting coefficients; The function symbol representing the radial basis functions; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The radial basis function support point to the first The distance between the radial basis function support points.

[0151] make:

[0152] (26)

[0153] (27)

[0154] (28)

[0155] in, This represents the displacement matrix of the boundary nodes of the structural mesh. Let represent the radial basis function matrix. Therefore, equation (25) can be simplified to:

[0156] (29)

[0157] The weight coefficient matrix can be obtained by solving formula (29). .

[0158] (2) Calculate the displacement values ​​of the boundary nodes of the CFD mesh.

[0159] Using the weight coefficient matrix obtained in the previous step, the CFD mesh boundary nodes are calculated based on their coordinates. Displacement:

[0160] (30)

[0161] In the formula, Indicates the first The displacements of the boundary nodes of the CFD mesh in the x, y, and z directions; The total number of nodes at the CFD mesh boundary; Represents the boundary of the structured mesh. The weighting coefficients of the radial basis function support points in the x, y, and z directions; The function symbol for the radial basis functions; Represents the first in the CFD boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the CFD mesh boundary The radial basis function support point is located in the structural boundary mesh. The distance between the radial basis function support points.

[0162] The displacement information of the boundary nodes of the CFD mesh can be calculated using formula (30). .

[0163] The second step is to update the internal mesh.

[0164] After obtaining the displacement information of the CFD mesh boundary nodes, the CFD spatial mesh needs to be updated based on this displacement information, i.e., the position of each boundary node in the CFD mesh needs to be updated. Similar to the process of transferring the displacement of the structural mesh boundary nodes to the CFD mesh boundary nodes in the first step, the displacement of the CFD mesh boundary nodes is again transferred to all points in the spatial CFD mesh using radial basis function interpolation (RBF). The specific process is as follows:

[0165] (1) Using the displacement of the CFD mesh boundary node obtained in the first step as a new known quantity, the RBF interpolation polynomial is reconstructed to obtain the weight coefficient matrix.

[0166] (31)

[0167] In the formula, This represents the displacement of the boundary nodes of the CFD mesh in the x, y, and z directions; The total number of nodes at the CFD mesh boundary; These are the weighting coefficients. Indicates and The corresponding weighting coefficients; The function symbol for the radial basis functions; Represents the first in the CFD boundary mesh The location of the radial basis function support points; Represents the first in the CFD boundary mesh The location of the radial basis function support points; Represents the first in the CFD mesh boundary The radial basis function support point to the first The distance between the radial basis function support points and their locations.

[0168] (2) Using the weighting coefficients obtained in the previous step The displacement values ​​of the nodes inside the CFD mesh are calculated based on the coordinates of the boundary nodes of the CFD mesh.

[0169] (32)

[0170] In the formula, Indicates the first The displacements of nodes inside the CFD mesh in the x, y, and z directions; This represents the total number of nodes within the CFD mesh. Represents the first CFD mesh boundary The weight coefficient corresponding to point number 1; The function symbol for the radial basis functions; Represents the first element in the CFD internal mesh. The location of the radial basis function support points; Represents the first in the CFD boundary mesh The location of the radial basis function support points; Represents the first element in the CFD internal mesh. The radial basis function support point is located in the CFD boundary mesh. The distance between the radial basis function support points.

[0171] The displacement information of the nodes inside the CFD mesh can be obtained using formula (32). .

[0172] After updating the node positions of the entire CFD mesh using the displacement information of all CFD mesh nodes obtained by interpolation, the deformed CFD mesh is obtained.

[0173] After obtaining the deformed CFD mesh, a second steady flow field solution was performed based on the Euler equations. Through this CFD calculation, the aerodynamic characteristics of the aircraft under static aeroelastic deformation state, including lift coefficient, moment coefficient, and pressure center position, can be directly obtained.

[0174] To verify the effectiveness of the method of the present invention (hereinafter referred to as the LPT-C method), the calculation results were compared with the high-precision CFD method and the traditional local flow piston theory method combined with CFD (LPT method).

[0175] Figure 2 The static aeroelastic deformation of the fuselage was calculated using both the traditional method combining CFD with the Local Flow Piston Theory (LPT method) and a high-precision CFD method. Figure 2 It can be seen that the static aeroelastic deformation obtained by the LPT method is very close to that obtained by the high-precision CFD method. Since the overall steps of the LPT method for calculating deformation are the same as those of the analysis method of the improved CFD and local flow piston theory coupling strategy of this invention (LPT-C method), the LPT-C method of this invention also has good accuracy in calculating static aerodynamic characteristic deformation.

[0176] The changes in the pressure center caused by static aeroelasticity have a significant impact on aircraft trim, so it is very important to accurately calculate the amount of movement of the pressure center. Figure 3 The results show that the displacement of the pressure core caused by static aeroelastic deformation at different angles of attack (Mach 3) was calculated using the LPT method, the LPT-C method of this invention, and the high-precision CFD method. The results indicate that when the angle of attack is above 10°, both the LPT method and the LPT-C method of this invention can obtain relatively accurate displacement of the pressure core. However, when the angle of attack is less than 10°, the displacement of the pressure core calculated by the LPT method is significantly higher than that calculated by the high-precision CFD method. Nevertheless, the displacement of the pressure core calculated by the LPT-C method of this invention still maintains a good agreement with that calculated by the high-precision CFD method.

[0177] Table 1 compares the computation time of the three methods. As can be seen from Table 1, compared with the high-precision CFD method, the computation time of the LPT-C method of this invention is reduced from 6 hours to 3.8 hours, a reduction of about 37%, while maintaining high core-pointing accuracy under small angle of attack conditions. Compared with the existing LPT method, although the LPT-C method of this invention increases the computation time slightly (by 1.6 hours), it achieves a significant improvement in core-pointing accuracy at small angles of attack, realizing an optimized balance between accuracy and efficiency.

[0178] Table 1 Comparison of calculation time using different methods

[0179]

[0180] The static aeroelasticity analysis method provided by this invention, based on an improved CFD and local flow piston theory coupling strategy, can be applied to the aeroelastic design and analysis of supersonic aircraft, and is particularly suitable for the static aeroelasticity analysis of complex-shaped aircraft such as fuselage-rudder assemblies. This invention significantly improves analysis efficiency while maintaining computational accuracy, and has industrial applicability.

[0181] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the scope of the technology disclosed in the present invention, and such modifications or substitutions should all be covered within the scope of protection of the present invention.

Claims

1. A static aeroelastic analysis method based on improved coupled CFD technology and local flow piston theory, characterized in that, Includes the following steps: Step 1: Establish the structural finite element model and aerodynamic surface element model of the target aircraft, perform modal analysis on the structural finite element model, and obtain the mode shapes of each order; Step 2: Establish a CFD mesh for the undeformed state of the aircraft, perform the first steady flow field solution, and extract the local flow parameters on the aircraft surface, including local pressure, local density, local sound speed, and local flow velocity vector; interpolate these local flow parameters from the CFD mesh onto the aerodynamic surface element model, and establish the aerodynamic stiffness matrix in the modal coordinate system based on the local flow piston theory. ; Step 3: Use the modal method to establish the static aeroelastic balance equation of the aircraft. Substitute the mode shape obtained in Step 1 and the aerodynamic stiffness matrix in the modal coordinate system established in Step 2 into the static aeroelastic balance equation to solve for the static aeroelastic deformation and balance amount of the aircraft. Step 4: Based on the static aeroelastic deformation obtained in Step 3, update the CFD mesh of the aircraft to obtain the deformed CFD mesh; use CFD technology to perform a second steady flow field solution on the deformed CFD mesh to obtain the aerodynamic characteristics of the deformed aircraft. In step 2 and step 4, only one steady flow field calculation based on CFD technology is performed each time, for a total of two steady flow field calculations based on CFD technology. In step 2, the local flow piston theory is only used to establish the aerodynamic stiffness matrix in the modal coordinate system to solve the static aeroelastic deformation of the aircraft in step 3. Furthermore, the local flow piston theory does not participate in the calculation of the aerodynamic characteristics of the deformed aircraft. The aerodynamic characteristics of the deformed aircraft are directly obtained from the second steady flow field solution in step 4.

2. The static aeroelasticity analysis method according to claim 1, characterized in that, In step 2, before interpolating the local flow parameters from the CFD mesh onto the aerodynamic surface element model, the expression for the pressure at each point on the aircraft surface when the aircraft is disturbed is first established based on the local flow piston theory: In the formula, Surface pressure; Local pressure; For local density; The washing speed is caused by the deformation of the object surface.

3. The static aeroelasticity analysis method according to claim 2, characterized in that, The washing speed Calculated using the following formula: In the formula, Let be the local flow velocity of the aerodynamic surface element on the object surface. This is the unit outward normal vector of the object surface in its undeformed state. It is the unit outward normal vector after the object surface has been deformed.

4. The static aeroelasticity analysis method according to claim 3, characterized in that, The aerodynamic stiffness matrix in the modal coordinate system established in step 2 satisfy: In the formula, For the overall modal aerodynamics of the aircraft; It is the modal aerodynamic force corresponding to the aerodynamic force in the undeformed state; It is the transpose of the centroidal mode shape matrix of the aerodynamic surface element; These are modal coordinates; It is the overall aerodynamic vector of the aircraft; It is the incoming flow density; It is the Mach number of the incoming flow; It is the incoming flow velocity; Furthermore, the aerodynamic stiffness matrix in the modal coordinate system elements Calculated using the following formula: In the formula, This represents the total number of aerodynamic surface units; The number of the aerodynamic surface element indicates the number of the... One aerodynamic surface unit; , , and For the introduced local dimensionless parameters, where It is the ratio of the local pressure of the aerodynamic surface element on the object surface to the incoming flow pressure; It is the ratio of the local density of the aerodynamic surface element on the object surface to the incoming flow density; It is the ratio of the local sound velocity of the aerodynamic surface element on the object surface to the incoming flow sound velocity; It is the ratio of the local Mach number of the aerodynamic surface element on the object surface to the incoming Mach number; Indicates the first Unit outward normal vector of aerodynamic surface element under first-order mode shape; For the first Displacement of the first-order mode aerodynamic surface element; It is the area of ​​the aerodynamic surface unit.

5. The static aeroelasticity analysis method according to claim 4, characterized in that, Step 2 also includes adjusting the aerodynamic stiffness matrix in the modal coordinate system based on the rigid body modes and elastic modes. The block-based step involves dividing the aerodynamic stiffness matrix in the modal coordinate system. The blocks are as follows: In the formula, Here is the aerodynamic stiffness matrix for rigid body modes. Here is the aerodynamic stiffness matrix under elastic modes. Let be the stiffness matrix of the rigid body modal aerodynamic forces induced by elastic modal displacements. is the stiffness matrix for the elastic aerodynamic forces induced by the rigid body modes.

6. The static aeroelasticity analysis method according to claim 5, characterized in that, In step 3, the static aeroelastic deformation of the aircraft is obtained by solving the following elastic deformation equation: In the formula, It is the structural stiffness matrix of the elastic modes; It is an elastic modal displacement; It is the transpose of the modal shape matrix; This is the partial derivative matrix of the distributed aerodynamic forces with respect to the balancing quantity at each structural node.

7. The static aeroelasticity analysis method according to claim 1, characterized in that, Step 3 involves solving the static aeroelastic balancing equation, which includes the following two scenarios: Scenario 1: Given the angle of attack and rudder deflection, the elastic modal displacement can be obtained directly using the elastic deformation solution equation. Finally, the elastic deformation vectors of the structural nodes are used. With elastic modal displacement The relational expression yields the elastic deformation vector of the structural node; Scenario 2: Given an acceleration trim amount, adjust the downwash velocity to cause a unit change in the aircraft's angle of attack or unit rudder deflection, thereby calculating the angle of attack and rudder deflection. Then, use the elastic deformation solution equation to obtain the elastic modal displacement. Finally, the elastic deformation vectors of the structural nodes are used. With elastic modal displacement The relational expression yields the elastic deformation vector of the structural node.

8. The static aeroelasticity analysis method according to claim 1, characterized in that, In step 4, the CFD mesh of the aircraft is updated using the radial basis function interpolation method, and the second steady flow field is solved using CFD technology to obtain the aerodynamic characteristics of the deformed aircraft. The specific process includes the following steps: Step 4.1, Update the CFD mesh: Step 1: Displacement boundary transfer; displacement of the boundary nodes of the structured mesh Given the quantities, the RBF interpolation polynomial is constructed for the first time: In the formula, This represents the displacement of the boundary nodes of the structural mesh in the x, y, and z directions; Indicates the total number of boundary nodes of the structured mesh; These are the weighting coefficients. Indicates and The corresponding weighting coefficients; The function symbol for the radial basis functions; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The radial basis function support point to the first The distance between the radial basis function support points; The RBF interpolation polynomial constructed in the first step is rewritten as follows: In other words: ， ， ； in, This represents the displacement matrix of the boundary nodes of the structural mesh. Represents the radial basis function matrix; The weight coefficient matrix is ​​obtained by solving the rewritten RBF interpolation polynomial. ; Using the weight coefficient matrix And the coordinates of the CFD mesh boundary nodes, to calculate the CFD mesh boundary nodes. The displacement is obtained by measuring the displacement of the boundary nodes of the CFD mesh. : In the formula, Indicates the first The displacements of the boundary nodes of the CFD mesh in the x, y, and z directions; The total number of nodes at the CFD mesh boundary; Represents the boundary of the structured mesh. The weighting coefficients of the radial basis function support points in the x, y, and z directions; The function symbol for the radial basis functions; Represents the first in the CFD boundary mesh The location of the radial basis function support points; Represents the first in the structural boundary mesh The location of the radial basis function support points; Represents the first in the CFD mesh boundary The radial basis function support point is located in the structural boundary mesh. The distance between the radial basis function support points; Step 2: Update the internal grid; The CFD spatial grid is updated based on the displacement information of the boundary nodes of the CFD grid, that is, the position of each boundary node in the CFD grid is updated. Using the same transfer process as in the first step, where the displacement of the structural mesh boundary nodes is transferred to the CFD mesh boundary nodes, the displacement of the CFD mesh boundary nodes is transferred to all CFD mesh points in space via radial basis function interpolation, thus obtaining the displacement information of the nodes inside the CFD mesh. ; Step 4.2: Using the displacement information of the boundary nodes of the CFD mesh and the displacement information of the internal nodes of the CFD mesh obtained in Step 4.1, update the node positions of the entire CFD mesh to obtain the deformed CFD mesh. Step 4.3: Use CFD technology to perform a second steady flow field solution on the deformed CFD mesh to directly obtain the aerodynamic characteristics of the deformed aircraft, including lift coefficient, moment coefficient and pressure center position.

9. The static aeroelasticity analysis method according to claim 8, characterized in that, The radial basis function interpolation method uses Wendland's C2 function, which has the following form: In the formula, The function symbol representing the radial basis functions; , It is the position vector of any point in the computational domain. It is the first The position vector of the radial basis function support point. Represents any point in the computational domain and the first... The distance between the support points of the radial basis functions; Let be the radius of action of the radial basis functions; outside the radius of action, i.e. hour, .

10. The static aeroelasticity analysis method according to claim 1, characterized in that, The target aircraft mentioned in step 1 is a complex-shaped aircraft consisting of a rotating fuselage and control surfaces.

Images