CN112861250B - Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance - Google Patents

Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance Download PDF

Info

Publication number
CN112861250B
CN112861250B CN202011532177.8A CN202011532177A CN112861250B CN 112861250 B CN112861250 B CN 112861250B CN 202011532177 A CN202011532177 A CN 202011532177A CN 112861250 B CN112861250 B CN 112861250B
Authority
CN
China
Prior art keywords
solution
equation
gliding
aircraft
formula
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011532177.8A
Other languages
Chinese (zh)
Other versions
CN112861250A (en
Inventor
朱恒伟
廖波
江庆平
郭胜鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
23 Units Of Chinese People's Liberation Army 96901 Force
Original Assignee
23 Units Of Chinese People's Liberation Army 96901 Force
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 23 Units Of Chinese People's Liberation Army 96901 Force filed Critical 23 Units Of Chinese People's Liberation Army 96901 Force
Priority to CN202011532177.8A priority Critical patent/CN112861250B/en
Publication of CN112861250A publication Critical patent/CN112861250A/en
Application granted granted Critical
Publication of CN112861250B publication Critical patent/CN112861250B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Fluid Mechanics (AREA)
  • Automation & Control Theory (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
  • Navigation (AREA)

Abstract

The invention belongs to the technical field of guidance, and particularly relates to a gliding trajectory calculation method. The technical scheme is as follows: the method for calculating the degradation solution of the gliding trajectory along with the energy change based on the attack angle and the resistance comprises the following steps: establishing a dynamic model of the gliding aircraft; establishing an approximate equation of a dynamic model of the gliding aircraft containing three-dimensional motion; an analytical solution of the first order approximation equation is solved. Compared with the traditional method, the method has the advantages that the method takes the attack angle and the resistance acceleration as control variables, establishes the three-dimensional gliding aircraft dynamics approximate equation, obviously improves the precision, can quickly solve the large-range and three-dimensional maneuvering trajectory, and can meet the requirements of on-line trajectory planning and quick guidance solving of the missile-borne computer in engineering realization.

Description

Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance
Technical Field
The invention belongs to the technical field of guidance, and particularly relates to a gliding trajectory energy-variation-along-gliding degradation calculation method based on an attack angle and resistance.
Background
The guidance planning of the high-speed gliding trajectory is complex and large in calculation amount, and the traditional numerical calculation method for solving the differential equation is long in calculation time consumption and difficult to meet the requirements of on-line trajectory planning and guidance quick calculation. Because ballistic planning needs to execute a large number of iterative operations, the analytic solutions have great value in the aspect of online ballistic planning, but most of glide ballistic analytic solutions in the prior art are difficult to meet the engineering application requirements of online ballistic planning and guidance quick solving due to the defects in the aspects of calculation speed and precision when processing three-dimensional maneuvering ballistics with large-range transverse motion, so that a cost reduction solution which can quickly solve the large-range three-dimensional maneuvering ballistics and can meet the engineering application requirements in precision is necessary to be researched.
Disclosure of Invention
The invention aims to provide a gliding trajectory degradation solution of an attack angle and a resistance acceleration along with energy change, the attack angle and the resistance acceleration are used as control variables, the requirements of high-speed gliding trajectory rapid guidance solution and online trajectory planning are met, compared with the traditional analysis solution, the gliding trajectory degradation solution can rapidly solve a large-range three-dimensional maneuvering trajectory, and the calculation speed and the calculation precision can meet the engineering application requirements of online guidance rapid solution.
The technical scheme adopted by the invention is as follows:
the gliding trajectory degradation solution calculation method along with the energy change based on the attack angle and the resistance comprises the following steps of:
step1: establishing a dynamic model of a gliding aircraft
Defining: initial longitude, initial latitude and initial azimuth angle of gliding aircraft in geocentric coordinate system are respectively lambda 0 、φ 0 And alpha 0 ；
Defining: the offset geocentric coordinate system rotates the geocentric coordinate system around the z, y and x axes by a rotating angle lambda respectively 0 、-φ 0 And alpha 0 Obtaining a coordinate system;
defining: the longitude, the latitude and the heading angle of the gliding aircraft under the offset geocentric coordinate system are respectively lambda, phi and psi;
introducing dimensionless heights
Figure GSB0000200765950000011
Wherein->
Figure GSB0000200765950000012
H m For the average flight height, R, of the gliding aircraft in the gliding flight section e Is the radius of the earth, r e The distance between the center of mass of the aircraft and the geocentric; let E 0 And E f The energy of the aircraft at the gliding starting point and the gliding ending point respectively;
establishing a dynamic model of the gliding aircraft with the energy E as an independent variable:
Figure GSB0000200765950000021
wherein A is D For resisting acceleration, L Dy Is the longitudinal lift-to-drag ratio, L Dz Is the transverse lift-to-drag ratio, omega e Is the angular velocity of rotation of the earth, mu is the gravitational constant of the earth;
Figure GSB0000200765950000022
C σ ≈2vω e (sinφ-sinψtanθcosφ)
C θ ≈-2vω e cosφcosψ
step2: approximate equation for establishing dynamics model of gliding aircraft based on Newton iteration method
Let x E = E, variable substitution for equation (1), and rewrite to y' (x) E )＝f(x E ，y)，
Wherein, the first and the second end of the pipe are connected with each other,
Figure GSB0000200765950000023
let y i Is y' (x) E )＝f(x E Y) and setting the i +1 th iteration value as y i+1 ＝y i +δy i Function y' (x) E )＝f(x E Y) is constructed as a newton iteration in function space, i.e.:
Figure GSB0000200765950000031
wherein:
Figure GSB0000200765950000032
x Eini 、/>
Figure GSB0000200765950000033
θ ini 、φ ini 、ψ ini 、λ ini initial conditions for compatibility;
knowing the value y of the ith iteration i Then, the solution y of the (i + 1) th step is obtained by the formula (2) i+1 ；
Defining an equation formed by using the formula (2) at the Nth time as an approximate equation of the formula (1) at the Nth time, and solving the approximate equation into an approximate solution of the Nth time; definition of
Figure GSB0000200765950000034
θ 0 ，φ 0 ，ψ 0 ，λ 0 And &>
Figure GSB0000200765950000035
θ 1 ，φ 1 ，ψ 1 ，λ 1 Respectively representing a zero-order approximate solution and a first-order approximate solution of formula (1);
when the offset geocentric system equation is adopted, the gliding aircraft approximately flies to the target along the equator of the offset geocentric coordinate system, and the gliding flight section has phi approximately equal to 0, theta approximately equal to 0, psi approximately equal to 0,
Figure GSB0000200765950000036
The characteristic of (1) is that the zero-order approximate solution is:
Figure GSB0000200765950000037
substituting the equation (3) into the equation (2) to obtain a first order approximation equation of the glide flight segment as follows:
Figure GSB0000200765950000038
wherein
Figure GSB0000200765950000039
Figure GSB00002007659500000310
As can be seen from the equation (4), the Jacobian matrix has the characteristic of block sparsity, and theta can be solved first 1 And
Figure GSB0000200765950000041
and phi 1 And psi 1 Then solve for λ 1 ；
Due to phi 1 And psi 1 Is the solution of the lateral motion variable,
Figure GSB0000200765950000042
and theta 1 Is a solution of a longitudinal motion variable, phi 1 And psi 1 To (X)>
Figure GSB0000200765950000043
And theta 1 Insensitive, equation (4) is further simplified as: />
Figure GSB0000200765950000044
Step3: analytic solution to solve first order approximation equation
Step31: using the attack angle and the resistance acceleration as control variables to solve lambda 0 Analytic solution of
Order to
Figure GSB0000200765950000045
Wherein, c 0 ，c 1 ，c 2 Are all controlled variables resistance acceleration A D The design parameters of (1); alpha is an attack angle, f (E) is an attack angle control function designed by taking energy as independent variable;
defining the flight altitude of the gliding aircraft as h, and respectively solving the first derivative of the aerodynamic drag acceleration and the aerodynamic drag speed with respect to energy to obtain:
Figure GSB0000200765950000046
finishing to obtain:
Figure GSB0000200765950000047
wherein, C D Is the drag coefficient of the gliding aircraft, h s The coefficients of an exponential formula being a function of atmospheric density;
neglecting the influence of the drag coefficient derivative, the approximate height h of the gliding aircraft * Comprises the following steps:
Figure GSB0000200765950000051
and (3) solving a second derivative related to energy for the resistance acceleration, and solving a longitudinal lift-drag ratio by using a speed dip angle theta equation in the joint type (1):
L Dy ＝a(A D ″-b)
wherein, the first and the second end of the pipe are connected with each other,
Figure GSB0000200765950000052
obtaining transverse lift-drag ratio by using pneumatic coefficient table or pneumatic interpolation function
Figure GSB0000200765950000053
/>
Wherein
Figure GSB0000200765950000054
The sign of which is determined by determining the azimuth deviation.
Formula (6) is substituted for formula (3) and from E 0 Accumulating to E to obtain:
Figure GSB0000200765950000055
step32: solution of phi 1 ，ψ 1 ，λ 1 Analytic solution of
The first two formulas of the united vertical type (5) are obtained
Figure GSB0000200765950000056
Order to
Figure GSB0000200765950000057
Left-multiplying equation (7) by M (E, E) 0 ) After arrangement, the differential multiplication rule is applied reversely, and then integration and arrangement are carried out to obtain:
Figure GSB0000200765950000058
using M (E, E) 0 ) Solving its inverse matrix [ M (E, E) 0 )] -1 Bringing into formula (8) order
Figure GSB0000200765950000061
Finish to obtain phi 1 And psi 1 The analytical formula (2):
Figure GSB0000200765950000062
lagrange interpolation polynomial interval [ E ] using root of Legendre of n times as sampling node f ，E 0 ]Approximate λ (E, E) on segment 0 )、cos(λ(x E ，E 0 ))m 3 (x E )、sin(λ(x E ，E 0 ))m 3 (x E ) Obtaining:
Figure GSB0000200765950000063
Figure GSB0000200765950000064
Figure GSB0000200765950000065
then solve to obtain phi 1 And psi 1 Analytic solution of (2):
Figure GSB0000200765950000066
c p0(k) 、c p1(k) 、c p2(k) is the coefficient of the interpolation polynomial; by calculated phi 1 And psi 1 Substitution of formula (5) to give lambda 1 Analytic solution of (2):
Figure GSB0000200765950000067
the invention has the beneficial effects that:
the invention establishes a gliding trajectory degradation solution of the attack angle and the resistance acceleration along with the energy change, realizes the rapid calculation of the hypersonic three-dimensional trajectory, has good calculation speed and precision in the aspect of calculating the large-range three-dimensional maneuvering trajectory compared with the traditional analytic solution, and can meet the application requirement of the online guidance planning engineering of the large-range three-dimensional gliding aircraft.
Drawings
None.
Detailed Description
The technical solution of the present invention is further specifically described below with reference to specific embodiments.
The method for calculating the degradation solution of the gliding trajectory along with the energy change based on the attack angle and the resistance comprises the following steps:
step1: establishing a dynamic model of a gliding aircraft
Defining: initial longitude, initial latitude and initial azimuth angle of gliding aircraft in geocentric coordinate system are respectively lambda 0 、φ 0 And alpha 0 ；
Defining: the offset geocentric coordinate system rotates the geocentric coordinate system around the z, y and x axes by a rotating angle lambda respectively 0 、-φ 0 And alpha 0 Obtaining a coordinate system;
defining: the longitude, the latitude and the heading angle of the gliding aircraft under the offset geocentric coordinate system are respectively lambda, phi and psi;
taking into account that the flying height during gliding flight is a small quantity relative to the radius of the earth, dimensionless heights are introduced
Figure GSB0000200765950000071
Wherein->
Figure GSB0000200765950000072
H m For the average flight height, R, of the gliding aircraft in the gliding flight section e Is the radius of the earth, r e The distance between the center of mass of the aircraft and the geocentric; let E 0 And E f The energy of the aircraft at the gliding starting point and the gliding ending point respectively;
establishing a dynamic model of the gliding aircraft with the energy E as an independent variable (the model introduces dimensionless altitude for the classical dynamic equation of the gliding aircraft)
Figure GSB0000200765950000073
To obtain): />
Figure GSB0000200765950000074
Wherein A is D For resisting acceleration, L Dy Is the longitudinal lift-to-drag ratio, L Dz Is the transverse lift-drag ratio, omega e Is the angular velocity of rotation of the earth, mu is the gravitational constant of the earth;
Figure GSB0000200765950000081
C σ ≈2vω e (sinφ-sinψtanθcosφ)
C θ ≈-2vω e cosφcosψ
step2: approximate equation for establishing dynamics model of gliding aircraft based on Newton iteration method
Let x be E = E, variable replacement write y' (x) to equation (1) E )＝f(x E ，y)，
Wherein the content of the first and second substances,
Figure GSB0000200765950000082
a Newton iteration method is applied to approximate the multi-dimensional nonlinear trajectory equation into several low-dimensional variable coefficient linear differential equations so as to realize high-precision flight trajectory fast calculation;
let y i Is y' (x) E )＝f(x E Y) and setting the i +1 th iteration value as y i+1 ＝y i +δy i Function y' (x) E )＝f(x E Y) is constructed as a newton iteration in function space, i.e.:
Figure GSB0000200765950000083
wherein:
Figure GSB0000200765950000084
x Eini 、/>
Figure GSB0000200765950000085
θ ini 、φ ini 、ψ ini 、λ ini initial conditions for compatibility;
knowing the value y of the ith iteration i Then, the solution y of the (i + 1) th step is obtained by the formula (2) i+1 ，y ini Initial conditions for compatibility;
to solve equation (2), a zero-order approximation solution of equation (1) needs to be provided before starting the iterative computation. The closer the zero solution is to the true solution, the fewer iterations are required to achieve the same accuracy requirement. For aircraft guidance problems, a large ballistic estimation error is usually allowed, if the zero-time solution is closer to the true solution, one solution can meet the requirement,
defining an equation formed by using the formula (2) at the Nth time as an approximate equation of the formula (1) at the Nth time, and solving the approximate equation into an approximate solution of the Nth time; definition of
Figure GSB0000200765950000091
θ 0 ，φ 0 ，ψ 0 ，λ 0 And &>
Figure GSB0000200765950000092
θ 1 ，φ 1 ，ψ 1 ，λ 1 Respectively representing a zero-order approximate solution and a first-order approximate solution of formula (1);
c used for calculating the Jacobian matrix in order to obtain the Jacobian matrix which is accurate and has the partitioned sparse characteristic of decoupling the equation θ 、C σ In the expression, only ω is ignored e 2 And J2, retention of ω e The first term of (a);
when the offset geocentric system equation is adopted, the gliding aircraft approximately flies to the target along the equator of the offset geocentric coordinate system, and the gliding flight section has phi approximately equal to 0, theta approximately equal to 0, psi approximately equal to 0,
Figure GSB0000200765950000093
The characteristic of (1) is that the zero-order approximate solution is:
Figure GSB0000200765950000094
substituting the equation (3) into the equation (2) to obtain a first order approximation equation of the glide flight segment as follows:
Figure GSB0000200765950000095
wherein
Figure GSB0000200765950000096
Figure GSB0000200765950000097
Description of the invention: the equation (4) reserves the basic characteristic of the equation (1) and contains a first term of the earth rotation rate, so that the accuracy of the 1-time approximate solution is obviously improved compared with the accuracy of the 0-time approximate solution, and the 1-time approximate solution is used for guidance calculation to obtain better guidance performance. Equation (2) can be used multiple times to obtain a more accurate solution if desired.
As can be seen from the equation (4), the Jacobian matrix has the characteristic of block sparsity, and theta can be solved first 1 And
Figure GSB0000200765950000098
and phi 1 And psi 1 Then solve for λ 1 ；
Due to phi 1 And psi 1 Is a solution to the lateral motion variable and,
Figure GSB0000200765950000101
and theta 1 Is a solution of a longitudinal motion variable, phi 1 And psi 1 Is paired and/or matched>
Figure GSB0000200765950000102
And theta 1 Insensitive, equation (4) is further simplified as:
Figure GSB0000200765950000103
step3: analytic solution to solve first order approximation equation
Step31: using the attack angle and the resistance acceleration as control variables to solve lambda 0 Analytic solution of
Order to
Figure GSB0000200765950000104
Wherein, c 0 ，c 1 ，c 2 Are all controlled variables resistance acceleration A D The design parameters of (1); alpha is an attack angle, f (E) is an attack angle control function designed by taking energy as independent variable;
defining the flight altitude of the gliding aircraft as h, and respectively solving the first derivative of the aerodynamic drag acceleration and the aerodynamic drag speed with respect to energy to obtain:
Figure GSB0000200765950000105
finishing to obtain:
Figure GSB0000200765950000106
wherein, C D Is the drag coefficient of the gliding aircraft, h s The coefficients of the exponential formula being a function of atmospheric density;
neglecting the influence of the drag coefficient derivative, the approximate height h of the gliding aircraft * Comprises the following steps:
Figure GSB0000200765950000111
solving a second derivative related to energy for the resistance acceleration, and solving a longitudinal lift-drag ratio by a speed dip angle theta equation in the joint type (1):
L Dy ＝a(A D ″-b)
wherein the content of the first and second substances,
Figure GSB0000200765950000112
obtaining transverse lift-drag ratio by using pneumatic coefficient table or pneumatic interpolation function
Figure GSB0000200765950000113
Wherein
Figure GSB0000200765950000114
The sign of which is determined by judging the azimuth deviation;
formula (6) is substituted for formula (3) and from E 0 Accumulating to E to obtain:
Figure GSB0000200765950000115
step32: solution of phi 1 ，ψ 1 ，λ 1 Analytic solution of
The first two formulas of the united vertical type (5) are obtained
Figure GSB0000200765950000116
Order to
Figure GSB0000200765950000117
Left-multiplying equation (7) by M (E, E) 0 ) After sorting, the differential multiplication principle is applied reversely, and then integration and sorting are carried out to obtain:
Figure GSB0000200765950000118
using M (E, E) 0 ) Solving its inverse matrix [ M (E, E) 0 )] -1 Bringing into formula (8) order
Figure GSB0000200765950000121
Finish to obtain phi 1 And psi 1 The analytical formula (2):
Figure GSB0000200765950000122
lagrange interpolation polynomial interval [ E ] using root of Legendre of n times as sampling node f ，E 0 ]Approximate λ (E, E) on segment 0 )、cos(λ(x E ，E 0 ))m 3 (x E )、sin(λ(x E ，E 0 ))m 3 (x E ) Obtaining:
Figure GSB0000200765950000123
Figure GSB0000200765950000124
/>
Figure GSB0000200765950000125
then solve to get phi 1 And psi 1 Analytic solution of (2):
Figure GSB0000200765950000126
c p0(k) 、c p1(k) 、c p2(k) is the coefficient of the interpolating polynomial; by calculated phi 1 And psi 1 Substitution of formula (5) to give lambda 1 Analytic solution of (2):
Figure GSB0000200765950000127
the above detailed description of the present invention is only used for illustrating the present invention and is not limited to the technical solutions described in the embodiments of the present invention, and it should be understood by those skilled in the art that the present invention can be modified or substituted equally to achieve the same technical effects; as long as the use requirements are met, the method is within the protection scope of the invention.

Claims (1)

1. The method for calculating the degradation solution of the gliding trajectory along with the change of energy based on the attack angle and the resistance is characterized by comprising the following steps of:
step1: establishing a dynamics model of a gliding aircraft
Defining: initial longitude, initial latitude and initial azimuth angle of gliding aircraft in geocentric coordinate system are respectively lambda 0 、φ 0 And alpha 0 ；
Defining: offset geocentric coordinate system is geocentric coordinateAre rotated about the z, y and x axes by an angle lambda, respectively 0 、-φ 0 And alpha 0 Obtaining a coordinate system;
defining: the longitude, the latitude and the heading angle of the gliding aircraft under the offset geocentric coordinate system are respectively lambda, phi and psi;
introducing dimensionless heights
Figure FSB0000200765940000011
Wherein->
Figure FSB0000200765940000012
r e ＝R e +H m ，H m Average flying height, R, of gliding aircraft in gliding flight section e Is the radius of the earth, r e The distance between the center of mass of the aircraft and the center of the earth; let E 0 And E f The energy of the aircraft at the glide starting point and the glide ending point respectively;
establishing a dynamic model of the gliding aircraft with the energy E as an independent variable:
Figure FSB0000200765940000013
wherein A is D For resisting acceleration, L Dy Is the longitudinal lift-to-drag ratio, L Dz Is the transverse lift-drag ratio, omega e Is the angular velocity of rotation of the earth, mu is the gravitational constant of the earth;
Figure FSB0000200765940000014
C σ ≈2vω e (sinφ-sinψtanθcosφ)
C θ ≈-2vω e cosφcosψ
step2: approximate equation for establishing dynamics model of gliding aircraft based on Newton iteration method
Let x E = E, variable substitution for equation (1), and rewrite to y' (x) E )＝f(x E ，y)，
Wherein the content of the first and second substances,
Figure FSB0000200765940000021
let y i Is y' (x) E )＝f(x E Y) and setting the i +1 th iteration value as y i+1 ＝y i +δy i Function y' (x) E )＝f(x E Y) is constructed as a newton iteration in function space, i.e.:
Figure FSB0000200765940000022
wherein:
Figure FSB0000200765940000023
x Eini 、/>
Figure FSB0000200765940000024
θ ini 、φ ini 、ψ ini 、λ ini initial conditions for compatibility;
knowing the value y of the ith iteration i Then, the solution y of the (i + 1) th step is obtained by the formula (2) i+1 ；
To solve equation (2), a zero-order approximate solution of equation (1) needs to be provided before starting the iterative computation. The closer the zero solution is to the true solution, the fewer iterations are required to achieve the same accuracy requirement. For aircraft guidance problems, a large ballistic estimation error is usually allowed, if the zero-time solution is closer to the true solution, one solution can meet the requirement,
defining an equation formed by using the formula (2) at the Nth time as an approximate equation of the formula (1) at the Nth time, and solving the approximate equation into an approximate solution of the Nth time; definition of
Figure FSB0000200765940000025
θ 0 ，φ 0 ，ψ 0 ，λ 0 And &>
Figure FSB0000200765940000026
θ 1 ，φ 1 ，ψ 1 ，λ 1 Respectively representing a zero-order approximate solution and a first-order approximate solution of formula (1);
when the offset geocentric system equation is adopted, the gliding aircraft approximately flies to the target along the equator of the offset geocentric coordinate system, and the gliding flight section has phi approximately equal to 0, theta approximately equal to 0, psi approximately equal to 0,
Figure FSB0000200765940000027
The characteristic of (1) is that the zero-order approximate solution is:
Figure FSB0000200765940000031
substituting the equation (3) into the equation (2) to obtain a first order approximation equation of the gliding flight section as follows:
Figure FSB0000200765940000032
wherein
Figure FSB0000200765940000033
Figure FSB0000200765940000034
As can be seen from the equation (4), the Jacobian matrix has the characteristic of block sparsity, and theta can be solved first 1 And
Figure FSB0000200765940000035
and phi 1 And psi 1 Then solve for λ 1 ；
Due to phi 1 And psi 1 Is the solution of the lateral motion variable,
Figure FSB0000200765940000036
and theta 1 Is a solution of a longitudinal motion variable, phi 1 And psi 1 To (X)>
Figure FSB0000200765940000037
And theta 1 Insensitive, equation (4) is further simplified as:
Figure FSB0000200765940000038
step3: analytic solution to solve first order approximation equation
Step31: using the attack angle and the resistance acceleration as control variables to solve lambda 0 Analytic solution of
Order to
Figure FSB0000200765940000041
Wherein, c 0 ，c 1 ，c 2 Are all controlled variable resistance acceleration A D The design parameters of (1); alpha is an attack angle, f (E) is an attack angle control function designed by taking energy as independent variable;
defining the flight height of the gliding aircraft as h, and respectively solving the first derivative of the aerodynamic drag acceleration and the aerodynamic drag speed with respect to energy to obtain:
Figure FSB0000200765940000042
finishing to obtain:
Figure FSB0000200765940000043
wherein, C D Is the drag coefficient of the gliding aircraft, h s System of exponential formula as function of atmospheric densityCounting;
neglecting the influence of the drag coefficient derivative, the approximate height h of the gliding aircraft * Comprises the following steps:
Figure FSB0000200765940000044
and (3) solving a second derivative related to energy for the resistance acceleration, and solving a longitudinal lift-drag ratio by using a speed dip angle theta equation in the joint type (1):
L Dy ＝a(A D ″-b)
wherein the content of the first and second substances,
Figure FSB0000200765940000045
obtaining transverse lift-drag ratio by using pneumatic coefficient table or pneumatic interpolation function
Figure FSB0000200765940000046
Wherein
Figure FSB0000200765940000051
The sign of which is determined by determining the azimuth deviation.
Formula (6) is substituted for formula (3) and from E 0 Accumulating to E to obtain:
Figure FSB0000200765940000052
step32: solution of phi 1 ，ψ 1 ，λ 1 Analytic solution of
The first two formulas of the united vertical type (5) are obtained
Figure FSB0000200765940000053
Order to
Figure FSB0000200765940000054
Left-multiplying equation (7) by M (E, E) 0 ) After arrangement, the differential multiplication rule is applied reversely, and then integration and arrangement are carried out to obtain:
Figure FSB0000200765940000055
using M (E, E) 0 ) Solving its inverse matrix [ M (E, E) 0 )] -1 Bringing into formula (8) order
Figure FSB0000200765940000056
Finish to obtain phi 1 And psi 1 The analytical formula (2):
Figure FSB0000200765940000057
lagrange interpolation polynomial interval [ E ] using root of Legendre of n times as sampling node f ，E 0 ]Approximate λ (E, E) on segment 0 )、cos(λ(x E ，E 0 ))m 3 (x E )、sin(λ(x E ，E 0 ))m 3 (x E ) And obtaining:
Figure FSB0000200765940000058
Figure FSB0000200765940000059
Figure FSB0000200765940000061
then solve to get phi 1 And psi 1 Analytic solution of (2):
Figure FSB0000200765940000062
c p0(k) 、c p1(k) 、c p2(k) is the coefficient of the interpolation polynomial; by calculated phi 1 And psi 1 Substitution of formula (5) to give lambda 1 Analytic solution of (2):
Figure FSB0000200765940000063
/>
CN202011532177.8A 2020-12-21 2020-12-21 Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance Active CN112861250B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011532177.8A CN112861250B (en) 2020-12-21 2020-12-21 Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011532177.8A CN112861250B (en) 2020-12-21 2020-12-21 Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance

Publications (2)

Publication Number Publication Date
CN112861250A CN112861250A (en) 2021-05-28
CN112861250B true CN112861250B (en) 2023-03-28

Family

ID=75996271

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011532177.8A Active CN112861250B (en) 2020-12-21 2020-12-21 Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance

Country Status (1)

Country Link
CN (1) CN112861250B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116992553B (en) * 2023-05-25 2024-02-06 中国人民解放军32804部队 Whole-course trajectory estimation method of boosting gliding aircraft

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107941087A (en) * 2017-10-18 2018-04-20 北京航空航天大学 A kind of superb steady gliding reentry guidance method of high lift-drag ratio based on resistance profiles
CN109740198A (en) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 Analytic prediction-based three-dimensional reentry guidance method for gliding aircraft
CN111306989A (en) * 2020-03-12 2020-06-19 北京航空航天大学 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107941087A (en) * 2017-10-18 2018-04-20 北京航空航天大学 A kind of superb steady gliding reentry guidance method of high lift-drag ratio based on resistance profiles
CN109740198A (en) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 Analytic prediction-based three-dimensional reentry guidance method for gliding aircraft
CN111306989A (en) * 2020-03-12 2020-06-19 北京航空航天大学 Hypersonic velocity reentry guidance method based on steady glide trajectory analytic solution

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
防区外投放制导炸弹滑翔段垂直面最优制导;侯明善;《兵工学报》;20070515(第05期);全文 *

Also Published As

Publication number Publication date
CN112861250A (en) 2021-05-28

Similar Documents

Publication Publication Date Title
CN111783307B (en) Hypersonic aircraft state estimation method
CN109614633A (en) A kind of composite rotor craft non-linear modeling method and Calculate Ways
CN112668104B (en) Online identification method for pneumatic parameters of hypersonic aircraft
CN109446582B (en) High-precision order-reduction steady gliding dynamics modeling method considering earth rotation
CN109240323B (en) Aerospace vehicle reentry guidance method capable of analyzing structure in real time
CN107330152B (en) Efficient pneumatic balancing method suitable for rotor craft
Zhu et al. Impact time and angle control guidance independent of time-to-go prediction
CN111290278B (en) Hypersonic aircraft robust attitude control method based on prediction sliding mode
CN109062055A (en) A kind of Near Space Flying Vehicles control system based on Back-stepping robust adaptive dynamic surface
CN112861250B (en) Method for calculating degradation solution of gliding trajectory along with energy change based on attack angle and resistance
CN104536448B (en) Backstepping based control method for unmanned-plane attitude system
CN111008488B (en) Propeller unmanned aerial vehicle launching process reaction torque modeling method
Ben-Asher et al. New proportional navigation law for ground-to-air systems
CN112198885A (en) Unmanned aerial vehicle control method capable of meeting autonomous landing requirement of maneuvering platform
CN112507467B (en) Method for calculating descending order solution of glide trajectory along with speed change based on resistance and lift-drag ratio
CN107796401B (en) Skip reentry vehicle linear pseudo-spectrum parameter correction transverse guidance method
CN112861249B (en) Method for calculating degradation solution of gliding trajectory along with speed change based on attack angle and resistance
CN109445283B (en) Control method for fixed-point tracking of under-actuated aerostat on plane
CN112507465B (en) Gliding trajectory energy change degradation solution calculation method based on resistance and lift-drag ratio
CN112861251B (en) Glide trajectory degradation solution calculation method with speed change considered by earth rotation rate
CN112507466B (en) Method for calculating glide trajectory energy-dependent degradation solution considering earth rotation rate
CN110750053A (en) Error analysis method for semi-physical simulation system of aircraft
CN102508819B (en) Angular-speed-based quaternion Legendre approximate output method during extreme flying of aircraft
CN104166401A (en) Method for calculating single-sliding-block moving mass controlled aircraft balance motion state
Kim et al. Flight dynamics analyses of a propeller-driven airplane (I): aerodynamic and inertial modeling of the propeller

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant