Design method of rarefied atmosphere pneumatic assisted rendezvous maneuver strategy

Abstract

The invention discloses a rarefied atmosphere pneumatic auxiliary rendezvous maneuver strategy design method, and belongs to the technical field of aerospace, and the method comprises the steps: building a polar coordinate dynamic model of a Mars probe flying in the Mars atmosphere and an external flight dynamic model; calculating an intersection position of different planes or an intersection phase of a given same orbit plane; roughly matching rendezvous phases; accurate matching of intersection phases; a pneumatic auxiliary rendezvous maneuvering track is obtained; according to the rarefied atmosphere pneumatic auxiliary rendezvous maneuvering strategy design method, pneumatic auxiliary rendezvous task planning of any number of turns can be achieved, the calculation speed is high, and the method can be applied to the task requirement of online trajectory planning; the trajectory result can be combined with an existing nonlinear programming algorithm, and the optimality of the trajectory is further improved.

Description

Technical Field

[0001] The present invention relates to the field of aerospace technology, and particularly to a method for designing a thin - atmosphere aerodynamic - assisted rendezvous maneuver strategy. Background Art

[0002] In the field of space exploration, the development of spacecraft orbit adjustment technology is crucial. As an advanced technology for adjusting the spacecraft orbit by using the planetary atmosphere, aerodynamic - assisted maneuver can significantly save fuel by utilizing aerodynamic force to attenuate excess energy compared with traditional space maneuvers. This is of great significance for improving the maneuverability of spacecraft and promoting scientific space exploration.

[0003] Aerodynamic - assisted orbit maneuvers cover a series of special technologies. Among them, aerodynamic - assisted orbit rendezvous and interception technologies have attracted much attention in recent research due to their significant advantages. However, due to the strong non - linear characteristics of aerodynamic problems and the extremely complex endpoint constraints of aerodynamic - assisted problems, most of the existing technologies and methods are only applicable to single - loop aerodynamic - assisted orbit transfer. However, single - loop aerodynamic orbit transfer has extremely high requirements for spacecraft thermal protection. In actual engineering tasks, single - loop aerodynamic - assisted orbit transfer tasks are often decomposed into multi - loop aerodynamic orbit transfers. By extending the aerodynamic - assisted time, the negative impact brought by thermal ablation can be effectively alleviated, thereby greatly improving the engineering feasibility and cost - effectiveness of aerodynamic orbit transfer.

[0004] However, the trajectory planning problem of multi - loop aerodynamic - assisted de - orbit involves many complex constraints. After discretizing this problem, the problem scale increases sharply, which brings great difficulties to directly solving the trajectory planning problem and leads to a significant extension of the solution time. In this context, there is an urgent need for an innovative method to break through this dilemma and achieve the rapid generation of aerodynamic - assisted rendezvous trajectories. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for designing a thin - atmosphere aerodynamic - assisted rendezvous maneuver strategy to solve the problems existing in the background art.

[0006] To achieve the above purpose, the present invention provides a method for designing a thin - atmosphere aerodynamic - assisted rendezvous maneuver strategy, including the following steps:

[0007] S1. Establish a polar - coordinate dynamics model of a Mars probe flying in the Mars atmosphere and an extra - atmospheric flight dynamics model;

[0008] S2. Calculate the out - of - plane rendezvous position or the given rendezvous phase in the same orbital plane;

[0009] S3. Rough matching of the rendezvous phase;

[0010] S4. Precise matching of the rendezvous phase;

[0011] S5. Obtain the aerodynamic - assisted rendezvous maneuver trajectory.

[0012] Preferably, the polar coordinate dynamics of the Mars probe flying in the Martian atmosphere in S1 are as follows:

[0013]

[0014]

[0015] where g r , g φ represent the radial and tangential gravitational accelerations respectively; their expressions are:

[0016]

[0017] where r is the distance between the center of mass of the aircraft and the center of the planet; V is the magnitude of the velocity; γ is the flight path angle or track angle of the aircraft, ψ is the heading angle of the aircraft; σ is the bank angle of the aircraft, which is the control variable; θ and φ are the longitude and latitude respectively; ω m is the angular velocity of the planet's rotation; μ is the gravitational constant of the planet, r M is the radius of the planet; J2 is the second-order spherical harmonic coefficient of the planet; the lift acceleration L and the drag acceleration D are:

[0018]

[0019] where S is the reference area of the aircraft; C L and C D are the lift coefficient and the drag coefficient respectively; the atmospheric density ρ = ρ0exp(-h / h s ), where ρ0 is the atmospheric density on the planet's surface; h s is the atmospheric density coefficient; m is the mass of the aircraft, h is the current altitude of the aircraft, h = r - rM.

[0020] Preferably, the specific content of S2 is as follows:

[0021] For the case of rendezvous in the same orbital plane (the difference in orbital inclination is less than 1 deg, and the difference in right ascension of the ascending node is less than 0.5 deg), it can be directly specified by the user according to the mission requirements. If not directly specified, the default of this method is to rendezvous at a phase of 45 deg. For the case of non-coplanar orbital rendezvous, the following calculation and processing are carried out.

[0022] Let the six orbital elements of the initial orbit of the aircraft be [a1, e1, ω1, Ω1, i1, θ1], which represent the semi-major axis, eccentricity, argument of perigee, right ascension of the ascending node, orbital inclination, and true anomaly respectively. Then the normal vector n1 of the orbital plane can be calculated according to the definition of the orbital elements:

[0023] n1 = [sini1sinΩ1, -sini1cosΩ1, cosi1](10)

[0024] Similarly, according to the target orbital elements [a2, e2, ω2, Ω2, i2, θ2], the normal vector n2 of the target orbital plane is calculated as follows:

[0025] n2 = [sin i2 sin Ω2, -sin i2 cos Ω2, cos i2] (11)

[0026] The cross product of the two normal vectors of the orbital planes gives the vector n3 of the intersection line of the two orbital planes:

[0027] n3 = n1 × n2 (12)

[0028] Then the rendezvous phase is located at the intersection point of this vector on the target orbit, and this intersection point is projected onto the planetary Cartesian inertial system. When the X-axis of the Cartesian inertial system is rotated in the right-hand direction around the positive Z-axis to coincide with this projection, the swept angle is the rendezvous phase. The rendezvous phase is defined in [0, 2π]. For any two skew orbits, there should be two calculations of the rendezvous phase, and any one of them can be taken as the rendezvous phase, or it can be selected according to the needs of the engineering mission. During specific calculations, assume

[0029] n3 = [x3, y3, z3] (13)

[0030] Assume the rendezvous phases obtained according to this definition are ψ1 and ψ2. When x3 ≥ 0 and y3 > 0,

[0031] ψ1 = arctan(y3 / x3) (14)

[0032] When x3 ≥ 0 and y3 ≤ 0,

[0033] ψ1 = 2π + arctan(y3 / x3) (15)

[0034] When x3 ≤ 0,

[0035] ψ1 = π + arctan(y3 / x3) (16)

[0036] When calculating ψ2, just consider n3 as the opposite vector and calculate again.

[0037] Preferably, the content of S3 is as follows:

[0038] First, the first-step integration obtains the deorbiting state, that is, integrating from the initial position of the aircraft to the current phase of the aircraft being equal to the target phase. Assume the initial state of the aircraft is x0 = [r0, v0], and assume the cut-off condition for the orbital integration is that the current phase is the target orbital phase. Assume the phase of the aircraft is φ a , and the initial time is t0. g(x) is the orbital dynamics equation. Then the integration to the deorbiting phase can be expressed as:

[0039]

[0040] In the above formula, x a = [r a , v a , t a is the time at this moment. The given off-orbit altitude h can be calculated through the off-orbit phase a0 = ||r a ||. Let the apoapsis altitude of the target orbit be h afaim . The relationship between the orbital elements and the periapsis velocity is:

[0041]

[0042] where the semi-major axis of the orbit a = h a + h p + 2R, the orbital eccentricity e = (h a - h p ) / (2a); μ is the planetary gravitational constant, h a is the apoapsis altitude, and h a0 is substituted.

[0043] When h p = h atm , the maximum value of the periapsis velocity is calculated according to the initial apoapsis altitude, and the minimum value of the periapsis velocity is obtained according to the apoapsis altitude of the target orbit. The minimum and maximum values of the velocity at the periapsis are V p,low , V p,up respectively, and the calculation formula is:

[0044]

[0045] where R represents the radius of the planet, and r M is substituted; h atm represents the atmospheric altitude.

[0046] A periapsis sequence [V p1 , V p2 ,..., V pN is constructed arithmetically between the minimum and maximum values of the periapsis altitude, which is expressed as:

[0047]

[0048] where ΔV p = (V p,up - V p,low ) / (N - 1) is the step size of the periapsis velocity sequence, and N is the given total number of orbits; according to the relationship formula (17) between the orbital elements and the periapsis velocity, the altitude sequence [h a1 , h a2 ,..., haN :

[0049]

[0050] where V pn = V p,low + (n - 1)ΔV p ; where R represents the planet radius, use r M and substitute it. Using this altitude sequence to solve the cooperative pulse quantity by the sequence bisection method, a full - course aerodynamic - assisted maneuver trajectory can be constructed.

[0051] Let the pulse applied at the apocenter after the n - th atmospheric crossing be ΔV n , and the initial - end pulse be ΔV0; r a0 , v a0 are respectively the initial - orbit position and velocity vector, given as the orbit state at the apocenter of an elliptical orbit; are respectively the apocenter position and velocity vector reached after flying out of the atmosphere after the (n - 1) - th atmospheric crossing; then the position - velocity vector after applying the pulse is the initial - apocenter position vector of the n - th circle, denoted as ΔV n and the pulse - application directions of ΔV0 are the same as direction, then the position - velocity vector satisfies the following formula:

[0052]

[0053] Let the cut - off condition of the orbit integration be Γ0 = ||r|| - r atm = 0. From the initial time t ini , the state vector after applying the pulse is integrated to the atmosphere inlet and expressed as:

[0054]

[0055] where are respectively the position vector, velocity vector, and time of the atmosphere inlet; g(r, v) is the general orbit - dynamics differential equation considering the two - body perturbation.

[0056] The transformation from the rv vector in the Cartesian inertial coordinate system to the atmosphere - flight state is set as the following function:

[0057]

[0058] Further, from the atmosphere - inlet state the outlet state is obtained by integrating through the atmosphere - flight dynamics equation, and the integration cut - off condition is Γ0 = r - r atm= 0, the control variable bank angle σ in dynamics is taken as a constant value; expressed as:

[0059]

[0060] where f(x ) is the atmospheric flight dynamics equation as shown in (1)-(6).

[0061] Further, through the function Φ x→rv the atmospheric exit position velocity vector in the Cartesian inertial coordinate system is obtained

[0062]

[0063] Finally, the apocenter state is obtained through orbit integration. The integration cut-off condition is Γ1 = r·v = 0, and the initial conditions are

[0064]

[0065] are respectively the position vector, velocity vector and time before applying the pulse at the apogee; let the apocenter height that can be reached when flying out of the atmosphere in the nth circle be satisfies the equation:

[0066]

[0067] Construct the pulse ΔV applied at the apocenter after the nth atmospheric crossing n or the initial pulse ΔV0, the apocenter position and velocity vectors after flying out of the atmosphere after the (n - 1)th atmospheric crossing and the apocenter position velocity vector when flying out of the atmosphere in the nth circle apocenter height The functional relationship is expressed by the equation as:

[0068]

[0069] The construction of the cooperative pulse and the non-linear mapping relationship of the exit apocenter is completed. For the first pulse ΔV0, the mapping is constructed similarly as follows:

[0070]

[0071] where r a0 , v a0 are the position and velocity vectors of the initial orbit position.

[0072] Let the output of equation (30) be When, its mapping is expressed as:

[0073]

[0074] The obtained altitude sequence, where the nth element is the target apoapsis altitude of the nth atmospheric crossing. Let the output of (31) be this target apoapsis altitude. When calculating for the nth orbit, it has been obtained from the calculation of the previous orbit, and a univariate non - linear equation about ΔV n is constructed:

[0075]

[0076] where h an is the constant value calculated in (22); then it is solved by the bisection method, and the initial solution boundary is taken as [- 10,10] m / s.

[0077] The obtained position vector of the apoapsis at the end of the aerodynamic - assisted maneuver is The corresponding rendezvous phase is φ b . Then the phase deviation can be expressed as Δφ = φ b - ψ1. Therefore, to reduce this phase deviation, the de - orbit phase fine - tuning should be used to adapt to this phase change. So let the de - orbit phase to be corrected be φ′ a . Its calculation formula is:

[0078] φ′ a = φ a - Δφ (34)

[0079] After calculating the new de - orbit phase, the continuously corrected de - orbit phase can be obtained by repeating the above steps according to the above steps. Set the stop condition of the iteration as the following formula.

[0080] |φ′ a - φ a |≤δ φ (35)

[0081] where δ φ is the iteration error limit; after the iteration satisfies the iteration condition, the rough phase matching is completed. The obtained trajectory satisfies the terminal altitude constraint and basically satisfies the position constraint. However, since the above trajectory generation does not consider the mission time, if the aircraft flies according to the above trajectory, when the aircraft reaches the interaction position, the target satellite may not have reached the target phase. Therefore, it is necessary to adjust the calculated aerodynamic - assisted maneuver trajectory so that the arrival time of the aircraft is exactly equal to the arrival time of the target satellite without affecting the end - point arrival position.

[0082] Preferably, the specific content of S4 is:

[0083] Assume that the number of loops of the generated aerodynamic - assisted maneuver is N, and the orbital altitude during the de - orbit maneuver is h a0, after each subsequent circle of atmospheric flight, the apocenter heights reached are successively [h a1 , h a2 ,..., h aN , and the perigee heights are successively [h p1 , h p2 ,..., h pN . Then the following calculates the feasible rendezvous time series of the target satellite. Assume the initial orbital state of the target satellite is x t0 = [r t0 , v t0 , and let the cut-off condition of the orbital integration be that the current phase is the target orbital phase. Let the phase of the target vehicle be φ ta , and the initial time be t0. Then the integration to the target rendezvous can be expressed as:

[0084]

[0085] According to the initial orbital elements [a2, e2, ω2, Ω2, i2, θ2] of the target satellite, the orbital period of the target satellite can be calculated as:

[0086]

[0087] Therefore, the time series when the target satellite can reach the target phase is:

[0088] t tar = t ta + kT t , k = 1, 2, 3,.... (38)

[0089] Since the influence of the apocenter height of the vehicle on the terminal rendezvous phase is relatively small, by slightly adjusting the apocenter height at the exit of the first atmospheric flight of the vehicle, the time adjustment of the aerodynamic assisted rendezvous can be achieved.

[0090] According to the apocenter and perigee sequences of the aerodynamic assisted maneuver before adjustment, the orbital period of the first circle can be calculated. It should be noted that here the first circle refers to starting from when the vehicle first reaches the perigee until the vehicle reaches the perigee for the second time. Since during a single atmospheric crossing process, the change in perigee height generally does not exceed 50 km, while the orbital apocenter height is usually greater than 500 km. Therefore, when calculating the flight time of the first circle, the perigee height can be directly regarded as the perigee of the first atmospheric flight, without separately calculating the two half-periods. The orbital period of the first circle before adjustment is:

[0091]

[0092] h p1 is the perigee height; let h a1 The adjusted apocenter height is h'a1 , and its feasible range is:

[0093] h a2 < h′ a1 < h a0 (40)

[0094] The orbital period of the first lap after adjustment is:

[0095]

[0096] Let the time when the aircraft reaches the apocenter position after flying out of the atmosphere for the last time be That is, the time when the aircraft is expected to reach the target phase. Then the conditions that the rendezvous needs to meet are:

[0097]

[0098] where t tar is the time for the target to reach the target phase; the time for the adjusted aircraft to reach the rendezvous phase is It can be calculated by the following formula:

[0099]

[0100] Then the conditions that the adjusted rendezvous needs to meet are:

[0101]

[0102] Substitute the t tar expression, and the apocenter orbital altitude of the first lap after adjustment can be derived as

[0103]

[0104] Therefore, a suitable h can be selected according to the feasible range of the orbital altitude a1 , thereby determining the entire aerodynamic assist maneuver trajectory. After correcting to obtain the new orbital altitude of the first lap, the adjusted apocenter altitude sequence can be obtained as [h′ a1 ,h a2 ,...,h aN . By solving equation (33) lap by lap, the entire aerodynamic assist maneuver trajectory can be obtained.

[0105] Therefore, the present invention adopts the above-mentioned method for designing a thin atmosphere aerodynamic assist rendezvous maneuver strategy, and has the following beneficial effects:

[0106] (1) This method can not only handle the aerodynamic mission planning of any number of laps, but also has the characteristics of fast calculation speed, and is very suitable for real-time online trajectory planning requirements;

[0107] (2) The generated initial trajectory provides a good basis for subsequent optimization, and by combining advanced non-linear programming techniques, better trajectory planning results can be achieved;

[0108] (3) The algorithm proposed in the present invention uses an analytically corrected feasible trajectory and has strong robustness.

[0109] The technical solution of the present invention will be further described in detail below with reference to the drawings and embodiments. Description of the Drawings

[0110] Figure 1 It is a schematic flow chart of a design method for a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy of the present invention;

[0111] Figure 2 It is a schematic diagram of the full-course trajectory obtained by a design method for a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy of the present invention. Detailed Embodiments

[0112] The following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0113] Please refer to Figure 1 , a design method for a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy,

[0114] Therefore, the present invention adopts the above-mentioned design method for a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy, including the following steps:

[0115] S1. Establish a polar coordinate dynamics model for the Mars probe flying in the Mars atmosphere and an extra-atmospheric flight dynamics model.

[0116] Preferably, the polar coordinate dynamics of the Mars probe flying in the Mars atmosphere in S1 is:

[0117]

[0118]

[0119] Among them, g r , g φ respectively represent the radial and tangential gravitational accelerations; their expressions are:

[0120]

[0121] Wherein, r is the distance between the center of mass of the aircraft and the center of the planet; V is the magnitude of the velocity; γ is the flight path angle or track angle of the aircraft, ψ is the heading angle of the aircraft; σ is the bank angle of the aircraft, which is the control variable; θ and φ are the longitude and latitude respectively; ω m is the angular velocity of the planet's rotation; μ is the gravitational constant of the planet, r M is the radius of the planet; J2 is the second-order spherical harmonic coefficient of the planet; the lift acceleration L and the drag acceleration D are:

[0122]

[0123] Wherein, S is the reference area of the aircraft; C L and C D are the lift coefficient and the drag coefficient respectively; the atmospheric density ρ = ρ0exp(-h / h s ), where ρ0 is the atmospheric density on the planet's surface; h s is the atmospheric density coefficient; m is the mass of the aircraft, h is the current altitude of the aircraft, h = r - r M .

[0124] S2. Calculate the off-plane rendezvous position or specify the rendezvous phase on the same orbital plane.

[0125] The specific content of S2 is as follows:

[0126] For the case of rendezvous on the same orbital plane (the difference in orbital inclination is less than 1 deg, and the difference in right ascension of the ascending node is less than 0.5 deg), it can be directly specified by the user according to the mission requirements. If not directly specified, the default of this method is to rendezvous at a phase of 45 deg. For off-plane orbital rendezvous, the following calculations are performed.

[0127] Let the six orbital elements of the initial orbit of the aircraft be [a1, e1, ω1, Ω1, i1, θ1], which represent the semi-major axis, eccentricity, argument of perigee, right ascension of the ascending node, orbital inclination, and true anomaly respectively. Then the normal vector n1 of the orbital plane can be calculated according to the definition of the orbital elements:

[0128] n1 = [sini1sinΩ1, -sini1cosΩ1, cosi1] (10)

[0129] Similarly, according to the target orbital elements [a2, e2, ω2, Ω2, i2, θ2], the normal vector n2 of the target orbital plane is calculated:

[0130] n2 = [sini2sinΩ2, -sini2cosΩ2, cosi2] (11)

[0131] The cross product of the normal vectors of the two orbital planes is used to obtain the vector n3 of the intersection line of the two orbital planes:

[0132] n3 = n1 × n2 (12)

[0133] Then the intersectable phase is located at the intersection of this vector on the target orbit, and project this intersection onto the planetary Cartesian inertial system. When the X-axis of the Cartesian inertial system is rotated in the right-hand direction around the positive direction of the Z-axis to coincide with this projection, the swept angle is the intersection phase. The intersection phase is defined in [0, 2π]. For any two skew orbits, there should be two calculations for their intersection phase, and any one of them can be taken as the intersection phase, or it can be selected according to the needs of the engineering mission. During specific calculation, let

[0134] n3 = [x3, y3, z3] (13)

[0135] Let the intersection phases obtained according to this definition be ψ1 and ψ2. When x3 ≥ 0 and y3 > 0,

[0136] ψ1 = arctan(y3 / x3) (14)

[0137] When x3 ≥ 0 and y3 ≤ 0,

[0138] ψ1 = 2π + arctan(y3 / x3) (15)

[0139] When x3 ≤ 0,

[0140] ψ1 = π + arctan(y3 / x3)(16)

[0141] When calculating ψ2, just consider n3 as the opposite vector and calculate again.

[0142] S3. Coarse matching of the intersection phase.

[0143] The content of S3 is as follows:

[0144] First, integrate in the first step to obtain the deorbiting state, that is, integrate from the initial position of the aircraft to the current phase of the aircraft equal to the target phase. Let the initial state of the aircraft be x0 = [r0, v0], and let the cut-off condition of the orbit integration be that the current phase is the target orbit phase. Let the phase of the aircraft be φ a , and the initial time be t0. g(x) is the orbit dynamics equation. Then the integration to the deorbiting phase can be expressed as:

[0145]

[0146] In the above formula, x a = [r a , v a , t a is the time at this moment. The given deorbiting orbit altitude h a0 = ||r a||Let the apoapsis altitude of the target orbit be h afaim . The relationship between the orbital elements and the periapsis velocity is as follows:

[0147]

[0148] where the semi-major axis of the orbit a = h a +h p +2R, and the orbital eccentricity e = (h a -h p ) / (2a); μ is the planetary gravitational constant, h a is the apoapsis altitude, and use h a0 to substitute.

[0149] When h p =h atm , calculate the maximum periapsis velocity according to the initial apoapsis altitude, and obtain the minimum periapsis velocity according to the apoapsis altitude of the target orbit. The minimum and maximum velocities at the periapsis are V p,low ,V p,up , and the calculation formula is:

[0150]

[0151]

[0152] where R represents the radius of the planet, and use r M to substitute; h atm represents the atmospheric altitude.

[0153] Construct an arithmetic sequence of periapsis points [V p1 ,V p2 ,...,V pN between the minimum and maximum periapsis altitudes, which is expressed as:

[0154]

[0155] where ΔV p =(V p,up -V p,low ) / (N - 1) is the step size of the periapsis velocity sequence, and N is the given total number of orbits; according to the relationship formula (17) between the orbital elements and the periapsis velocity, solve the altitude sequence [h a1 ,h a2 ,...,h aN :

[0156]

[0157] where V pn =V p,low +(n - 1)ΔV p; R represents the radius of the planet, using r M Substitute it; Using this altitude sequence to solve the cooperative pulse quantity by sequence dichotomy, a full-course aerodynamic-assisted maneuver trajectory can be constructed.

[0158] Let the pulse applied at the apocenter after the nth atmospheric crossing be ΔV n , and the initial pulse be ΔV0; r a0 , v a0 are the initial orbital position and velocity vector respectively, given as the orbital state at the apocenter of an elliptical orbit; are the apocenter position and velocity vector reached after flying out of the atmosphere after the (n - 1)th atmospheric crossing respectively; then the position and velocity vector after applying the pulse is the initial apocenter position vector of the nth orbit, denoted as ; ΔV n and the pulse application directions of ΔV0 are the same as the direction, then the position and velocity vector satisfies the following equation:

[0159]

[0160] Let the cut-off condition of the orbit integral be Γ0 = ||r|| - r atm = 0. From the initial time t ini , the state vector after applying the pulse is integrated to the atmosphere entry and expressed as:

[0161]

[0162] where are the position vector, velocity vector, and time of the atmosphere entry respectively; g(r, v) is the general orbit dynamics differential equation considering the two-body perturbation.

[0163] Converting from the rv vector in the Cartesian inertial coordinate system to the atmospheric flight state is set as the following function:

[0164]

[0165] Furthermore, from the atmosphere entry state the exit state is obtained by integrating through the atmospheric flight dynamics equation The integration cut-off condition is Γ0 = r - r atm = 0, and the control variable bank angle σ in the dynamics is taken as a constant value; it is expressed as:

[0166]

[0167] where f(x) is the atmospheric flight dynamics equation as shown in (1)-(6). Further through the function Φ x→rvObtain the velocity vector of the atmospheric exit position in the Cartesian inertial coordinate system

[0168]

[0169] Finally, the apocenter state is obtained through orbital integration. The integration cut-off condition is Γ1 = r·v = 0, and the initial conditions are

[0170]

[0171] They are respectively the position vector, velocity vector, and time before applying the pulse at the apogee; Let the apocenter altitude that can be reached when flying out of the atmosphere in the nth orbit be Satisfy the equation:

[0172]

[0173] Construct the pulse ΔV applied at the apocenter after the nth atmospheric crossing n Or the initial pulse ΔV0, the apocenter position and velocity vectors after flying out of the atmosphere after the (n - 1)th atmospheric crossing And the apocenter position and velocity vectors when flying out of the atmosphere in the nth orbit Apocenter altitude The functional relationship, expressed by the equation as:

[0174]

[0175] The construction of the non-linear mapping relationship between the cooperative pulse and the exit apocenter is completed. For the first pulse ΔV0, the mapping is constructed as follows in the same way:

[0176]

[0177] Where r a0 , v a0 Are the position and velocity vectors of the initial orbit position.

[0178] Let the output of equation (30) be When, its mapping is expressed as:

[0179]

[0180] Through the obtained altitude sequence, the nth element is the target apocenter altitude of the nth atmospheric crossing. Let the output of (31) be this target apocenter altitude. When calculating in the nth orbit, Has been obtained from the calculation of the previous orbit. Construct a single-variable non-linear equation about ΔV n As:

[0181]

[0182] Here, h an is the constant value calculated in (22); then, the bisection method is used to solve, and the initial solution boundary is taken as [-10, 10] m / s.

[0183] The position vector of the apocenter at the end of the aerodynamic-assisted maneuver obtained by solving is The corresponding rendezvous phase is φ b . Then the phase deviation can be expressed as Δφ = φ b - ψ1. Therefore, to reduce this phase deviation, the off-orbit phase fine-tuning should be used to adapt to this phase change. Thus, it is assumed that the off-orbit phase should be corrected to φ' a . Its calculation formula is:

[0184] φ' a = φ a - Δφ (34)

[0185] After calculating the new off-orbit phase, the continuously corrected off-orbit phase can be obtained by repeating the above steps. The stopping condition for iteration is set as the following formula.

[0186] |φ' a - φ a | ≤ δ φ (35)

[0187] where δ φ is the iteration error limit. After the iteration satisfies the iteration condition, the rough phase matching is completed. The obtained trajectory satisfies the terminal altitude constraint and basically satisfies the position constraint. However, since the above trajectory generation does not consider the mission time, if the aircraft flies according to the above trajectory, when the aircraft reaches the interaction position, the target satellite may not have reached the target phase yet. Therefore, it is necessary to adjust the calculated aerodynamic-assisted maneuver trajectory so that the arrival time of the aircraft is exactly equal to the arrival time of the target satellite without affecting the end arrival position.

[0188] S4. Precise matching of the rendezvous phase.

[0189] The specific content of S4 is:

[0190] Assume that the number of loops of the generated aerodynamic-assisted maneuver is N, and the orbital altitude during the off-orbit maneuver is h a0 , and after each loop of atmospheric flight, the apocenter altitudes reached are [h a1 , h a2 ,..., h aN , and the perigee altitudes are [h p1 , h p2 ,..., h pN , then the following calculates the feasible rendezvous time series of the target satellite. Assume that the initial orbital state of the target satellite is xt0 = [r t0 , v t0 , set the cut-off condition of the orbital integral as the current phase being the target orbital phase. Let the phase of the target vehicle be φ ta , and the initial time be t0. Then the integration to the target rendezvous can be expressed as:

[0191]

[0192] According to the initial orbital elements [a2, e2, ω2, Ω2, i2, θ2] of the target satellite, the orbital period of the target satellite can be calculated as:

[0193]

[0194] Therefore, the time series when the target satellite can reach the target phase is:

[0195] t tar = t ta + kT t , k = 1, 2, 3,.... (38)

[0196] Since the apocenter altitude of the vehicle has little influence on the terminal rendezvous phase, by slightly adjusting the apocenter altitude at the exit of the first atmospheric flight of the vehicle, the time adjustment of the aerodynamic-assisted rendezvous can be achieved.

[0197] According to the apocenter and perigee sequences of the aerodynamic-assisted maneuver before adjustment, the orbital period of the first circle can be calculated. It should be noted that here the first circle refers to starting from the moment when the vehicle first reaches the perigee until the vehicle reaches the perigee for the second time. Due to the single atmospheric crossing process, the perigee altitude generally does not change by more than 50 km, while the orbital apocenter altitude is usually greater than 500 km. Therefore, when calculating the flight time of the first circle, the perigee altitude can be directly regarded as the perigee of the first atmospheric flight, without calculating separately for two half-periods. The orbital period of the first circle before adjustment is:

[0198]

[0199] h p1 is the perigee altitude; let h a1 The adjusted apocenter altitude is h' a1 , and its feasible range is:

[0200] h a2 <h' a1 <h a0 (40)

[0201] The orbital period of the first circle after adjustment is:

[0202]

[0203] Let the time when the vehicle reaches the apocenter position after flying out of the atmosphere for the last time be That is, the time when the vehicle is expected to reach the target phase. Then the conditions that the rendezvous needs to meet are:

[0204]

[0205] where t tar is the time when the target reaches the target phase. The time when the adjusted vehicle reaches the rendezvous phase is which can be calculated by the following formula:

[0206]

[0207] Then the conditions that the adjusted rendezvous needs to meet are:

[0208]

[0209] Substituting the t tar expression, the adjusted apocenter orbit altitude in the first orbit can be derived as

[0210]

[0211] Therefore, an appropriate h can be selected according to the feasible range of the orbit altitude a1 , thereby determining the entire aerodynamic assist maneuver trajectory. After correcting to obtain the new orbit altitude in the first orbit, the adjusted apocenter altitude sequence can be obtained as [h′ a1 , h a2 ,..., h aN . By solving equation (33) orbit by orbit, the entire aerodynamic assist maneuver trajectory can be obtained.

[0212] S5. Obtain the aerodynamic assist rendezvous maneuver trajectory.

[0213] The specific embodiments are as follows:

[0214] The vehicle model adopts NASA's latest "Orion" multi-purpose deep space manned spacecraft (MPCV) with a mass of 10387 kg and a nominal lift-drag ratio of 0.27. The gravitational constant of Mars μ = 42828 m 3 / s 2 , the radius of Mars is r M = 3396 km, the second-order spherical harmonic coefficient J2 = -1.95545×10 -3 , the atmospheric altitude is 128 km, the Mars atmospheric density index parameter h s = 8805.7 m, ρ0 = 0.01474 kg / m 3 . The initial orbital elements are

[0215] [3×10 4 m, 0.3333, 20°, 45°, 0°, 180°]. The initial orbital elements of the target satellite are taken as those of Phobos, which are [0.94×10 4 m, 0.015, 1.1°, 169.2°, 189°, 180°]. The calculated full - course trajectory diagram is as shown in Figure 2 the figure. Table 1 shows the de - orbit pulse of the aerodynamic - assisted maneuver and the apoapsis maneuver pulse per orbit.

[0216] Table 1 De - orbit pulse of the aerodynamic - assisted maneuver and the apoapsis maneuver pulse per orbit

[0217] Number of turns 0 (off-rail) 1 2 3 (on-rail) Maneuvering pulse (m / s) -895.9265 -0.6203 0.4565 909.2842

[0218] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit them. Although the present invention has been described in detail with reference to the preferred embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions of the present invention or make equivalent replacements, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy design method, characterized in that: The following steps are involved: S1. Establish polar coordinate dynamics model of the Mars probe flying in the Martian atmosphere and flight dynamics model outside the atmosphere; S2, calculate the out-of-plane rendezvous position or the rendezvous phase of a given co-orbital plane; S3, rough matching of intersection phase; S4, precise matching of intersection phase; S5. Obtain an aerodynamically assisted rendezvous maneuver trajectory.

2. The method for designing a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy according to claim 1, characterized in that: The polar coordinate dynamics of the Mars probe flying in the Martian atmosphere in S1 is: Among them, g r , g φ Represent the radial and tangential gravitational acceleration respectively; their expressions are: Among them, r is the distance between the center of mass of the aircraft and the center of the planet; V is the speed; γ is the flight path angle or track angle of the aircraft, ψ is the heading angle of the aircraft; σ is the roll angle of the aircraft, which is the control quantity; θ and φ are longitude and latitude respectively; ω m is the planet's rotation angular velocity; μ is the planet's gravitational constant, r M is the planet radius; J2 is the planet second-order spherical harmonic coefficient; lift acceleration L and drag acceleration D are: Where S is the reference area of ​​the aircraft; C L and C D are lift coefficient and drag coefficient respectively; atmospheric density ρ=ρ0exp(-h / h s ), where ρ0 is the atmospheric density on the planet surface; h s is the atmospheric density coefficient; m is the mass of the aircraft, h is the current altitude of the aircraft, h = rr M .

3. The method for designing a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy according to claim 2, characterized in that: The specific contents of S2 are: For the out-of-plane orbit rendezvous, let the six elements of the initial orbit of the spacecraft be [a1, e1, ω1, Ω1, i1, θ1], which represent the semi-major axis, eccentricity, perigee argument, ascending node right ascension, orbit inclination, and true anomaly respectively; then calculate the normal vector n1 of the orbital plane according to the definition of the orbital elements: n1=[sini1sinΩ1,-sini1cosΩ1,cosi1] (10) Similarly, according to the target orbital elements [a2, e2, ω2, Ω2, i2, θ2], the target orbital plane normal vector n2 is calculated: n2=[sini2sinΩ2,-sini2cosΩ2,cosi2] (11) The cross product of the two orbital plane normal vectors gives the vector n3 of the intersection line of the two orbital planes: n3=n1×n2 (12) Then the rendezvous phase is located at the intersection of the vector on the target orbit, and the intersection is projected to the planetary Cartesian inertial system; when the X-axis of the Cartesian inertial system is rotated in the right-hand direction around the positive direction of the Z-axis until it coincides with the projection, the angle swept is the rendezvous phase; the rendezvous phase is defined in [0,2π]; for any two skew orbits, there are two rendezvous phases to be calculated, and any one of them is taken as the rendezvous phase. In the specific calculation, assume n3=[x3,y3,z3] (13) Assume that the intersection phases obtained according to this definition are ψ1 and ψ2. When x3≥0 and y3>0, ψ1=arctan(y3 / x3) (14) When x3≥0 and y3≤0, ψ1=2π+arctan(y3 / x3) (15) When x3≤0, ψ1=π+arctan(y3 / x3) (16) When calculating ψ2, n3 is regarded as an opposite vector and the calculation is performed again.

4. The method for designing a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy according to claim 3, characterized in that: The S3 content is as follows: First, the first step is to integrate and get the de-orbit state, integrating from the initial position of the spacecraft to the current phase of the spacecraft equal to the target phase; let the initial state of the spacecraft be x0 = [r0, v0], let the cutoff condition of the orbital integration be that the current phase is the target orbital phase; let the phase of the spacecraft be φ a , the initial time is t0, g(x) is the orbital dynamics equation; then the integral to the de-orbit phase is expressed as: Where x a =[r a ,v a ],t a For this time, the given de-orbit orbit height h is calculated by the de-orbit phase. a0 =||r a ||; Assume the height of the apocenter of the target orbit is h afaim , the relationship between the orbital element and the pericentric velocity is: Where, the semi-major axis of the orbit a = h a +h p +2R, orbital eccentricity e=(h a -h p ) / (2a); μ is the planetary gravitational constant, h a is the apocenter height, use h a0 Substitution; When h p =h atm When the maximum value of the pericenter velocity is calculated according to the initial apocenter height, the minimum value of the pericenter velocity is obtained according to the apocenter height of the target orbit. The minimum and maximum values ​​of the pericenter velocity are V p,low ,V p,up , the calculation formula is: Where R represents the radius of the planet, and r M Substitute; h atm Indicates the atmospheric height; Construct the pericentric point sequence [V p1 ,V p2 ,...,V pN ], expressed as: Where ΔV p =(V p,up -V p,low ) / (N-1) is the step length of the pericentric velocity sequence, and N is the given total number of laps. According to the relationship between the orbital root number and the pericentric velocity, formula (17) is used to inversely solve the height sequence [h a1 ,h a2 ,...,h aN ]: Where V pn =V p,low +(n-1)ΔV p ; R represents the radius of the planet, use r M Substitution; Suppose the pulse applied at the apocenter after the nth atmospheric crossing is ΔV n , the starting pulse is ΔV0; r a0 ,v a0 are the initial orbital position and velocity vector respectively, given as the orbital state at the apocenter of an elliptical orbit; are the apocenter position and velocity vector after the n-1th atmospheric crossing and after leaving the atmosphere; the position velocity vector after the pulse is applied is the initial apocenter position vector of the nth circle, recorded as ΔV n The pulse application direction of ΔV0 is The direction is the same, then the position velocity vector after the pulse is applied Satisfy the following formula: Assume that the cutoff condition of the orbital integral is Γ0=||r||-r atm =0; from the initial time t ini , the state vector after the pulse is applied The integration to the atmospheric inlet is expressed as: in are the position vector, velocity vector and time of the atmospheric inlet respectively; g(r,v) is the orbital dynamics differential equation considering two-body perturbation; Transformation of the rv vector from the Cartesian inertial coordinate system to atmospheric flight conditions Set it as the following function: Further by the atmospheric inlet state The exit state is obtained by integrating the atmospheric flight dynamics equations The integral cutoff condition is Γ0 = rr atm = 0, the control amount in the dynamics, the tilt angle σ, is taken as a constant; it is expressed as: Where f(x) is the atmospheric flight dynamics equation as shown in (1)-(6); further through the function Φ x→rv Get the velocity vector of the atmospheric outlet position in the Cartesian inertial coordinate system Finally, the state of the apocenter is obtained through orbital integration. The integral cutoff condition is Γ1=r·v=0, and the initial condition is are the position vector, velocity vector and time before the pulse is applied at the apogee respectively; let the apocenter height that can be reached by the nth circle flying out of the atmosphere be Satisfies the equation: Construct the pulse ΔV applied at the apocenter after the nth atmospheric crossing n Or the initial pulse ΔV0, the apocenter position and velocity vector reached after the n-1th atmospheric crossing and after leaving the atmosphere The velocity vector of the apocenter position reached by the nth circle flying out of the atmosphere Apocenter height The functional relationship is expressed as: The nonlinear mapping relationship between the coordinated pulse and the outlet apocenter is constructed. For the first pulse ΔV0, the mapping is constructed similarly as follows: where r a0 ,v a0 are the position and velocity vectors of the initial orbital position; Assume that the output of formula (30) is When , its mapping is expressed as: The nth element of the height sequence obtained by solving is the target apocenter height of the nth atmospheric crossing. Let the output of (31) be the target apocenter height. In the calculation of the nth circle, It has been obtained from the calculation of the previous cycle, and the ΔV n A univariate nonlinear equation for : h an is the constant calculated in (22); then the solution is obtained by the bisection method, and the initial solution boundary is taken as [-10,10]m / s; The position vector of the apocentric point of the terminal rendezvous of the aerodynamic-assisted maneuver is obtained by solving: The corresponding intersection phase is φ b ; then the phase deviation is expressed as Δφ=φ b -ψ1, assuming that the de-orbit phase should be corrected to φ′ a , and its calculation formula is: f′ a =φ a -Df (34) After calculating the new off-orbit phase, repeat the above steps to obtain the continuously corrected off-orbit phase; set the iterative stop condition to the following formula: |φ′ a -f a |≤δ φ (35) where δ φ is the iteration error limit; after the iteration meets the iteration condition, the phase rough matching is completed.

5. The method for designing a rarefied atmosphere aerodynamic-assisted rendezvous maneuver strategy according to claim 4, characterized in that: The specific contents of S4 are: Assume that the number of generated aerodynamic assist maneuvers is N, and the orbital height during the deorbit maneuver is h a0 After each subsequent atmospheric flight, the apocenter height reached is [h a1 ,h a2 ,...,h aN ], the height of the pericenter is [h p1 ,h p2 ,...,h pN ], then the feasible rendezvous time sequence of the target satellite is calculated below; let the initial orbital state of the target satellite be x t0 =[r t0 ,v t0 ], assuming that the cutoff condition of orbital integration is that the current phase is the target orbital phase; set up The target aircraft phase is φ ta , the initial time is t0, then the integral to the target intersection is expressed as: According to the initial orbital elements of the target satellite [a2, e2, ω2, Ω2, i2, θ2], the orbital period of the target satellite is calculated as: The time series of the target satellite's reachable target phase is: t tar =t ta +kT t ,k=1,2,3,.... (38) Fine-tune the exit apocenter height of the aircraft's first atmospheric flight to achieve time adjustment for aerodynamic-assisted rendezvous; The orbital period of the first circle is calculated based on the apocenter and pericenter sequence of the aerodynamic assist maneuver before adjustment; the orbital period of the first circle before adjustment is: h p1 is the pericentric height; h a1 The adjusted apocenter height is h′ a1 , the feasible range is: h a2 <h′ a1 <h a0 (40) The orbital period of the first circle after adjustment is: Suppose the time it takes for the aircraft to reach the apocenter position after it leaves the atmosphere for the last time is The time when the aircraft is expected to arrive at the target phase, the conditions for rendezvous are: where t tar is the time when the target reaches the target phase; the time when the adjusted aircraft reaches the rendezvous phase is Calculated by the following formula: The conditions satisfied by the adjusted intersection are: t tar Substituting the expression into the equation, we can deduce that the adjusted orbital height of the apocenter point of the first circle is: Choose an appropriate h based on the feasible range of orbital altitudes a1 , determines the entire aerodynamic-assisted maneuver trajectory; after correcting the new first orbit height, the adjusted apocenter height sequence is [h′ a1 ,h a2 ,...,h aN ]; By solving equation (33) circle by circle, the full aerodynamic-assisted maneuver trajectory is obtained.

Images