Helicopter blade dynamic load identification method based on blade profile bending moment measurement

Abstract

This paper proposes a new method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurements. Based on the force conditions of the blade unit, an inverse transfer matrix is ​​established. The measured blade bending moment is substituted into the matrix to calculate the blade displacement response. A dynamic model of the rotating blade is then established, and the blade's natural frequency and mode shape are calculated. The blade's generalized coordinates are calculated based on the blade's displacement response and modal displacement. Finally, the blade's aerodynamic loads are calculated by substituting these into the aerodynamic force identification equation derived from the blade's coupled motion equations.

Description

Technical Field

[0001] The invention belongs to the field of helicopter blade dynamic load identification, and specifically proposes a helicopter blade dynamic load identification method based on blade section bending moment measurement. Background Art

[0002] The main sources of helicopter vibration include rotor vibration loads, rotor mass and aerodynamic imbalance, engine and gearbox vibrations, and tail rotor vibration loads. Among them, rotor vibration loads and vibrations caused by rotor mass and aerodynamic imbalance are the main factors causing helicopter vibration. The rotor blades are in a complex, periodically changing aerodynamic environment. There are multiple aerodynamic, inertial, structural, and geometric couplings between the various degrees of freedom of the blades, and there is also a complex coupling relationship between the rotating rotor and the fuselage. This causes serious vibration problems for the helicopter, greatly limiting its application. Excessive vibration levels can affect ride comfort, increase the pilot's workload, reduce the structural fatigue life of the structure, reduce the reliability of the helicopter, and increase the use and maintenance costs of the helicopter.

[0003] The primary source of vibration in helicopters is rotor loads. The alternating loads on the rotor blades are the most important loads in determining the fatigue life of the rotor itself and other helicopter components. Rotor load identification determines blade structural loads and aerodynamic load distribution by measuring blade strain during flight. This provides a basis for various helicopter research areas, including aeroelastic analysis of helicopter blades, vibration reduction research, and component life prediction. Summary of the Invention

[0004] Purpose of the invention: In view of the shortcomings of the existing technology, the present invention proposes a method for identifying rotor dynamic loads based on the bending moment of each blade section, which has high accuracy and is easy to implement.

[0005] Technical solution: A rotor dynamic load identification method based on the bending moment of each blade section mainly includes the following steps:

[0006] Step 1: Discretize the rotor blade into several elastic blade units, use sensors to obtain the bending moment and torque of the rotor blade cross section, and perform Fourier series expansion on the bending moment and torque.

[0007] Step 2: Perform force balance analysis on the blade unit and derive the inverse transfer matrix, specifically:

[0008] (1) Force analysis in the swinging direction

[0009] Take the nth blade unit to analyze the force in the flapping direction. The force is as follows: Figure 2 As shown:

[0010] From the moment balance we can get:

[0011]

[0012] Where M y，n represents the bending moment on the right end of the nth blade unit; W represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q y，n It represents the inertia moment of the nth blade unit in the flapping direction.

[0013] From the relationship between the force on the blade unit and the rotation angle at both ends, we can get:

[0014]

[0015] Where, EI y，n represents the flapping stiffness of the nth blade unit; W′ n represents the deflection of the right side of the nth blade unit; represents the deflection of the left side of the nth blade unit;

[0016] From the relationship between the blade unit force and the deflection at both ends, we can get:

[0017]

[0018] (2) Force analysis in the swing direction

[0019] Take the nth blade unit to perform the force balance analysis in the swing direction, and the force is as follows: Figure 3 shown.

[0020] From the moment balance we can get:

[0021]

[0022] Where M z，n represents the bending moment on the right end of the nth blade unit; V represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q z，n It represents the inertia moment of the nth blade unit in the swing direction.

[0023] From the relationship between the force on the blade unit and the rotation angle at both ends, we can get:

[0024]

[0025] Where, EI z，n represents the shimmying stiffness of the nth blade unit; V′ n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit;

[0026] From the relationship between the force on the blade unit and the deflection at both ends, we can get:

[0027]

[0028] (3) Twist direction

[0029] Take the nth blade unit to perform torsional force balance analysis, the force is as follows Figure 4 shown.

[0030] From the unit torsion angle balance, we can get:

[0031]

[0032] Where k A Indicates the radius of gyration of the blade section; represents the torsion angle of the left end of the nth blade unit; φ n GJ represents the torsion angle of the right end of the nth blade unit; n represents the torsional stiffness of the nth blade unit; θ′ is the first-order derivative of the pitch angle with respect to x.

[0033] Step 3: Analyze the load relationship between adjacent elements of the blade and derive the conversion equation, which is:

[0034] Take the left side section of the nth blade unit and the right side section of the n+1th blade unit, and the angle between the two units is Δθ n , the force is as Figure 5 The load on the right side of the blade unit in segment n+1 is calculated based on the load on the left side of the blade unit in segment n and the angle between the units.

[0035] Step 4: Establish the inverse transfer matrix and calculate the blade displacement response, specifically:

[0036] Combining steps 2 and 3, a transfer matrix between adjacent units of the blade can be established, and then an inverse transfer matrix of the entire blade can be established.

[0037] The inverse transfer matrix equation can be expressed as:

[0038] [A n ]{X n}=[B n ]{X n+1}+{C n}

[0039] Among them, the profile state vector is

[0040]

[0041] Where, O c represents the amplitude of the cosine term of the state vector; s Represents the sine term of the state vector. [A n ] and [B n ] is a 14×14 matrix, {C n} is a 14×1 column vector.

[0042] After obtaining the distribution of blade bending moment and torque along the span direction of the blade, as long as the flapping angle, swing angle, blade root damping torque and displacement boundary conditions at the blade root are measured and substituted into the inverse transfer matrix, the state vector of each blade section can be calculated from the blade root to the blade tip, and the blade displacement response and shear force distribution along the span direction can be obtained.

[0043] Step 5: Calculate the natural modes of the rotating blades, specifically:

[0044] Discretize the blade into n units along the radial direction, and the Hamilton principle can be used to obtain:

[0045]

[0046] Where j represents the j-th blade unit; U represents strain energy; T represents kinetic energy; and W represents virtual work of external force.

[0047] The blade is divided into n segments, each segment is equivalent to a two-node blade element, and each node contains six degrees of freedom, namely tensile displacement u, flapping displacement w, flapping angle w′, shimmy displacement v, shimmy angle v′, and torsion angle φ. The unit node vector can be expressed as:

[0048]

[0049] Using the unit node vector and shape function to represent the displacement at the blade unit node, and substituting it into the energy expression, we can get:

[0050]

[0051] Where, [M] j represents the mass matrix of the jth unit; [C] j represents the damping matrix of the jth unit; [K] j represents the stiffness matrix of the jth element; [F] j represents the load of the jth element.

[0052] The total mass matrix and total stiffness matrix can be obtained by assembling each unit according to the relationship between the nodes of each unit.

[0053] In the case of free motion, the equation can be simplified to:

[0054]

[0055] The natural frequency matrix and natural mode matrix of the blade can be obtained according to the characteristic equation.

[0056] Step 6: Generalized coordinate calculation

[0057] The actual displacement q(r, t) of the blade section can be expressed as the product of each modal displacement φ(r, j) and the corresponding generalized coordinate g(t, j).

[0058]

[0059] Where r represents the section radius, and j represents the j-th mode.

[0060] Substituting the bending moment and torque after Fourier series expansion, the generalized coordinates can be obtained as:

[0061] [g]=([φ] T [φ]) -1 φ T [q]

[0062] Step 7: Aerodynamic load calculation

[0063] Under the action of alternating aerodynamic loads, the blades undergo forced vibration response. When damping is neglected, the corresponding equation for the blade forced vibration is:

[0064]

[0065] The actual displacement is expressed using modal displacement and generalized coordinates, and the Fourier series expansion is performed. After sorting, we can get:

[0066]

[0067]

[0068] Where m j is the rotating blade m j The modal mass, ω j is the j-order natural frequency of the rotating blade; are the j-order generalized coordinates k-order harmonic cosine component and k-order sine component, are the kth harmonic cosine component and kth sine component of the jth generalized aerodynamic load.

[0069] Beneficial effects: The present invention establishes an aerodynamic load identification model for rotor blades under the blade flapping and torsional coupling motion state. Compared with the existing technology, this method avoids the problem of ill-conditioned matrix, improves the calculation accuracy, and contains blade modal information, which is beneficial to subsequent aeroelastic research and vibration reduction research. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 This is a flow chart of a rotor dynamic load identification method based on the bending moment of each blade section;

[0071] Figure 2 This is the force analysis diagram of the blade unit in the flapping direction;

[0072] Figure 3 It is the force analysis diagram of the blade unit in the swing direction;

[0073] Figure 4 It is the force analysis diagram of the blade unit in the torsional direction;

[0074] Figure 5 It is a schematic diagram of the forces acting on adjacent units of the blade;

[0075] Figure 6 It is the aerodynamic load in the flapping direction of 0.75R;

[0076] Figure 7 It is 0.88R flapping direction aerodynamic load;

[0077] Figure 8 It is 0.97R aerodynamic load in the flapping direction. DETAILED DESCRIPTION

[0078] The present invention relates to a rotor dynamic load identification method based on the bending moment of each blade section, which mainly includes the following steps:

[0079] Step 1: Discretize the rotor blade into several elastic blade units, use sensors to obtain the bending moment and torque of the rotor blade cross section, and perform Fourier series expansion on the bending moment and torque.

[0080] Step 2: Perform force balance analysis on the blade unit and derive the inverse transfer matrix, specifically:

[0081] (1) Force analysis in the swinging direction

[0082] Take the nth blade unit to analyze the force in the flapping direction. The force is as follows: Figure 2 As shown:

[0083] From the moment balance we can get:

[0084]

[0085] Where M y，nrepresents the bending moment on the right end of the nth blade unit; W represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q y，n It represents the inertia moment of the nth blade unit in the flapping direction.

[0086] From the relationship between the force on the blade unit and the rotation angle at both ends, we can get:

[0087]

[0088] Where, EI y，n represents the flapping stiffness of the nth blade unit; W′ n represents the deflection of the right side of the nth blade unit; represents the deflection of the left side of the nth blade unit;

[0089] From the relationship between the blade unit force and the deflection at both ends, we can get:

[0090]

[0091] (2) Force analysis in the swing direction

[0092] Take the nth blade unit to perform the force balance analysis in the swing direction, and the force is as follows: Figure 3 shown.

[0093] From the moment balance we can get:

[0094]

[0095] Where M z，n represents the bending moment on the right end of the nth blade unit; V represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q z，n It represents the inertia moment of the nth blade unit in the swing direction.

[0096] From the relationship between the force on the blade unit and the rotation angle at both ends, we can get:

[0097]

[0098] Where, EI z，n represents the shimmying stiffness of the nth blade unit; V′ nrepresents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit;

[0099] From the relationship between the force on the blade unit and the deflection at both ends, we can get:

[0100]

[0101] (3) Twist direction

[0102] Take the nth blade unit to perform torsional force balance analysis, the force is as follows Figure 4 shown.

[0103] From the unit torsion angle balance, we can get:

[0104]

[0105] Where k A Indicates the radius of gyration of the blade section; represents the torsion angle of the left end of the nth blade unit; φ n GJ represents the torsion angle of the right end of the nth blade unit; n represents the torsional stiffness of the nth blade unit; θ′ is the first-order derivative of the pitch angle with respect to x.

[0106] Step 3: Analyze the load relationship between adjacent elements of the blade and derive the conversion equation, which is:

[0107] Take the left side section of the nth blade unit and the right side section of the n+1th blade unit, and the angle between the two units is Δθ n , the force is as Figure 5 The load on the right side of the blade unit in segment n+1 is calculated based on the load on the left side of the blade unit in segment n and the angle between the units.

[0108] Step 4: Establish the inverse transfer matrix and calculate the blade displacement response, specifically:

[0109] Combining steps 2 and 3, a transfer matrix between adjacent units of the blade can be established, and then an inverse transfer matrix of the entire blade can be established.

[0110] The inverse transfer matrix equation can be expressed as:

[0111] [A n ]{X n}=[B n ]{X n+1}+{C n} (8)

[0112] Among them, the profile state vector is

[0113]

[0114] Where, O c represents the amplitude of the cosine term of the state vector; s Represents the sine term of the state vector. [A n ] and [B n ] is a 14×14 matrix, {C n} is a 14×1 column vector.

[0115] After obtaining the distribution of blade bending moment and torque along the span direction of the blade, as long as the flapping angle, swing angle, blade root damping torque and displacement boundary conditions at the blade root are measured and substituted into the inverse transfer matrix, the state vector of each blade section can be calculated from the blade root to the blade tip, and the blade displacement response and shear force distribution along the span direction can be obtained.

[0116] Step 5: Calculate the natural modes of the rotating blades, specifically:

[0117] Discretize the blade into n units along the radial direction, and the Hamilton principle can be used to obtain:

[0118]

[0119] Where j represents the j-th blade unit; U represents strain energy; T represents kinetic energy; and W represents virtual work of external force.

[0120] The blade is divided into n segments, each segment is equivalent to a two-node blade element, and each node contains six degrees of freedom, namely tensile displacement u, flapping displacement w, flapping angle w′, shimmy displacement v, shimmy angle v′, and torsion angle φ. The unit node vector can be expressed as:

[0121]

[0122] Using the unit node vector and shape function to represent the displacement at the blade unit node, and substituting it into the energy expression, we can get:

[0123]

[0124] Where, [M] j represents the mass matrix of the jth unit; [C] j represents the damping matrix of the jth unit; [K] j represents the stiffness matrix of the jth element; [F] j represents the load of the jth element.

[0125] The total mass matrix and total stiffness matrix can be obtained by assembling each unit according to the relationship between the nodes of each unit.

[0126] In the case of free motion, the equation can be simplified to:

[0127]

[0128] The natural frequency matrix and natural mode matrix of the blade can be obtained according to the characteristic equation.

[0129] Step 6: Generalized coordinate calculation

[0130] The actual displacement q(r, t) of the blade section can be expressed as the product of each modal displacement φ(r, j) and the corresponding generalized coordinate g(t, j).

[0131]

[0132] Where r represents the section radius, and j represents the j-th mode.

[0133] Substituting the bending moment and torque after Fourier series expansion, the generalized coordinates can be obtained as:

[0134] [g]=([φ] T [φ]) -1 φ T [q] (15)

[0135] Step 7: Aerodynamic load calculation

[0136] Under the action of alternating aerodynamic loads, the blades undergo forced vibration response. When damping is neglected, the corresponding equation for the blade forced vibration is:

[0137]

[0138] The actual displacement is expressed using modal displacement and generalized coordinates, and the Fourier series expansion is performed. After sorting, we can get:

[0139]

[0140]

[0141] Where m j is the rotating blade m j The modal mass, ω j is the j-order natural frequency of the rotating blade; are the j-order generalized coordinates k-order harmonic cosine component and k-order sine component, is the kth harmonic cosine component and kth sine component of the jth generalized aerodynamic load. The synthesis can be used to calculate the aerodynamic loads along the blade span and azimuth.

[0142] Figure 6-Figure 8This chart compares the aerodynamic load test values ​​for a certain helicopter blade in the flapping direction with the identified values ​​calculated using the dynamic load identification method proposed in this invention, based on parameters such as the dynamic bending moment and dynamic torque of the helicopter blade. The proposed method improves calculation accuracy and includes blade modal information, facilitating subsequent aeroelastic and vibration reduction research.

Claims

1. A method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement, characterized in that: The following steps are involved: Step 1: Discretize the rotor blade into several elastic blade units, use sensors to obtain the bending moment and torque of the rotor blade section, and perform Fourier series expansion on the bending moment and torque; Step 2: Perform force balance analysis on the blade unit and derive equations for aerodynamic forces and section loads; Step 3: Analyze the load relationship between two adjacent blade units and derive the conversion equation; Step 4: Based on the equation obtained in step 2 and the conversion equation obtained in step 3, an inverse transfer matrix between adjacent elements of the blade is established, and the blade displacement response is obtained according to the inverse transfer matrix; Step 5: Calculate the natural modes of the rotor blades; Step 6: Calculate the generalized coordinates based on the blade displacement response obtained in step 4 and the rotor blade modal uniqueness obtained in step 5; Step 7: Calculate the aerodynamic load of the blade based on the rotor blade natural mode obtained in step 5 and the generalized coordinates obtained in step 6.

2. The method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement according to claim 1, characterized in that: The step 2 is specifically as follows: Step 2.1, perform a force analysis on the blade in the flapping direction. Take the nth blade unit and perform the force analysis in the flapping direction. From the bending moment balance, we can obtain: Where M y,n represents the bending moment on the right end of the nth blade unit; W represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q y，n It represents the inertia moment of the nth blade unit in the flapping direction; According to the relationship between the force on the blade unit and the rotation angles at both ends, we can get: Where, EI y，n represents the flapping stiffness of the nth blade unit; W′ n represents the deflection of the right side of the nth blade unit; represents the deflection of the left side of the nth blade unit; According to the relationship between the force on the blade unit and the deflection at both ends, we can get: Step 2.2: Force analysis in the swing direction. Take the nth blade unit and perform force balance analysis in the swing direction. According to the bending moment balance, we can get: Where M z，n represents the bending moment on the right end of the nth blade unit; V represents the bending moment on the left end of the nth blade unit; n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; l n Indicates the length of the nth blade unit; Q z，n It represents the inertia moment of the nth blade unit in the swing direction; According to the relationship between the force on the blade unit and the rotation angle at both ends, we can get: Where, EI z，n represents the shimmying stiffness of the nth blade unit; V′ n represents the deflection of the right end of the nth blade unit; represents the deflection of the left end of the nth blade unit; According to the relationship between the force on the blade unit and the deflection at both ends, we can get: Step 2.3, torsional direction analysis, take the nth blade unit and perform torsional direction force balance analysis, according to the unit torsional angle balance, we can get: Where k A Indicates the radius of gyration of the blade section; represents the torsion angle of the left end of the nth blade unit; φ n GJ represents the torsion angle of the right end of the nth blade unit; n represents the torsional stiffness of the nth blade unit; θ′ is the first-order derivative of the pitch angle with respect to x.

3. The method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement according to claim 2, characterized in that: The step three is specifically as follows: Set the left side section of the nth blade unit and the right side section of the n+1th blade unit, and the angle between the two units is Δθ n ; Calculate the load on the right side of the blade unit in the nth segment according to the load on the left side of the blade unit and the angle between the units; The conversion equation is derived based on the load relationship.

4. The method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement according to claim 2, characterized in that: The step 4 is specifically as follows: According to steps 2 and 3, a transfer matrix between two adjacent blade units is established, and then an inverse transfer matrix of the entire blade is established; The inverse transfer matrix equation is expressed as: [A n ]{X n }=[B n ]{X n+1 }+{C n } Among them, the profile state vector is Where, () c Represents the amplitude of the cosine term of the state vector; () s represents the sine term of the state vector, [A n ] and [B n ] is a 14×14 matrix, {C n } is a 14×1 column vector; After obtaining the distribution of blade bending moment and torque along the span direction of the blade, as long as the flapping angle, swing angle, blade root damping torque and displacement boundary conditions at the blade root are measured and substituted into the inverse transfer matrix, the state vector of each blade section can be calculated from the blade root to the blade tip, and the blade displacement response and shear force distribution along the span direction can be obtained.

5. The method for identifying dynamic loads of helicopter blades based on blade profile bending moment measurement according to claim 2, characterized in that: Step 5: Calculate the natural modes of the rotor blades. This is characterized by discretizing the blades into n segments along the radial direction. The Hamilton principle can be used to obtain: Where, j represents the jth blade unit; U represents strain energy; T represents kinetic energy; W represents the virtual work of external force; The blade is divided into n segments, each segment is equivalent to a two-node blade unit, and each node contains six degrees of freedom, namely, tensile displacement u, flapping displacement w, flapping angle w′, shimmy displacement v, shimmy angle v′, and torsion angle φ; the unit node vector can be expressed as: The displacement at the blade unit node is expressed by the unit node vector and the shape function, and substituted into the energy expression to obtain: Where, [M] j represents the mass matrix of the jth unit; [C] j represents the damping matrix of the jth unit; [K] j represents the stiffness matrix of the jth element; [F] j represents the load of the jth unit; The total mass matrix and the total stiffness matrix are obtained by assembling each unit according to the relationship between the nodes of each unit. In the case of free motion, the motion equation of the blade is simplified to: Where [M] is the total mass matrix; [K] is the total stiffness matrix; The natural frequency matrix and natural mode matrix of the blade can be obtained according to the characteristic equation.

6. The method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement according to claim 2, characterized in that: Step 6: Generalized coordinate calculation The actual displacement q(r, t) of the blade section can be expressed as the product of each modal displacement φ(r, j) and the corresponding generalized coordinate g(t, j); Where r represents the section radius, j represents the jth mode; Substituting the bending moment and torque after Fourier series expansion, the generalized coordinates can be obtained as: [g]=([φ] T [φ])-1φ T [q]。 7. The method for identifying dynamic loads on helicopter blades based on blade profile bending moment measurement according to claim 2, characterized in that: Step 7: Aerodynamic load calculation Under the action of alternating aerodynamic loads, the blade undergoes forced vibration response. When damping is neglected, the corresponding equation for the blade forced vibration is: The actual displacement is expressed using modal displacement and generalized coordinates, and the Fourier series expansion is performed. After sorting, we can get: Where m j is the rotor blade m j The modal mass, ω j is the j-order natural frequency of the rotor blade; are the j-order generalized coordinates k-order harmonic cosine component and k-order sine component, are the kth harmonic cosine component and kth sine component of the jth generalized aerodynamic load, and g is an element in the generalized coordinate matrix.

Images