Dynamic modeling method for attitude control system under dynamic fault

Abstract

The invention relates to a dynamic modeling method for an attitude control system under a dynamic fault. The method comprises the following steps: determining a modeling hypothesis of a nonlinear dynamic model; establishing a coordinate system used for modeling; rocket configuration parameters and engine related parameters are obtained; building displacement and acceleration descriptions of the three types of structures on the rocket; establishing a rigid body translation equation, a rigid body rotation equation, a rocket body elastic vibration equation, a propellant shaking equation and a sensitive element measurement equation; summarizing the five established equations to form a nonlinear dynamic equation set; and performing small-disturbance linearization expansion on the nonlinear dynamic equation set on the basis of the nominal reconstructed trajectory to obtain a linearized small-deviation equation set suitable for controller design and stability analysis. According to the method, the influence of the thrust drop fault is comprehensively considered, and the model can be used for carrying out relatively accurate dynamic response analysis and reconstruction scheme design under the fault.

Description

Technical Field

[0001] The invention belongs to the technical fields of overall design of launch vehicles, dynamics and control of launch vehicles, and relates to a dynamics modeling method for an attitude control system under power failure. Background Art

[0002] In order to meet the development needs of the aerospace industry, countries around the world are vigorously developing a new generation of high-thrust launch vehicles. Due to the limited thrust of a single engine, most bare rockets or bundled rockets currently use more than a dozen or even dozens of engines. The number of engines has increased significantly compared to traditional rockets, and the probability of failure, thrust reduction or even shutdown has greatly increased.

[0003] Thrust reduction will have many effects on the launch vehicle: (1) Thrust reduction will affect the accuracy of orbit insertion and even require trajectory reconstruction and orbit degradation; (2) Thrust reduction will cause the center of mass shift and overall parameter changes; (3) Thrust reduction will produce unbalanced interference torque; (4) The dynamic characteristics of liquid sloshing and the elastic mode of the entire rocket will change; (5) The control ability of the attitude control system will decrease, etc. Obviously, after the thrust reduction, the dynamic characteristics of the launch vehicle have deviated far from the nominal design conditions, and some of the assumptions used in the derivation of the traditional dynamic model have become invalid. The dynamic model derived based on normal conditions and the linearized small deviation model can no longer be used to describe the flight state under faults. Summary of the invention

[0004] The technical problem solved by the present invention is: to overcome the problem of inaccurate modeling of existing dynamic models under power system failures, and to propose a dynamic modeling method for an attitude control system under power failures, thereby achieving comprehensive consideration of the impact of thrust drop failures. The model can be used to carry out more accurate dynamic response analysis and reconstruction scheme design under failures.

[0005] The technical solution of the present invention is: a method for dynamic modeling of an attitude control system under power failure, comprising the following steps:

[0006] Establish the coordinate system used for modeling, including the launch coordinate system a, the rocket body coordinate system b, the semi-speed coordinate system h, the elastic coordinate system E, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks ;

[0007] Based on the launch vehicle object description, obtain the rocket configuration parameters and engine related parameters;

[0008] Based on the launch coordinate system a, the rocket body coordinate system b, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks, establish the displacement and acceleration description of three types of structures on the rocket, including the displacement and acceleration description of any point on the rocket body structure under the consideration of the elastic deformation of the rocket, the displacement and acceleration description of the swaying particle when the propellant sways, and the displacement and acceleration description of the center of mass of the nozzle of the kth booster sth engine;

[0009] Based on the semi-velocity coordinate system h, the rigid body translation equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0010] Based on the rocket body coordinate system b, the rigid body rotation equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0011] Based on the elastic coordinate system E, the elastic vibration equation of the rocket body is established according to the rocket configuration parameters and engine related parameters;

[0012] Based on the rocket body coordinate system b, the propellant sloshing equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0013] According to the output of the rigid body rotation equation and the mode at the installation position of the inertial device, the measurement equation of the sensitive element is established;

[0014] Summarize the rigid body translation equation, rigid body rotation equation, rocket body elastic vibration equation, propellant sloshing equation, and sensitive element measurement equation to form a nonlinear dynamic equation group;

[0015] The nonlinear dynamic equations are linearized with small perturbations based on the nominal reconstructed trajectory to obtain a linearized small deviation equation set suitable for controller design and stability analysis.

[0016] Furthermore, before conducting dynamic modeling, the modeling assumptions of the nonlinear dynamic model are clarified as follows:

[0017] Assume that the strain and rotation angle caused by the elastic deformation of the rocket body are 1% of the rocket body structure size;

[0018] Assume that the liquid surface sway angle is less than 10°;

[0019] Assume that the engine swing angle is less than 10°;

[0020] Ignore the effect of the Earth's rotation;

[0021] Only the reaction force generated by the ejection of propellant combustion products is considered and included in the thrust, while the additional Coriolis force, additional Coriolis moment and additional relative moment are ignored.

[0022] Furthermore, the arrow body coordinate system b is defined as follows:

[0023] The origin of the rocket body coordinate system is the projection point O1 of the rocket instantaneous center of mass O2 on the core stage longitudinal axis. The O1X1 axis is the core stage longitudinal axis, pointing to the core stage head. The O1Y1 axis is located in the main symmetry plane of the rocket. The O1X1Y1 plane coincides with the OXY plane of the launch coordinate system at the moment of launch. The O1Z1 axis is determined according to the right-hand rule.

[0024] The nozzle coordinate system N ks is defined as follows:

[0025] Nozzle coordinate system N ks Connected to the nozzle of the kth booster sth engine, origin O δks At this engine swing point, O δ ks X δks Along the symmetry axis of the nozzle shell, it points to the head of the arrow body, that is, the thrust P ks For a single pendulum motor, the positive direction of the servo direction angle is defined as O δks Y δks Axis, O δks Z δks The axis is determined according to the right-hand rule. For the orthogonal double-pendulum engine, after the positive polarity of the servo direction swing angle is determined, the unique O δks Y δks Axis and O δks Z δks The positive direction of the axis.

[0026] Furthermore, the rocket configuration parameters and engine related parameters include:

[0027] The rocket configuration parameters include: the number of boosters N zt , the total number of engines on the kth booster is The total number of tanks is When k = 0, it represents the core stage rocket body, that is, the total number of core stage rocket engines is The total number of tanks is Booster installation radius L k k represents the distance between the longitudinal axis of the booster and the longitudinal axis of the core stage, and the booster installation azimuth μ k represents the angle between the plane formed by the k-boost longitudinal axis and the core stage longitudinal axis and the O1X1Z1 plane;

[0028] The engine related parameters include: each module engine installation position X kR represents the distance between the swing point of the kth booster engine and the theoretical vertex; the installation radius L of the sth engine in the module δks , represents the distance between the engine swing point and the longitudinal axis of the module; δ qks The engine installation azimuth angle represents the angle between the line between the engine swing point and the module longitudinal axis and the O1Y1 axis; is the engine swing azimuth, indicating the rotation from the rocket body coordinate system O1Y1 axis to the nozzle coordinate system O δks Y δks angle.

[0029] Furthermore, the description of displacement and acceleration of any point on the rocket body structure includes:

[0030] The radius vector of any point on the rocket body structure:

[0031]

[0032] The acceleration of the structure particle is:

[0033]

[0034] In the formula, represents the absolute derivative in the emission coordinate system a, It represents the velocity of the origin O1 of the rocket body coordinate system relative to the launch coordinate system. represents the local derivative in the rocket body coordinate system b, In the local elastic moving coordinate system e k(x) The local derivative in u k(x) is the origin o of the local elastic moving coordinate system after the projectile is deformed ek(x) The elastic translation displacement, R c.m is the position vector of the origin O1 of the rocket body coordinate system in the launch coordinate system a, R k(x) is the radius vector from O1 to the projection point of any point of the structure on the longitudinal axis of the k-th booster before the k-th booster is deformed, ω1 is the angular velocity of the rigid body rotation of the rocket relative to the launch coordinate system a, ω ek(x) is the local kinematic elastic coordinate system e k(x) The elastic rotation angular velocity relative to the rocket body coordinate system b, In the coordinate plane o ek(x) y ek(x) z ek(x) origin o ek(x) The position radius vector to any structural point, is the position radius vector of any structural point when the rocket is regarded as a rigid body;

[0035] The description of the displacement and acceleration of the shaking particles when the propellant shakes includes:

[0036] When the rocket deforms and the propellant sloshes, the sloshing mass m kp The position vector in the launch coordinate system a can be expressed as:

[0037]

[0038] The acceleration of the shaking particle is:

[0039]

[0040] In the formula, is the radius of the sway center of the rocket before deformation in the rocket body coordinate system b, The sway center is located in the local elastic follower coordinate system e k(x) The translational displacement, h kp is the swaying mass in the local elastic moving coordinate system e k(x) The sway displacement;

[0041] The displacement and acceleration description of the center of mass of the nozzle of the k-th booster s-th engine includes:

[0042] In the launch coordinate system a, the position radius vector of the kth booster sth engine nozzle center of mass is:

[0043]

[0044] The nozzle center of mass acceleration is:

[0045]

[0046] In the formula, Swing center O for the kth booster s engine δks The position vector of the projection point on the boost longitudinal axis in the rocket body coordinate system b is: To establish the elastic translation displacement of the origin of the local elastic servo coordinate system at this projection point, ρ rks is the local elastic moving coordinate system e k(x) The swing point radius vector, ρ Rks is the nozzle coordinate system N ks The radius vector of the nozzle center of mass, ρ r1ks is the radius vector of the swing point when the rocket is regarded as a rigid body, ρ R1ks is the radius vector of the nozzle mass center when the engine is not swaying, is the local elastic follower coordinate system Angular velocity of rotation relative to the rocket body coordinate system b, ω δks is the nozzle relative to the local elastic coordinate system e k(x) Angular velocity of rotation.

[0047] Furthermore, the establishment of the rigid body translation equation includes:

[0048] According to Newton's second law, the rigid body translation equation of the entire rocket system is:

[0049]

[0050] In the formula, mg is gravity, P is engine thrust, F is aero is the aerodynamic force, F wind is the wind disturbance force, Fstructure is the structural disturbance force, is the volume fraction, For point Structure density at m kp is the sloshing mass of the p-th order sloshing of the k-th boost, m Rks is the swing mass of the k-th booster s-th engine, N zt is the number of boosters, is the total number of tanks, is the total number of engines on the kth booster;

[0051] In the half-velocity coordinate system h, the coefficientized rigid body translation equation is expanded as follows:

[0052]

[0053] At the same time, the apparent acceleration in the rocket body coordinate system is established as:

[0054]

[0055] Among them, g is the acceleration due to gravity, V is the flight speed, θ is the ballistic inclination, σ is the ballistic deviation angle, α is the ballistic attack angle, α wp ——steady wind angle of attack, α wq ——Shear wind angle of attack, β——Ballistic sideslip angle, β wp ——Side slip angle in steady wind, β wq ——shear wind sideslip angle, ——Rocket body pitch angular velocity, ——Yaw angular velocity of the rocket body, q i ——elastic generalized displacement, ——elastic generalized displacement velocity, δ Yks , are the swing angle and swing angular acceleration of the swing engine in the Y direction, δ Zks , are the swing angle and swing angular acceleration of the swing engine in the Z direction, c X_aero ——axial force coefficient, ——pitch channel aerodynamic coefficient, ——Yaw channel aerodynamic coefficient, ——pitch channel aerodynamic damping force coefficient, ——Yaw channel aerodynamic damping force coefficient, ——the i-th order elastic additional aerodynamic coefficient of the pitch channel, ——is the i-th order elastic additional aerodynamic coefficient of the yaw channel, ——the aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, ——The aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, c 1P ——axial thrust coefficient, are the swing control force coefficients of the pitch channel engine in the Y and Z directions, are the swing control force coefficients of the yaw channel engine in the Y and Z directions, are the swing inertia force coefficients of the pitch channel engine in the Y and Z directions, are the swing inertia force coefficients of the yaw channel engine in the Y and Z directions, ——Pitch channel sway cross-link force coefficient, ——Yaw channel sway cross-linking force coefficient, ——the i-th order elastic additional thrust coefficient of the pitch channel, ——the i-th order elastic additional thrust coefficient of the yaw channel, ——Interference force coefficient of pitch channel structure, ——Yaw channel structure disturbance force coefficient, ——Coordinate transformation matrix from rocket body coordinate system to half-speed coordinate system, is the projection of the sway acceleration on the O1Y1 axis of the rocket body coordinate system, is the projection of the sway acceleration on the O1Z1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1X1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Y1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Z1 axis of the rocket body coordinate system.

[0056] Furthermore, the establishment of the rigid body rotation equation includes:

[0057] The rigid body rotation equation of the entire rocket system is:

[0058]

[0059] Where M CM_g is the moment of gravity of each part on the origin O1 of the arrow body coordinate system, M CM_P is the torque of the engine thrust on the origin O1 of the rocket body coordinate system, M CM_aero is the moment of aerodynamic force on the origin O1 of the rocket coordinate system, M CM_wind is the moment of the wind disturbance force on the origin O1 of the arrow body coordinate system, M CM_structure is the moment of the structural disturbance force on the origin O1 of the rocket body coordinate system;

[0060] After expanding and substituting the external torque term, the coefficientized full rigid body rotation equation in the arrow body coordinate system is obtained:

[0061]

[0062]

[0063] In the formula, ——Rocket body rolling angular velocity, ——Roll channel yaw angle acceleration coupling coefficient, ——Roll channel pitch angle acceleration coupling coefficient, d PX ——Fault rolling channel mass center lateral displacement involvement coefficient, d1——Roll damping moment coefficient, d 3Yks ——k booster s engine Y-axis swing and roll control torque coefficient, d″ 3Yks ——Y-axis swing inertia moment coefficient of rolling channel, d 3Zks ——Z-axis swing and rolling control force coefficient, d″ 3Zks ——Z-axis swing inertia moment coefficient of rolling channel, d hykp —— Y-axis sloshing displacement of tank propellant to rolling channel cross-linking coefficient, d″ hykp ——Y-axis sway acceleration to rolling channel cross-linking coefficient, d hzkp ——Z-axis sway displacement to rolling channel cross-linking coefficient, d″ hzkp ——Z-axis sway acceleration to rolling channel cross-linking coefficient, d 1i ——The i-th order elastic aerodynamic damping moment coefficient of the rolling channel, d 2i ——the i-th order elastic moment coefficient of the rolling channel, ——interference torque coefficient of rolling channel structure, ——Yaw channel rolling angular acceleration coupling coefficient, ——Yaw channel pitch angle acceleration coupling coefficient, ——The lateral displacement coefficient of the center of mass of the faulty yaw channel, ——Yaw channel aerodynamic damping moment coefficient, ——Yaw channel aerodynamic moment coefficient, ——Yaw control torque coefficient of engine Y-axis swing, ——Y-axis swing inertia moment coefficient of the engine, ——Engine Z-axis swing yaw control torque coefficient, ——Inertia moment coefficient of engine Z-axis swing, ——Yaw channel sway moment cross-coefficient, ——Yaw channel sway center of mass cross-connection coefficient, ——The i-th order elastic additional aerodynamic damping moment coefficient of the yaw channel, ——The i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the yaw channel, ——engine thrust yaw channel moment, ——Yaw channel structure disturbance torque, ——Roll angular acceleration coupling coefficient of pitch channel, ——pitch channel yaw angle acceleration coupling coefficient, ——the lateral displacement coefficient of the center of mass of the fault pitch channel, ——pitch channel aerodynamic damping moment coefficient, ——pitch channel aerodynamic moment coefficient, ——Y-axis pitch control torque coefficient of the engine, The Y-axis pitch inertia moment coefficient of the engine, Engine Z-axis pitch control moment coefficient, ——engine Z-axis pitch inertia moment coefficient, ——the coupling coefficient of the pitch channel sway moment, ——pitch channel sway center of mass cross-connection coefficient, ——the i-th order elastic additional aerodynamic damping moment coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the pitch channel, ——engine thrust pitch channel moment, ——Pitch channel structure disturbance moment, Y kp ——Slosh normal displacement, Z kp ——Sway lateral displacement.

[0064] Furthermore, the establishment of the elastic vibration equation of the rocket body includes:

[0065]

[0066] In the formula, q i ——generalized displacement of the i-th order elastic vibration, ξ i ——damping ratio of the i-th order mode, ω i ——i-th order elastic frequency, D 4i ——Additional aerodynamic generalized force coefficient of rolling channel, ——the i-th order generalized aerodynamic damping force coefficient of the yaw channel, ——the i-th order generalized aerodynamic damping force coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the yaw channel, D 3iYks ——Generalized control force coefficient of engine Y-axis swing, D″ 3iYks ——Generalized inertia force coefficient of engine Y-axis swing, D 3iZks ——Generalized control force coefficient of engine Z-axis swing, D″ 3iZks ——Generalized inertia force coefficient of engine Z-axis swing, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing acceleration of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the yaw channel and the i-th order elastic coupling, ——Generalized force coefficient caused by the sloshing acceleration of the pth tank in the yaw channel and the i-th order elastic coupling, R ij ——generalized force coefficient of the j-th order elastic to i-th order elastic coupling, R′ ij ——generalized damping force coefficient of the j-th order elastic to i-th order elastic coupling, The generalized force generated by the thrust of the sth engine in the kth module is, Generalized force due to structural disturbance, Y kp ——Slosh normal displacement, Z kp —— lateral displacement of sway, ——sway normal acceleration, ——Sway lateral acceleration.

[0067] Furthermore, the establishment of the propellant sloshing equation includes:

[0068] The propellant sloshing equation in the rocket body coordinate system b is:

[0069]

[0070] are the sloshing displacement, velocity, and acceleration of the sloshing mass in the rocket body coordinate system, respectively. m kp is the propellant mass in case of thrust reduction failure, is the coupling coefficient matrix of rocket body rotation and liquid sloshing under thrust reduction failure, is the coupling term caused by the displacement of the center of mass of the rocket body due to the thrust reduction failure, is the equivalent damping coefficient matrix, is the equivalent stiffness coefficient matrix, is the apparent acceleration of the system, is the elastic motion-liquid sloshing coupling coefficient matrix, is the angular acceleration of the arrow body, is the elastic generalized acceleration;

[0071] Expand the above formula and derive the liquid sloshing equation as follows:

[0072]

[0073] ζ kpy , kpz are the normal and lateral sway damping ratios, Ω kp ——Sway frequency, Z Lkp——Z position of the p tank of the kth booster in the rocket body coordinate system, Y Lkp ——Y position of the p tank of the kth booster in the rocket body coordinate system, E kpz ——the distance from the center of mass to the center of mass, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the yaw channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the yaw channel and the p-th tank sloshing, ——Rocket body pitch acceleration, ——Yaw angular acceleration of the rocket body, ——Rolling angular acceleration of the rocket body, n——Number of non-zero natural frequencies.

[0074] Furthermore, the establishment of the sensitive element measurement equation includes:

[0075] The inertial measurement equation is:

[0076]

[0077] In the formula, φ 0ix (x gz ) is the torsional vibration mode of the core stage rocket, φ 0iy (x gz ) and φ 0ix (x gz ) is the vibration mode slope at the core stage rocket, ——Inertial measurement value of pitch attitude angle considering elasticity, ——Rocket body pitch attitude angle, ψ gz ——Yaw attitude angle measured by the inertial group considering elasticity, ψ——Yaw attitude angle of the rocket body, γ gz ——Measured value of the rocket body rolling attitude angle considering elasticity, γ——rocket body rolling attitude angle;

[0078] Assuming that the rate gyro is installed on the kth booster, its measurement equation is:

[0079]

[0080] In the formula, φ kix (x st ) is the torsional vibration mode of the kth booster rate gyro installation, φ kiy (x st ) and φ kiz (x st ) is the vibration mode slope at the kth booster upper rate gyro installation location, ωZ1_st ——Rocket body pitch angular velocity rate gyro measurement value, ——Rocket body pitch angular velocity, ——Measurement value of the rocket body yaw angular velocity rate gyro, ——Yaw angular velocity of the rocket body, ——Rocket body rolling angular velocity rate gyro measurement value, ——Rocket body rolling angular velocity, ——Elastic generalized velocity.

[0081] Furthermore, the linearized small deviation equation group includes:

[0082] 1) Small deviation equation of center of mass motion:

[0083]

[0084] ——dynamic coefficient of pitch channel, ——Yaw channel dynamic coefficient, ——gravity coefficient of pitch channel, ——Yaw channel gravity coefficient, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the yaw channel, Δ represents the deviation perturbation of the state quantity based on the nominal reconstructed trajectory;

[0085] 2) Small deviation equation of rigid body rotation:

[0086]

[0087]

[0088] 3) Linearization results of elastic motion equations:

[0089]

[0090] 4) The small deviation equation of shaking is:

[0091]

[0092] The beneficial effects of the present invention compared with the prior art are:

[0093] The present invention establishes a nonlinear dynamic model of a large-scale bundled launch vehicle that considers rigid-sloshing coupling under any multiple engine thrust reduction failures, and the dynamic model includes a rigid body translation equation, a rigid body rotation equation, a full rocket elastic vibration equation, a liquid sloshing equation, and a sensitive element measurement equation. The established nonlinear dynamic model has a certain universality, and can not only reflect the rocket mass eccentricity, the change of the moment of inertia and the product of inertia, the unbalanced interference torque caused by thrust asymmetry, the center of mass and the coupling characteristics of the motion around the center of mass that may be caused by any multiple engine thrust reductions. It is also applicable to normal flight bundled launch vehicles with general any multiple boosters and asymmetric distribution of engines. BRIEF DESCRIPTION OF THE DRAWINGS

[0094] Figure 1 This is a flow chart of a method for dynamic modeling of an attitude control system under power failure;

[0095] Figure 2 Define schematic diagrams for the rocket body coordinate system and the elastic coordinate system;

[0096] Figure 3 (a) is a schematic diagram of the engine installation orientation. Figure 3 (b) Schematic diagram of the engine swing orientation;

[0097] Figure 4 Schematic diagram of radial vectors at different mass point positions;

[0098] Figure 5 This is the force analysis diagram of the interface between the center body and the substructure; DETAILED DESCRIPTION

[0099] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0100] The present invention proposes a method for dynamic modeling of an attitude control system under power failure, such as Figure 1 As shown, the following steps are included:

[0101] (1) Clarify the modeling assumptions of the nonlinear dynamics model;

[0102] (2) Establish the coordinate systems used for modeling, including the launch coordinate system a, the rocket body coordinate system b, the semi-velocity coordinate system h, the elastic coordinate system E, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks ;

[0103] (3) Based on the launch vehicle object description, obtain the rocket configuration parameters and engine related parameters;

[0104] (4) Based on the launch coordinate system a, the rocket body coordinate system b, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks, establish the displacement and acceleration description of three types of structures on the rocket, including the displacement and acceleration description of any point on the rocket body structure under the consideration of the elastic deformation of the rocket, the displacement and acceleration description of the swaying particle when the propellant sways, and the displacement and acceleration description of the center of mass of the nozzle of the kth booster sth engine;

[0105] (5) Based on the semi-velocity coordinate system h, the rigid body translation equation is established according to the rocket configuration parameters and engine-related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0106] (6) Based on the rocket body coordinate system b, the rigid body rotation equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0107] (7) Based on the elastic coordinate system E, the elastic vibration equation of the rocket body is established according to the rocket configuration parameters and engine related parameters;

[0108] (8) Based on the rocket body coordinate system b, the propellant sloshing equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket;

[0109] (9) Establishing the measurement equation of the sensitive element based on the output of the rigid body rotation equation and the mode at the installation position of the inertial device;

[0110] (10) Summarize the rigid body translation equation, rigid body rotation equation, rocket body elastic vibration equation, propellant sloshing equation, and sensitive element measurement equation to form a nonlinear dynamic equation group;

[0111] (11) The nonlinear dynamic equations are linearized with small perturbations based on the nominal reconstructed trajectory to obtain a linearized small deviation equation set suitable for controller design and stability analysis.

[0112] Example 1

[0113] Step 1: Clarify the modeling assumptions of the nonlinear dynamic model.

[0114] When deriving the nonlinear dynamics model, the following assumptions are made:

[0115] (1) Assuming that the elastic deformation of the rocket body is a small deformation, the strain and rotation angle caused by the elastic deformation are very small relative to the original size of the rocket body, about 1% of the size of the rocket body structure;

[0116] (2) Assuming that the liquid sloshes with a small amplitude, it is generally believed that the assumption is met when the liquid surface sloshing angle is less than 10°;

[0117] (3) Assume that the rocket body rotates slowly around its center of mass (the angular velocity ω1 is small, and the angular velocity of a general launch vehicle satisfies the assumed conditions);

[0118] (4) Assuming that the engine swing angle is small, it is considered that the assumption is satisfied when the swing angle is less than 10°;

[0119] (5) Ignore the effect of the Earth's rotation;

[0120] (6) Only the reaction force generated by the ejection of propellant combustion products is considered and included in the thrust, while the additional Coriolis force, additional Coriolis moment and additional relative moment generated are ignored.

[0121] Step 2: Establish the coordinate system used for modeling.

[0122] The following coordinate systems are used in the modeling process:

[0123] (1) Emission coordinate system a: O-XYZ;

[0124] (2) Rocket body coordinate system b: O1-X1Y1Z1;

[0125] (3) Velocity coordinate system c: O2-X c Y c Z c ;

[0126] (4) Half-speed coordinate system h: O2-X h Y h Z h ;

[0127] (5) Elastic coordinate system E:O E -X E Y E Z E ;

[0128] (6) Local elastic moving coordinate system e k(x) :O ek(x) -X ek(x) Y ek(x) Z ek(x) ;

[0129] (7) Nozzle coordinate system N ks :O δks -X δks Y δks Z δks .

[0130] The coordinate system defines the arrow body coordinate system b and the nozzle coordinate system N ks The definitions of other coordinate systems are the same as the traditional definitions. The rocket body coordinate system b and the nozzle coordinate system N ks The detailed definition is as follows:

[0131] Rocket body coordinate system b: The origin is the projection point O1 of the rocket instantaneous center of mass O2 (referring to the center of mass of the rocket body without considering liquid sloshing and engine swing) on ​​the longitudinal axis of the core stage (note that points O1 and O2 do not necessarily coincide, such as Figure 2 As shown in the figure, O1X1 axis is the longitudinal axis of the core stage, pointing to the core stage head, O1Y1 axis is located in the main symmetry plane of the rocket, O1X1Y1 plane coincides with the OXY plane of the launch coordinate system at the moment of launch, and O1Z1 axis is determined by the right-hand rule.

[0132] Nozzle coordinate system N ks : Nozzle coordinate system N ks Connected to the nozzle of the kth booster sth engine, origin O δks At this engine swing point, O δks X δks Along the symmetry axis of the nozzle shell, it points to the head of the arrow body, that is, the thrust P ks For a single pendulum motor, the positive direction of the servo direction angle is defined as O δks Y δks Axis, O δks Z δks For an orthogonal double-pendulum engine, after the positive polarity of the servo direction swing angle is determined, the unique O can be determined according to the right-hand rule. δks Y δks Axis and O δks Z δks The positive direction of the axis. Figure 3 (a) Figure 3 (b) as shown.

[0133] Step 3: Based on the launch vehicle object description, obtain the rocket configuration parameters and engine related parameters.

[0134] The rocket configuration parameters include: the number of boosters N zt , the total number of engines on the kth booster is The total number of tanks is In particular, when k = 0, it represents the core stage rocket body, that is, the total number of core stage rocket engines is The total number of tanks is Booster installation radius L k k represents the distance between the longitudinal axis of the booster and the longitudinal axis of the core stage, and the booster installation azimuth μ k It represents the angle between the plane formed by the longitudinal axis of the k-boost and the core stage and the O1X1Z1 plane. Taking the tail view as a reference, the clockwise rotation of the O1X1Z1 plane is positive.

[0135] The engine related parameters include: each module engine installation position X kR represents the distance between the swing point of the kth booster engine and the theoretical vertex; the installation radius L of the sth engine in the module δks , i.e. the distance between the engine swing point and the longitudinal axis of the module; δ qksThe engine installation azimuth angle describes the angle between the line between the engine swing point and the module longitudinal axis and the O1Y1 axis. With the tail view as a reference, O1Y1 rotates clockwise as positive; is the engine swing azimuth, indicating the rotation from the rocket body coordinate system O1Y1 axis to the nozzle coordinate system O δks Y δks The angle is measured with the tail view as reference and clockwise is positive.

[0136] Step 4: Based on the launch coordinate system a, the rocket body coordinate system b, and the local elastic follower coordinate system e k(x) , establish the displacement and acceleration description of the three types of structures on the arrow. See the schematic diagram of the radial vector at different mass point positions. Figure 4 .

[0137] 1) Description of displacement and acceleration of any point on the rocket structure:

[0138] The radius vector of any point on the rocket body structure:

[0139]

[0140] The acceleration of the structure particle is:

[0141]

[0142] In the formula, represents the absolute derivative in the emission coordinate system a, It represents the velocity of the origin O1 of the rocket body coordinate system relative to the launch coordinate system, denoted by V1 below. represents the local derivative in the rocket body coordinate system b, In the local elastic moving coordinate system e k(x) The local derivative in u k(x) is the origin o of the local elastic moving coordinate system after the projectile is deformed ek(x) The elastic translation displacement, R c.m is the position vector of the origin O1 of the rocket body coordinate system in the launch coordinate system a, R k(x) is the radius vector from O1 to the projection point of any point of the structure on the longitudinal axis of the k-th booster before the k-th booster is deformed, ω1 is the angular velocity of the rigid body rotation of the rocket relative to the launch coordinate system a, ω ek(x) is the local kinematic elastic coordinate system e k(x) The elastic rotation angular velocity relative to the rocket body coordinate system b, Represented in the coordinate plane origin o ek(x) The position radius vector to any structural point, is the position radius vector of any structural point when the rocket is regarded as a rigid body.

[0143] 2) Description of the displacement and acceleration of the shaking particle when the rocket is deformed and the propellant is shaking:

[0144] When the rocket deforms and the propellant sloshes, the sloshing mass m kp The position vector in the launch coordinate system a can be expressed as:

[0145]

[0146] The acceleration of the shaking particle is

[0147]

[0148] In the formula, is the radius of the sway center of the rocket before deformation in the rocket body coordinate system b, The sway center is located in the local elastic follower coordinate system e k(x) The translational displacement, h kp is the swaying mass in the local elastic moving coordinate system e k(x) The shaking displacement.

[0149] 3) Description of the displacement and acceleration of the center of mass of the nozzle of the kth booster sth engine:

[0150] In the launch coordinate system a, the position radius vector of the kth booster sth engine nozzle center of mass is:

[0151]

[0152] The nozzle center of mass acceleration is:

[0153]

[0154] In the formula, Swing center O for the kth booster s engine δks The position vector of the projection point on the boost longitudinal axis in the rocket body coordinate system b is: To establish the elastic translation displacement of the origin of the local elastic servo coordinate system at this projection point, ρ rks is the local elastic moving coordinate system e k(x) The swing point radius vector, ρ Rks is the nozzle coordinate system N ks The radius vector of the nozzle center of mass, ρ r1ks is the radius vector of the swing point when the rocket is regarded as a rigid body, ρ R1ks is the radius vector of the nozzle mass center when the engine is not swaying, is the local elastic follower coordinate system Angular velocity of rotation relative to the rocket body coordinate system b, ω δks is the nozzle relative to the local elastic coordinate system e k(x) Angular velocity of rotation.

[0155] Step 5: Based on the semi-velocity coordinate system h, the rigid body translation equation is established according to the rocket configuration parameters, engine-related parameters, and the displacement and acceleration descriptions of the three types of structures on the rocket.

[0156] According to Newton's second law, the rigid body translation equation of the entire rocket system is:

[0157]

[0158] In the formula, mg is gravity, P is engine thrust (including control force), F aero is the aerodynamic force (including aerodynamic force and aerodynamic damping force), F wind is the wind disturbance force, F structure is the structural disturbance force, is the volume fraction, For point Structure density at m kp is the sloshing mass of the p-th order sloshing of the k-th boost, m Rks is the swing mass of the sth engine in the kth boost, N zt is the number of boosters, is the total number of tanks, is the total number of engines on the kth booster.

[0159] In the half-velocity coordinate system h, the coefficientized rigid body translation equation is expanded as follows:

[0160]

[0161] At the same time, the apparent acceleration in the rocket body coordinate system is established as:

[0162]

[0163] Among them, g is the acceleration due to gravity, V is the flight speed, θ is the ballistic inclination, σ is the ballistic deviation angle, α is the ballistic attack angle, α wp ——steady wind angle of attack, α wq ——Shear wind angle of attack, β——Ballistic sideslip angle, β wp ——Side slip angle in steady wind, β wq ——shear wind sideslip angle, ——Rocket body pitch angular velocity, ω Y1 ——Yaw angular velocity of the rocket body, q i ——elastic generalized displacement, ——elastic generalized displacement velocity, δ Yks , are the swing angle and swing angular acceleration of the swing engine in the Y direction, δ Zks , are the swing angle and swing angular acceleration of the swing engine in the Z direction, cX_aero ——axial force coefficient, ——pitch channel aerodynamic coefficient, ——Yaw channel aerodynamic coefficient, ——pitch channel aerodynamic damping force coefficient, ——Yaw channel aerodynamic damping force coefficient, ——the i-th order elastic additional aerodynamic coefficient of the pitch channel, ——is the i-th order elastic additional aerodynamic coefficient of the yaw channel, ——the aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, ——The aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, c 1P ——axial thrust coefficient, are the swing control force coefficients of the pitch channel engine in the Y and Z directions, are the swing control force coefficients of the yaw channel engine in the Y and Z directions, are the swing inertia force coefficients of the pitch channel engine in the Y and Z directions, are the swing inertia force coefficients of the yaw channel engine in the Y and Z directions, ——Pitch channel sway cross-link force coefficient, ——Yaw channel sway cross-linking force coefficient, ——the i-th order elastic additional thrust coefficient of the pitch channel, ——the i-th order elastic additional thrust coefficient of the yaw channel, ——Interference force coefficient of pitch channel structure, ——Yaw channel structure disturbance force coefficient, ——Coordinate transformation matrix from rocket body coordinate system to half-speed coordinate system, is the projection of the sway acceleration on the O1Y1 axis of the rocket body coordinate system, is the projection of the sway acceleration on the O1Z1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1X1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Y1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Z1 axis of the rocket body coordinate system.

[0164] Step 6: Based on the rocket body coordinate system b, establish the rigid body rotation equation according to the rocket configuration parameters, engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket.

[0165] According to D'Alembert's principle, the sum of the moments of the resultant force on the launch vehicle about the origin O1 of the rocket body coordinate system is zero. The rigid body rotation equation of the entire rocket system is:

[0166]

[0167] Where M CM_g is the moment of gravity of each part on the origin O1 of the arrow body coordinate system, M CM_P M is the torque of the engine thrust (including control force) on the origin O1 of the rocket body coordinate system, CM_aero is the moment of aerodynamic force (including aerodynamic moment and aerodynamic damping moment) on the origin O1 of the rocket body coordinate system, M CM_wind is the moment of the wind disturbance force on the origin O1 of the arrow body coordinate system, M CM_structure is the moment of the structural disturbance force on the origin O1 of the rocket body coordinate system.

[0168] After expanding and substituting the external torque term, the coefficientized full rigid body rotation equation in the arrow body coordinate system can be obtained:

[0169]

[0170]

[0171] In the formula, ——Rocket body rolling angular velocity, ——Roll channel yaw angle acceleration coupling coefficient, ——Roll channel pitch angle acceleration coupling coefficient, d PX ——Fault rolling channel mass center lateral displacement involvement coefficient, d1——Roll damping moment coefficient, d 3Yks ——k booster s engine Y-axis swing and roll control torque coefficient, d″ 3Yks ——Y-axis swing inertia moment coefficient of rolling channel, d 3Zks ——Z-axis swing and rolling control force coefficient, d″ 3Zks ——Z-axis swing inertia moment coefficient of rolling channel, d hykp —— Y-axis sloshing displacement of tank propellant to rolling channel cross-linking coefficient, d″ hykp ——Y-axis sway acceleration to rolling channel cross-linking coefficient, d hzkp ——Z-axis sway displacement to rolling channel cross-linking coefficient, d″ hzkp ——Z-axis sway acceleration to rolling channel cross-linking coefficient, d 1i ——The i-th order elastic aerodynamic damping moment coefficient of the rolling channel, d 2i ——the i-th order elastic moment coefficient of the rolling channel, ——interference torque coefficient of rolling channel structure, ——Yaw channel rolling angular acceleration coupling coefficient, ——Yaw channel pitch angle acceleration coupling coefficient, ——The lateral displacement coefficient of the center of mass of the faulty yaw channel, ——Yaw channel aerodynamic damping moment coefficient, ——Yaw channel aerodynamic moment coefficient, ——Yaw control torque coefficient of engine Y-axis swing, ——Y-axis swing inertia moment coefficient of the engine, ——Engine Z-axis swing yaw control torque coefficient, ——Inertia moment coefficient of engine Z-axis swing, ——Yaw channel sway moment cross-coefficient, ——Yaw channel sway center of mass cross-connection coefficient, ——The i-th order elastic additional aerodynamic damping moment coefficient of the yaw channel, ——The i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the yaw channel, ——engine thrust yaw channel moment, ——Yaw channel structure disturbance torque, ——Roll angular acceleration coupling coefficient of pitch channel, ——pitch channel yaw angle acceleration coupling coefficient, ——the lateral displacement coefficient of the center of mass of the fault pitch channel, ——pitch channel aerodynamic damping moment coefficient, ——pitch channel aerodynamic moment coefficient, ——Y-axis pitch control torque coefficient of the engine, The engine Y-axis pitch inertia moment coefficient, Engine Z-axis pitch control moment coefficient, ——engine Z-axis pitch inertia moment coefficient, ——the coupling coefficient of the pitch channel sway moment, ——pitch channel sway center of mass cross-connection coefficient, ——the i-th order elastic additional aerodynamic damping moment coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the pitch channel, ——engine thrust pitch channel moment, ——Pitch channel structure disturbance moment, Y kp ——Slosh normal displacement, Z kp ——Sway lateral displacement.

[0172] Step 7: Based on the elastic coordinate system E, the elastic vibration equation of the rocket body is established according to the rocket configuration parameters and engine related parameters.

[0173] Take the interface force analysis between a single nozzle and a single slosh bomb and the center body as an example Figure 5Based on the force and moment balance equations of the elastic coordinate system, considering the force and displacement coordination conditions, the full arrow vibration equation considering liquid sloshing and engine swing can be simplified as:

[0174]

[0175] The undamped free vibration equation of the rocket is:

[0176]

[0177] Where M is the central body mass matrix after considering liquid sloshing and engine swing, is the node elastic vibration displacement vector (including elastic deformation linear displacement and angular displacement); C and K are the damping matrix and stiffness matrix respectively, F is the node external load vector acting on the central body (including active forces such as gravity, thrust, aerodynamic force and interference force), is the inertial force and torque vector of the rigid body motion acting on the spring oscillator (equivalent mechanical model of liquid sloshing), are the rigid body inertia force and moment vector acting on the nozzle substructure, Represents the additional inertial force and torque vector for rigid body motion.

[0178] Solution The generalized eigenvalues ​​of the system can be obtained by obtaining the n non-zero natural frequencies ω1, ω2, …Ω, ω n and the corresponding mode column vector make is called the modal matrix.

[0179] According to the modal superposition method, the elastic vibration displacement It can be expanded using non-zero frequency modes, that is,

[0180]

[0181] Where q=[q1,q2,q3,…,q n ] T is the generalized coordinate matrix. Substitute the above formula into the full arrow vibration equation and use Φ T Multiply each term first, and we get

[0182]

[0183] make

[0184]

[0185] Then the elastic vibration equation can be simplified as

[0186]

[0187] In the formula, is the modal mass matrix, is the modal stiffness matrix, is the modal damping matrix, Q A is the generalized force of the external force (main force) acting on the central body, is the generalized force corresponding to the inertial force of rigid body motion acting on the central body, is the generalized force corresponding to the inertia force of rigid body motion acting on the spring oscillator, is the generalized force corresponding to the inertia force of rigid body motion acting on the nozzle.

[0188] If all tank liquid sloshing and all engine oscillating nozzles are considered, the whole arrow vibration equation can be written as

[0189]

[0190] That is, the full arrow vibration equation in generalized coordinate form is obtained, where the modal mass matrix The effects of liquid sloshing in all tanks and the swinging nozzles of all engines must be considered.

[0191] Since the mode is orthogonal to the mass matrix and the stiffness matrix, that is,

[0192]

[0193] In the formula, and are the generalized mass (also called modal mass) and generalized stiffness (also called modal stiffness) corresponding to the i-th order mode, so the generalized mass matrix and the generalized stiffness matrix are both diagonal matrices, that is,

[0194]

[0195] In the formula, ω i is the i-th order modal frequency, When the damping matrix C is Rayleigh damping, that is, C = αM + βK, the modal damping matrix is:

[0196]

[0197] In the formula, ξ i is the i-th order mode damping ratio.

[0198] By expanding the generalized force corresponding to the inertial force and the generalized force of the main force, the elastic vibration equation in expanded form can be obtained by deduction:

[0199]

[0200] In the formula, q i ——Generalized displacement of the i-th order elastic vibration, ω i ——i-th order elastic frequency, D4i ——Additional aerodynamic generalized force coefficient of rolling channel, ——the i-th order generalized aerodynamic damping force coefficient of the yaw channel, ——the i-th order generalized aerodynamic damping force coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the yaw channel, D 3iYks ——Generalized control force coefficient of engine Y-axis swing, D″ 3iYks ——Generalized inertia force coefficient of engine Y-axis swing, D 3iZks ——Generalized control force coefficient of engine Z-axis swing, D″ 3iZks ——Generalized inertia force coefficient of engine Z-axis swing, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing acceleration of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the yaw channel and the i-th order elastic coupling, ——Generalized force coefficient caused by the sloshing acceleration of the pth tank in the yaw channel and the i-th order elastic coupling, R ij ——generalized force coefficient of the j-th order elastic to i-th order elastic coupling, R′ ij ——the generalized damping force coefficient of the j-th order elastic to i-th order elastic coupling, The generalized force generated by the thrust of the sth engine in the kth module is, generalized forces arising from structural disturbances;

[0201] Step 8: Based on the rocket body coordinate system b, the propellant sloshing equation is established according to the rocket configuration parameters, engine related parameters, and the displacement and acceleration descriptions of the three types of structures on the rocket;

[0202] The propellant sloshing equation in the rocket body coordinate system b is written as:

[0203]

[0204] Compared with the existing liquid sloshing equation, are the sloshing displacement, velocity and acceleration of the sloshing mass in the rocket body coordinate system, where m kp is the propellant mass in case of thrust reduction failure; The coupling coefficient matrix of rocket body rotation and liquid sloshing under thrust reduction failure; is the coupling term caused by the displacement of the center of mass of the rocket body due to the thrust reduction failure, C kp is the equivalent damping coefficient matrix, is the equivalent stiffness coefficient matrix, is the apparent acceleration of the system, is the elastic motion-liquid sloshing coupling coefficient matrix, is the angular acceleration of the arrow body, is the elastic generalized acceleration.

[0205] Expanding the above formula and deriving it, the liquid sloshing equation can be obtained as follows:

[0206]

[0207] Where Y kp , Z kp ——respectively the normal and lateral displacements of the slosh, ζ kpy , kpz ——respectively normal and lateral sway damping ratios, Ω kp ——Sway frequency, Z Lkp ——Z position of the p tank of the kth booster in the rocket body coordinate system, Y Lkp ——Y position of the p tank of the kth booster in the rocket body coordinate system, E kpz ——the distance from the center of mass to the center of mass, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the yaw channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the yaw channel and the p-th tank sloshing, ——Rocket body pitch acceleration, ——Yaw angular acceleration of the rocket body, ——Rolling angular acceleration of the rocket body, n——Number of non-zero natural frequencies.

[0208] Step 9: Establish a sensitive element measurement equation based on the output of the rigid body rotation equation and the vibration mode and vibration mode slope at the elastic node position of the inertial device.

[0209] The inertial group is installed on the core stage rocket. The elastic deformation of the core stage will affect its measurement value, so the inertial group measurement equation is:

[0210]

[0211] In the formula, φ 0ix (x gz ) is the torsional vibration mode of the core stage rocket, φ 0iy (x gz ) and φ 0ix (x gz) is the vibration mode slope at the core stage rocket, both of which are in the elastic coordinate system o E -x E y E z E Given in; ——Inertial measurement value of pitch attitude angle considering elasticity, ——Rocket body pitch attitude angle, ψ gz ——Yaw attitude angle measured by the inertial group considering elasticity, ψ——Yaw attitude angle of the rocket body, γ gz ——Measured value of the rocket body rolling attitude angle considering elasticity, γ——rocket body rolling attitude angle.

[0212] Assuming that the rate gyro is installed on the kth booster, its measurement equation is:

[0213]

[0214] In the formula, φ kix (x st ) is the torsional vibration mode of the kth booster rate gyro installation, φ kiy (x st ) and φ kiz (x st ) is the vibration mode slope of the kth booster upper rate gyro installation, both of which are in the elastic coordinate system o E -x E y E z E Given in; ——Rocket body pitch angular velocity rate gyro measurement value, ——Rocket body pitch angular velocity, ——Measurement value of the rocket body yaw angular velocity rate gyro, ——Yaw angular velocity of the rocket body, ——Rocket body rolling angular velocity rate gyro measurement value, ——Rocket body rolling angular velocity, ——Elastic generalized velocity.

[0215] Step 10: Summarize the rigid body translation equation, rigid body rotation equation, rocket body elastic vibration equation, propellant sloshing equation, and sensitive element measurement equation to form a group of nonlinear dynamic equations.

[0216] Step 11: Perform a small perturbation linear expansion of the nonlinear dynamics equations based on the nominal reconstructed trajectory to obtain a linearized small deviation equation group suitable for controller design and stability analysis.

[0217] The linearized small deviation equations include:

[0218] 1) Small deviation equation of center of mass motion:

[0219]

[0220] ——dynamic coefficient of pitch channel, ——Yaw channel dynamic coefficient, ——Pitch channel gravity coefficient, ——Yaw channel gravity coefficient, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the yaw channel, Δ represents the deviation perturbation of the state quantity based on the nominal reconstructed trajectory.

[0221] 2) Small deviation equation of rigid body rotation:

[0222]

[0223]

[0224] 3) Linearization results of elastic motion equations:

[0225]

[0226] 4) The small deviation equation of shaking is:

[0227]

[0228] The present invention proposes a universal large-scale liquid carrier nonlinear dynamic modeling method, which is suitable for the dynamic description of different configurations of rockets under normal conditions and power system failures, ensuring the accuracy of the model under failures, and providing a basis for analyzing the impact of thrust reduction on rocket characteristics and dynamic behavior and formulating fault-tolerant and reconstruction control strategies for attitude control systems.

[0229] Although the present invention has been disclosed as above in the form of a preferred embodiment, it is not intended to limit the present invention. Any person skilled in the art may make possible changes and modifications to the technical solution of the present invention by using the methods and technical contents disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention shall fall within the protection scope of the technical solution of the present invention.

[0230] The contents not described in detail in the specification of the present invention belong to the common knowledge of those skilled in the art.

Claims

1. A dynamic modeling method for an attitude control system under power failure, characterized in that: The following steps are involved: Establish the coordinate system used for modeling, including the launch coordinate system a, the rocket body coordinate system b, the semi-speed coordinate system h, the elastic coordinate system E, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks ; Based on the launch vehicle object description, obtain the rocket configuration parameters and engine related parameters; Based on the launch coordinate system a, the rocket body coordinate system b, and the local elastic follower coordinate system e k(x) , nozzle coordinate system N ks , establish the displacement and acceleration description of three types of structures on the rocket, including the displacement and acceleration description of any point on the rocket body structure under the consideration of the elastic deformation of the rocket, the displacement and acceleration description of the swaying particle when the propellant sways, and the displacement and acceleration description of the center of mass of the nozzle of the kth booster sth engine; Based on the semi-velocity coordinate system h, the rigid body translation equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket; Based on the rocket body coordinate system b, the rigid body rotation equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket; Based on the elastic coordinate system E, the elastic vibration equation of the rocket body is established according to the rocket configuration parameters and engine related parameters; Based on the rocket body coordinate system b, the propellant sloshing equation is established according to the rocket configuration parameters and engine related parameters, and the displacement and acceleration description of the three types of structures on the rocket; According to the output of the rigid body rotation equation and the mode at the installation position of the inertial device, the measurement equation of the sensitive element is established; Summarize the rigid body translation equation, rigid body rotation equation, rocket body elastic vibration equation, propellant sloshing equation, and sensitive element measurement equation to form a nonlinear dynamic equation group; The nonlinear dynamic equations are linearized with small perturbations based on the nominal reconstructed trajectory to obtain a linearized small deviation equation set suitable for controller design and stability analysis.

2. The method for dynamic modeling of attitude control system under power failure according to claim 1 is characterized in that: Before conducting kinetic modeling, the modeling assumptions of the nonlinear kinetic model are as follows: Assume that the strain and rotation angle caused by the elastic deformation of the rocket body are 1% of the rocket body structure size; Assume that the liquid surface sway angle is less than 10°; Assume that the engine swing angle is less than 10°; Ignore the effect of the Earth's rotation; Only the reaction force generated by the ejection of propellant combustion products is considered and included in the thrust, while the additional Coriolis force, additional Coriolis moment and additional relative moment are ignored.

3. The method for dynamic modeling of attitude control system under power failure according to claim 1 is characterized in that: The definition of the arrow body coordinate system b is as follows: The origin of the rocket body coordinate system is the projection point O1 of the rocket instantaneous center of mass O2 on the core stage longitudinal axis. The O1X1 axis is the core stage longitudinal axis, pointing to the core stage head. The O1Y1 axis is located in the main symmetry plane of the rocket. The O1X1Y1 plane coincides with the OXY plane of the launch coordinate system at the moment of launch. The O1Z1 axis is determined according to the right-hand rule. The nozzle coordinate system N ks is defined as follows: Nozzle coordinate system N ks Connected to the nozzle of the kth booster sth engine, origin O δks Located at the engine swing point, O δks X δks Along the symmetry axis of the nozzle shell, it points to the head of the arrow body, that is, the thrust P ks For a single pendulum motor, the positive direction of the servo direction angle is defined as O δks Y δks Axis, O δks Z δks The axis is determined according to the right-hand rule. For the orthogonal double-pendulum engine, after the positive polarity of the servo direction swing angle is determined, the unique O δks Y δks Axis and O δks Z δks The positive direction of the axis.

4. The method for dynamic modeling of attitude control system under power failure according to claim 3 is characterized in that: The rocket configuration parameters and engine related parameters include: The rocket configuration parameters include: the number of boosters N zt , the total number of engines on the kth booster is The total number of tanks is When k = 0, it represents the core stage rocket body, that is, the total number of core stage rocket engines is The total number of tanks is Booster installation radius L k k represents the distance between the longitudinal axis of the booster and the longitudinal axis of the core stage, and the booster installation azimuth μ k represents the angle between the plane formed by the k-boost longitudinal axis and the core stage longitudinal axis and the O1X1Z1 plane; The engine related parameters include: each module engine installation position X kR represents the distance between the swing point of the kth booster engine and the theoretical vertex; the installation radius L of the sth engine in the module δks , represents the distance between the engine swing point and the longitudinal axis of the module; δ qks The engine installation azimuth angle represents the angle between the line between the engine swing point and the module longitudinal axis and the O1Y1 axis; is the engine swing azimuth, indicating the rotation from the rocket body coordinate system O1Y1 axis to the nozzle coordinate system O δks Y δks angle.

5. The method for dynamic modeling of attitude control system under power failure according to claim 4 is characterized in that: The description of displacement and acceleration of any point on the rocket body structure includes: The radius vector of any point on the rocket body structure: The acceleration of the structure particle is: In the formula, represents the absolute derivative in the emission coordinate system a, It represents the velocity of the origin O1 of the rocket coordinate system relative to the launch coordinate system. represents the local derivative in the rocket body coordinate system b, In the local elastic moving coordinate system e k(x) The local derivative in u k(x) is the origin o of the local elastic moving coordinate system after the projectile is deformed ek(x) The elastic translation displacement, R c.m is the position vector of the origin O1 of the rocket body coordinate system in the launch coordinate system a, R k(x) is the radius vector from O1 to the projection point of any point of the structure on the longitudinal axis of the k-th booster before the k-th booster is deformed, ω1 is the angular velocity of the rigid body rotation of the rocket relative to the launch coordinate system a, ω ek(x) is the local kinematic elastic coordinate system e k(x) The elastic rotation angular velocity relative to the rocket body coordinate system b, In the coordinate plane o ek(x) y ek(x) z ek(x) origin o ek(x) The position radius vector to any structural point, is the position radius vector of any structural point when the rocket is regarded as a rigid body; The description of the displacement and acceleration of the shaking particles when the propellant shakes includes: When the rocket deforms and the propellant sloshes, the sloshing mass m kp The position vector in the launch coordinate system a can be expressed as: The acceleration of the shaking particle is: In the formula, is the radius of the sway center of the rocket before deformation in the rocket body coordinate system b, The sway center is located in the local elastic follower coordinate system e k(x) The translational displacement, h kp is the swaying mass in the local elastic moving coordinate system e k(x) The sway displacement; The displacement and acceleration description of the center of mass of the nozzle of the k-th booster s-th engine includes: In the launch coordinate system a, the position radius vector of the kth booster sth engine nozzle center of mass is: The nozzle center of mass acceleration is: In the formula, Swing center O for the kth booster s engine δks The position vector of the projection point on the boost longitudinal axis in the rocket body coordinate system b is: To establish the elastic translation displacement of the origin of the local elastic servo coordinate system at this projection point, ρ rks is the local elastic moving coordinate system e k(x) The swing point radius vector, ρ Rks is the nozzle coordinate system N ks The radius vector of the nozzle center of mass, ρ r1ks is the radius vector of the swing point when the rocket is regarded as a rigid body, ρ R1ks is the radius vector of the nozzle mass center when the engine is not swaying, is the local elastic follower coordinate system Angular velocity of rotation relative to the rocket body coordinate system b, ω δks is the nozzle relative to the local elastic coordinate system e k(x) Angular velocity of rotation.

6. The method for dynamic modeling of attitude control system under power failure according to claim 5 is characterized in that: The step of establishing a rigid body translation equation comprises: According to Newton's second law, the rigid body translation equation of the entire rocket system is: In the formula, mg is gravity, P is engine thrust, F is aero is the aerodynamic force, F wind is the wind disturbance force, F structure is the structural disturbance force, is the volume fraction, For point Structure density at m kp is the sloshing mass of the p-th order sloshing of the k-th boost, m Rks is the swing mass of the k-th booster s-th engine, N zt is the number of boosters, is the total number of tanks, is the total number of engines on the kth booster; In the half-velocity coordinate system h, the coefficientized rigid body translation equation is expanded as follows: At the same time, the apparent acceleration in the rocket body coordinate system is established as: Among them, g is the acceleration due to gravity, V is the flight speed, θ is the ballistic inclination, σ is the ballistic deviation angle, α is the ballistic attack angle, α wp ——steady wind angle of attack, α wq ——Shear wind angle of attack, β——Ballistic sideslip angle, β wp ——Side slip angle in steady wind, β wq ——shear wind sideslip angle, ——Rocket body pitch angular velocity, ——Yaw angular velocity of the rocket body, q i ——elastic generalized displacement, ——elastic generalized displacement velocity, δ Yks , are the swing angle and swing angular acceleration of the swing engine in the Y direction, δ Zks , are the swing angle and swing angular acceleration of the swing engine in the Z direction, c X_aero ——axial force coefficient, ——pitch channel aerodynamic coefficient, ——Yaw channel aerodynamic coefficient, ——pitch channel aerodynamic damping force coefficient, ——Yaw channel aerodynamic damping force coefficient, ——the i-th order elastic additional aerodynamic coefficient of the pitch channel, ——is the i-th order elastic additional aerodynamic coefficient of the yaw channel, ——the aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, ——The aerodynamic coefficient of the i-th order elastic additional damping of the pitch channel, c 1P ——axial thrust coefficient, are the swing control force coefficients of the pitch channel engine in the Y and Z directions, are the swing control force coefficients of the yaw channel engine in the Y and Z directions, are the swing inertia force coefficients of the pitch channel engine in the Y and Z directions, are the swing inertia force coefficients of the yaw channel engine in the Y and Z directions, ——Pitch channel sway cross-link force coefficient, ——Yaw channel sway cross-linking force coefficient, ——the i-th order elastic additional thrust coefficient of the pitch channel, ——the i-th order elastic additional thrust coefficient of the yaw channel, ——Interference force coefficient of pitch channel structure, ——Yaw channel structure disturbance force coefficient, ——Coordinate transformation matrix from rocket body coordinate system to half-speed coordinate system, is the projection of the sway acceleration on the O1Y1 axis of the rocket body coordinate system, is the projection of the sway acceleration on the O1Z1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1X1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Y1 axis of the rocket body coordinate system, ——Apparent acceleration of the O1Z1 axis of the rocket body coordinate system.

7. The method for dynamic modeling of attitude control system under power failure according to claim 6 is characterized in that: The step of establishing a rigid body rotation equation comprises: The rigid body rotation equation of the entire rocket system is: Where M CM_g is the moment of gravity of each part on the origin O1 of the arrow body coordinate system, M CM_P is the torque of the engine thrust on the origin O1 of the rocket body coordinate system, M CM_aero is the moment of aerodynamic force on the origin O1 of the rocket coordinate system, M CM_wind is the moment of the wind disturbance force on the origin O1 of the arrow body coordinate system, M CM_structure is the moment of the structural disturbance force on the origin O1 of the rocket body coordinate system; After expanding and substituting the external torque term, the coefficientized full rigid body rotation equation in the arrow body coordinate system is obtained: In the formula, ——Rocket body rolling angular velocity, ——Roll channel yaw angle acceleration coupling coefficient, ——Roll channel pitch angle acceleration coupling coefficient, d PX ——Fault rolling channel mass center lateral displacement involvement coefficient, d1——Roll damping moment coefficient, d 3Yks ——k booster s engine Y-axis swing and roll control torque coefficient, d3″ Yks ——Y-axis swing inertia moment coefficient of rolling channel, d 3Zks ——Z-axis swing and rolling control force coefficient, d3″ Zks ——Z-axis swing inertia moment coefficient of rolling channel, d hykp —— Y-axis sloshing displacement of tank propellant to rolling channel cross-linking coefficient, d h ″ ykp ——Y-axis sway acceleration to rolling channel cross-linking coefficient, d hzkp ——Z-axis sway displacement to rolling channel cross-linking coefficient, d h ″ zkp ——Z-axis sway acceleration to rolling channel cross-linking coefficient, d 1i ——The i-th order elastic aerodynamic damping moment coefficient of the rolling channel, d 2i ——the i-th order elastic moment coefficient of the rolling channel, ——interference torque coefficient of rolling channel structure, ——Yaw channel rolling angular acceleration coupling coefficient, ——Yaw channel pitch angle acceleration coupling coefficient, ——The lateral displacement coefficient of the center of mass of the faulty yaw channel, ——Yaw channel aerodynamic damping moment coefficient, ——Yaw channel aerodynamic moment coefficient, ——Yaw control torque coefficient of engine Y-axis swing, ——Y-axis swing inertia moment coefficient of the engine, ——Engine Z-axis swing yaw control torque coefficient, ——Inertia moment coefficient of engine Z-axis swing, ——Yaw channel sway moment cross-coefficient, ——Yaw channel sway center of mass cross-connection coefficient, ——The i-th order elastic additional aerodynamic damping moment coefficient of the yaw channel, ——The i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the yaw channel, ——engine thrust yaw channel moment, ——Yaw channel structure disturbance torque, ——Roll angular acceleration coupling coefficient of pitch channel, ——pitch channel yaw angle acceleration coupling coefficient, ——the lateral displacement coefficient of the center of mass of the fault pitch channel, ——pitch channel aerodynamic damping moment coefficient, ——pitch channel aerodynamic moment coefficient, ——Y-axis pitch control torque coefficient of the engine, The Y-axis pitch inertia moment coefficient of the engine, Engine Z-axis pitch control moment coefficient, ——engine Z-axis pitch inertia moment coefficient, ——the coupling coefficient of the pitch channel sway moment, ——the pitch channel sway center of mass cross-connection coefficient, ——the i-th order elastic additional aerodynamic damping moment coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic moment and additional thrust moment coefficient of the pitch channel, ——engine thrust pitch channel moment, ——Pitch channel structure disturbance moment, Y kp ——Slosh normal displacement, Z kp ——Sway lateral displacement.

8. The method for dynamic modeling of attitude control system under power failure according to claim 7 is characterized in that: The step of establishing the elastic vibration equation of the rocket body comprises: In the formula, q i ——generalized displacement of the i-th order elastic vibration, ξ i ——damping ratio of the i-th order mode, ω i ——i-th order elastic frequency, D 4i ——Additional aerodynamic generalized force coefficient of rolling channel, ——the i-th order generalized aerodynamic damping force coefficient of the yaw channel, ——the i-th order generalized aerodynamic damping force coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the pitch channel, ——the i-th order generalized aerodynamic coefficient of the yaw channel, D 3iYks ——Generalized control force coefficient of engine Y-axis swing, D3″ iYks ——Generalized inertia force coefficient of engine Y-axis swing, D 3iZks ——Generalized control force coefficient of engine Z-axis swing, D3″ iZks ——Generalized inertia force coefficient of engine Z-axis swing, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing acceleration of the pth tank in the pitch channel and the i-th order elastic coupling, ——the generalized force coefficient caused by the sloshing displacement of the pth tank in the yaw channel and the i-th order elastic coupling, ——Generalized force coefficient caused by the sloshing acceleration of the pth tank in the yaw channel and the i-th order elastic coupling, R ij ——generalized force coefficient of the j-th order elastic to i-th order elastic coupling, R′ ij ——the generalized damping force coefficient of the j-th order elastic to i-th order elastic coupling, The generalized force generated by the thrust of the sth engine in the kth module is, Generalized force due to structural disturbance, Y kp ——Slosh normal displacement, Z kp —— lateral displacement of sway, ——sway normal acceleration, ——Sway lateral acceleration.

9. The method for dynamic modeling of attitude control system under power failure according to claim 8 is characterized in that: The establishing of the propellant sloshing equation comprises: The propellant sloshing equation in the rocket body coordinate system b is: are the sloshing displacement, velocity, and acceleration of the sloshing mass in the rocket body coordinate system, respectively. m kp is the propellant mass in case of thrust reduction failure, is the coupling coefficient matrix of rocket body rotation and liquid sloshing under thrust reduction failure, is the coupling term caused by the displacement of the center of mass of the rocket body due to the thrust reduction failure, is the equivalent damping coefficient matrix, is the equivalent stiffness coefficient matrix, is the apparent acceleration of the system, is the elastic motion-liquid sloshing coupling coefficient matrix, is the angular acceleration of the arrow body, is the elastic generalized acceleration; Expand the above formula and derive the liquid sloshing equation as follows: ζ kpy , kpz are the normal and lateral sway damping ratios, Ω kp ——Sway frequency, Z Lkp ——Z position of the p tank of the kth booster in the rocket body coordinate system, Y Lkp ——Y position of the p tank of the kth booster in the rocket body coordinate system, E kpz ——the distance from the center of mass to the center of mass, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation displacement of the yaw channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the pitch channel and the p-th tank sloshing, ——Force coefficient of coupling between the i-th order elastic deformation acceleration of the yaw channel and the p-th tank sloshing, ——Rocket body pitch acceleration, ——Yaw angular acceleration of the rocket body, ——Rolling angular acceleration of the rocket body, n——Number of non-zero natural frequencies.

10. The method for dynamic modeling of attitude control system under power failure according to claim 9, characterized in that: The establishment of the sensitive element measurement equation includes: The inertial measurement equation is: In the formula, φ 0ix (x gz ) is the torsional vibration mode of the core stage rocket, φ 0iy (x gz ) and φ 0ix (x gz ) is the vibration mode slope at the core stage rocket, ——Inertial measurement of pitch attitude angle considering elasticity, ——Rocket body pitch attitude angle, ψ gz ——Yaw attitude angle measured by the inertial group considering elasticity, ψ——Yaw attitude angle of the rocket body, γ gz ——Measured value of the rocket body rolling attitude angle considering elasticity, γ——rocket body rolling attitude angle; Assuming that the rate gyro is installed on the kth booster, its measurement equation is: In the formula, φ kix (x st ) is the torsional vibration mode of the kth booster rate gyro installation, φ kiy (x st ) and φ kiz (x st ) is the vibration mode slope at the installation location of the kth booster upper rate gyro, ——Rocket body pitch angular velocity rate gyro measurement value, ——Rocket body pitch angular velocity, ——Measurement value of the rocket body yaw angular velocity rate gyro, ——Yaw angular velocity of the rocket body, ——Rocket body rolling angular velocity rate gyro measurement value, ——Rocket body rolling angular velocity, ——Elastic generalized velocity.

11. The method for dynamic modeling of attitude control system under power failure according to claim 10, characterized in that: The linearized small deviation equation group includes: 1) Small deviation equation of center of mass motion: ——dynamic coefficient of pitch channel, ——Yaw channel dynamic coefficient, ——gravity coefficient of pitch channel, ——Yaw channel gravity coefficient, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the pitch channel, ——the i-th order elastic additional aerodynamic force and thrust coefficient of the yaw channel, Δ represents the deviation perturbation of the state quantity based on the nominal reconstructed trajectory; 2) Small deviation equation of rigid body rotation: 3) Linearization results of elastic motion equations: 4) The small deviation equation of shaking is:

Images