CN113935176B - Efficient dynamics modeling method for electrodynamic force rope derailing device - Google Patents

Efficient dynamics modeling method for electrodynamic force rope derailing device Download PDF

Info

Publication number
CN113935176B
CN113935176B CN202111224605.5A CN202111224605A CN113935176B CN 113935176 B CN113935176 B CN 113935176B CN 202111224605 A CN202111224605 A CN 202111224605A CN 113935176 B CN113935176 B CN 113935176B
Authority
CN
China
Prior art keywords
unit
axis
electrodynamic
matrix
force rope
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
CN202111224605.5A
Other languages
Chinese (zh)
Other versions
CN113935176A (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202111224605.5A priority Critical patent/CN113935176B/en
Publication of CN113935176A publication Critical patent/CN113935176A/en
Application granted granted Critical
Publication of CN113935176B publication Critical patent/CN113935176B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a high-precision and high-efficiency dynamic modeling method for an electrodynamic force rope, and belongs to the field of dynamics and control of orbits and postures of spacecrafts. The implementation method of the invention comprises the following steps: in the electrodynamic force rope system, the main star and the sub-star are regarded as mass points, and the conductive rope is scattered into a plurality of rigid rod units, namely a scattered body model hinged by a plurality of rigid rods; a discrete electrodynamic force rope model considering the coupling effect of multiple physical fields is established by combining a Kane equation and a multi-body dynamics recursive algorithm, the characteristic that the computational complexity of the recursive algorithm is in linear relation with the freedom degree of an electrodynamic force rope system is utilized, the problem of numerical simulation calculation amount increase caused by the increase of discrete units in the discrete electrodynamic force rope model is reduced as much as possible on the premise of ensuring the unchanged precision, and the high-precision and high-efficiency dynamics modeling of the electrodynamic force rope is realized. The method is applied to the field of dynamics and control of the orbit and the attitude of the spacecraft, and can solve the technical problem of the related engineering of the orbit and the attitude of the spacecraft.

Description

Efficient dynamic modeling method for electrodynamic force rope derailing device
Technical Field
The invention relates to a high-precision and high-efficiency dynamic modeling method for an electrodynamic force rope, relates to system attitude control and orbit transfer, and belongs to the field of dynamics and control of orbits and attitudes of spacecrafts.
Background
The number of space debris on the near-earth orbit has shown a tendency to rise significantly in recent years. Space debris not only occupies limited orbital resources, but also can potentially threaten the on-orbit operation of the spaceflight and affect the subsequent spaceflight activities, so that the space debris must be actively removed. Among the numerous active scavenging techniques, electrodynamic ropes are of interest because of their low fuel consumption and low cost advantages.
The electrodynamic force rope is a typical strong nonlinear multi-field coupling system, and with the gradual and deep related research, researchers are concerned about how to establish a high-precision electrodynamic force rope dynamics model in order to more accurately describe the dynamics characteristics of the electrodynamic force rope system during track transfer. Because the electric power rope system has long track transfer time and is highly dependent on the space environment, a high-precision environment model is required to be established, and the method mainly comprises the following steps: an earth gravitational field model, a geomagnetic field model, an atmospheric density model, and a plasma density model. In addition, the conductive tether model is also the key for establishing a high-precision dynamic model, and the model mainly comprises the following components: dumbbell models, discrete body models, and continuum models. Because the dumbbell model and the continuum model are respectively over-simplified or complex, and the discrete body model can further improve the model accuracy by increasing discrete units, the discrete body model is adopted for modeling in most of the current researches, more discrete units are needed for improving the model accuracy, the numerical simulation calculation amount and the simulation time length are inevitably increased, and the research and development period is prolonged.
Disclosure of Invention
In order to solve the problem of high-precision high-efficiency dynamic modeling of the electrodynamic force rope, the invention aims to provide a high-precision high-efficiency dynamic modeling method of the electrodynamic force rope. The method is applied to the field of dynamics and control of the orbit and the attitude of the spacecraft, and can solve the technical problem of the related engineering of the orbit and the attitude of the spacecraft.
The dynamics and control field of the spacecraft orbit and the attitude comprises the stability analysis of the attitude of an electric power rope system, the attitude angle control of the electric power rope system and the orbit transfer of the electric power rope system.
The purpose of the invention is realized by the following technical scheme.
The invention discloses a high-precision and high-efficiency dynamic modeling method of an electrodynamic force rope, which comprises the following steps:
the method comprises the following steps: in the electrodynamic force tether system, a main satellite and a sub-satellite are taken as mass points, a conductive tether is discretized into a plurality of rigid rod units, namely a discretization model hinged by a plurality of rigid rods is adopted, the main satellite is recorded as a unit 0, a conductive tether unit close to the main satellite is recorded as a unit 1, then the number of the conductive tether units is sequentially increased to N, the sub-satellite is recorded as a unit N +1, the units are connected through a spherical hinge, and a relative attitude angle between the units is selected as a generalized coordinate.
Step two: and (3) deriving a kinetic equation of a single unit by using a Kane equation, wherein the single unit refers to the main star, the conductive tether unit and the subsatellite defined in the step one.
In modeling of dynamics of an electrodynamic tether system, coordinate systems OXYZ, OX ' Y ' Z ', OX are definedoyozo，oxuyuzu，oxdydzd，okxkykzkRepresenting the inertial system in turn
Figure BDA0003310739940000024
Earth-following coordinate system, orbit coordinate system
Figure BDA0003310739940000025
North-East-Up (NEU for short) coordinate System
Figure BDA0003310739940000026
North-East-Down (NED for short) coordinate system
Figure BDA0003310739940000027
Unit k body coordinate system
Figure BDA0003310739940000028
. Wherein the inertia system
Figure BDA0003310739940000029
The origin O is located at the center of mass of the earth, the OX axis points to the spring equinox, the OZ axis is consistent with the rotation axis of the earth, and the OY axis and the other two axes form a right-hand system; earth follow-up coordinate system
Figure BDA00033107399400000210
The origin is O, the OZ 'axis is consistent with the OZ axis, and the OX' axis and the OY 'axis rotate around the OZ' axis at the same time at the earth rotation angular speed; orbital coordinate system
Figure BDA00033107399400000211
The origin o is located at the centroid of the main star, oxoThe axis pointing along the earth's centre O to O, ozoThe axis is vertical to the orbital plane and is consistent with the orbital angular momentum direction of the main satellite, oyoThe shaft and the other two shafts form a right-hand system; NEU coordinate system
Figure BDA00033107399400000212
The origin is o, oxuThe axis pointing along point O, oyuThe axis is vertical to the meridian plane of the local area and points to the east, and the oz' axis and the other two axes form a right-handed system; NED coordinate system
Figure BDA00033107399400000213
The origin is o, ozdThe axis pointing along point O, oydThe axis is perpendicular to the meridian plane of the local region and points to east, oxdThe shaft and the other two shafts form a right-hand system; conductive tether unit k body coordinate system
Figure BDA00033107399400000214
Origin okAnd when the whole conductive tether is consistent with the local vertical line of the main star, the unit k body coordinate system is superposed with the track system. Selecting an attitude angle theta between the unitsk,k-1＝[θk,k-1z θk,k-1yθk,k-1x]TAs a generalized variable:
q＝[θ1,o θ2,1 … θN,N-1]T
and defines the attitude angle of the unit k with respect to the unit k-1 as thetak,k-1zAnd thetak,k-1yWherein thetak,k-1zIs xkAxis is at okxkykProjection of plane and xk-1Angle of axis thetak,k-1yIs xkAxis is at okxkykProjection of plane and xkAngle of axis thetak,k-1xFor the last attitude angle and specifying
Figure BDA00033107399400000215
To
Figure BDA00033107399400000216
The rotation sequence of the star is 3-2-1, and the subsatellite is fixedly connected with the conductive tether unit N.
Based on the algorithm of the vector array, dynamic modeling is carried out on the single unit by adopting a Kane equation, and the translation and rotation dynamic equations of the conductive tether unit k are obtained through derivation and arrangement:
Figure BDA0003310739940000021
in the formula mkMass of electrically conductive tether element k, JkIs unit k relative to okAt moment of inertia of
Figure BDA00033107399400000217
Description of (1), SkIs unit k relative to okStatic moment of
Figure BDA00033107399400000218
"x" represents a diagonally symmetric matrix of vectors, akIs a unit k at
Figure BDA00033107399400000219
Medium translational acceleration, ωkIs unit k is opposite
Figure BDA00033107399400000222
At an angular velocity of
Figure BDA00033107399400000223
The component (b) of (a) is,
Figure BDA0003310739940000022
in order to be able to take into account the angular acceleration,
Figure BDA0003310739940000023
is composed of
Figure BDA00033107399400000220
To
Figure BDA00033107399400000221
Of the rotation matrix fk e、fk,ck、fk+1,ckRespectively represent the environmental external force acting on the unit k, the constraint force of the hinge k to the unit k and the constraint force of the hinge k +1 to the unit kThe external force is shown in
Figure BDA00033107399400000311
In, Tk eRepresenting the ambient external moment, T, acting on unit kk+1,ckRepresents the constraint moment of the hinge k +1 to the unit k, which are all represented in
Figure BDA00033107399400000312
In (1).
The external forces acting on the electrodynamic tether system in low-rail environments are mainly composed of earth's gravity, lorentz forces, and atmospheric drag. Using the gravity acceleration at the unit k mass center, the geomagnetic field magnetic induction intensity and the unit k mass center relative
Figure BDA00033107399400000313
The relative speed of the unit k replaces the corresponding quantity of the rest points of the unit k, the calculation complexity of the environment external force is reduced, and the numerical simulation efficiency is improved.
High-order earth gravitational acceleration in
Figure BDA00033107399400000314
Is gu,k＝[gθ,k gλ,k gr,k]TThen the gravity acting on the unit k is at
Figure BDA00033107399400000315
Component (b):
Figure BDA0003310739940000031
wherein
Figure BDA0003310739940000032
Is composed of
Figure BDA00033107399400000316
To
Figure BDA00033107399400000317
The rotation matrix of (2).
The Lorentz force acting on each unit is obtained according to a calculation formula of the Lorentz force
Figure BDA00033107399400000318
Component (b):
Figure BDA0003310739940000033
wherein
Figure BDA0003310739940000034
Is composed of
Figure BDA00033107399400000319
To
Figure BDA00033107399400000320
Rotation matrix of, ItTo insulate the current on the conductive tether, Bk,cIs the magnetic induction at the centroid of the cell k.
According to the calculation formula of the atmospheric resistance, the atmospheric resistance acting on the unit k is
Figure BDA00033107399400000321
The components in (a) are:
Figure BDA0003310739940000035
where ρ isaIs the atmospheric density at the centroid of cell k, CdIs a coefficient of resistance, SkThe area of the wind-facing surface is,
Figure BDA0003310739940000036
is composed of
Figure BDA00033107399400000322
To
Figure BDA00033107399400000323
Of the rotation matrix，vS_I′Is the speed of movement of the center of mass of cell k relative to atmosphere.
The integration of (2), (3), (4) allows all environmental external forces acting on unit k:
fk e＝fk g+fk B+fk d (5)
from the relationship between force and moment, all the ambient external moments acting on unit k are derived:
Figure BDA0003310739940000037
wherein
Figure BDA0003310739940000038
Is composed of
Figure BDA00033107399400000324
To
Figure BDA00033107399400000325
The kinetic equation of the unit k is rewritten into a matrix form by combining the formulas (1), (5) and (6):
Figure BDA0003310739940000039
wherein:
Figure BDA00033107399400000310
a quality matrix for cell k;
Figure BDA0003310739940000041
is the acceleration vector of unit k;
Figure BDA0003310739940000042
is a non-linear term of the cell;
Figure BDA0003310739940000043
is the restraining force vector of the hinge k to the unit k;
Figure BDA0003310739940000044
is the restraining force vector of hinge k +1 to unit k.
Step three: based on the dynamic equation of the single unit obtained by derivation in the second step, a discrete electrodynamic force rope model considering the coupling effect of multiple physical fields is established by combining a Kane equation and a multi-body dynamic recursive algorithm, and the characteristic that the calculation complexity of the recursive algorithm is in linear relation with the freedom degree of an electrodynamic force rope system is utilized, so that on the premise of ensuring the unchanged precision, the increase of discrete units in the discrete electrodynamic force rope model is reduced as much as possible, the high-precision and high-efficiency dynamic modeling of the electrodynamic force rope is realized, and the numerical simulation calculation amount and time consumption are reduced.
The dynamic equation of the electrodynamic force rope system is deduced by a multi-body dynamics recursion algorithm and mainly comprises the following two steps: unit kinematics recursion and unit dynamics recursion. Firstly, performing kinematic recursion of a unit k, wherein the velocity and acceleration recursion relations of two adjacent units k and k-1 of the electrodynamic force rope system are respectively as follows:
Figure BDA0003310739940000045
Figure BDA0003310739940000046
wherein:
Figure BDA0003310739940000047
is the velocity vector of unit k, is the translation velocity, ωkIs the rotational speed;
Figure BDA0003310739940000048
for the transfer matrix of cell k-1 to cell k, I3×3Is a 3 × 3 unit matrix, 03×3Is a 3 x 3 zero matrix and is characterized in that,
Figure BDA0003310739940000049
is composed of
Figure BDA00033107399400000412
To
Figure BDA00033107399400000413
The rotation matrix of (a);
Ub,kprojecting a matrix for hinge k;
Figure BDA00033107399400000410
in order to rotate the non-linear term,
Figure BDA00033107399400000411
for the angular velocity of unit k relative to unit k-1 at
Figure BDA00033107399400000414
Of (c) is used.
And (4) obtaining the speed and acceleration recurrence relation of all the units from the main satellite to the sub-satellites by using the speed recurrence relation (8) and the acceleration recurrence relation (9).
After the kinematics recursion of the electrodynamic force tether system is completed, the dynamics recursion of the unit k-1 is carried out, and the inertia matrix and nonlinear term recursion relations of two adjacent units k and k-1 of the conductive tether are respectively as follows:
Figure BDA0003310739940000051
Figure BDA0003310739940000052
wherein: mk-1Is the inertia matrix when cell k-1 is not updated,
Figure BDA0003310739940000053
updating the inertia matrix for the unit k-1;
Figure BDA0003310739940000054
Figure BDA0003310739940000055
an updated inertial matrix for unit k, "T" represents the matrix transpose;
Figure BDA0003310739940000056
Figure BDA0003310739940000057
is composed of
Figure BDA0003310739940000058
The inverse matrix of (d);
Figure BDA0003310739940000059
for the non-linear term when cell k-1 is not updated,
Figure BDA00033107399400000510
the updated nonlinear term for the unit k-1;
Figure BDA00033107399400000511
updated non-linear terms for unit k
The inertia matrix and nonlinear terms of all units from the subsatellite to the main star are derived by using the equations (10) and (11), and the kinetic equation (7) of the unit k is rewritten as:
Figure BDA00033107399400000512
and specifies that:
Figure BDA00033107399400000513
when k is equal to 0, the acceleration of the main star is obtained
Figure BDA00033107399400000514
Expression (c):
Figure BDA00033107399400000515
the angular acceleration of the unit k relative to the unit k-1 is known according to a recursion algorithm derivation process
Figure BDA00033107399400000516
The relative angular acceleration between all the units is derived by combining the equations (9), (10), (11), (13), (14):
Figure BDA00033107399400000517
the obtained formula (15) is an electrodynamic force rope system dynamics model, the electrodynamic force rope system dynamics model utilizes the characteristic that the calculation complexity of a recursion algorithm is in a linear relation with the system degree of freedom, and on the premise that the accuracy is not changed, the problem that numerical simulation calculation amount is increased due to the fact that discrete units in the discrete electrodynamic force rope model are increased is reduced as much as possible, namely the high-accuracy and high-efficiency dynamics modeling of the electrodynamic force rope is achieved.
Step four: and (4) applying the dynamic model of the electrodynamic force rope system established in the step three to the field of dynamics and control of the orbit and the attitude of the spacecraft, and solving the technical problem of the related engineering of the orbit and the attitude of the spacecraft.
The dynamics and control field of the spacecraft orbit and the attitude comprises the stability analysis of the attitude of an electric power rope system, the attitude angle control of the electric power rope system and the orbit transfer of the electric power rope system.
Has the advantages that:
1. the invention discloses a high-precision high-efficiency dynamic modeling method for an electrodynamic force rope, which is characterized in that a discrete electrodynamic force rope model considering the coupling effect of multiple physical fields is established by combining a Kane equation and a multi-body dynamic recursive algorithm, the characteristic that the computational complexity of the recursive algorithm is in linear relation with the system freedom degree is utilized, the problem of numerical simulation calculation amount increase caused by the increase of discrete units in the discrete electrodynamic force rope model is reduced as much as possible on the premise of ensuring the unchanged precision, and the high-precision high-efficiency dynamic modeling of the electrodynamic force rope is realized. The high-precision and high-efficiency dynamic modeling method of the electrodynamic force rope is applied to the field of dynamics and control of orbits and postures of spacecrafts, and the technical problems of relevant engineering are solved.
2. The invention discloses a high-precision high-efficiency dynamic modeling method for an electrodynamic force rope, which considers the influence of a high-order gravitational field, a geomagnetic field, atmospheric density and plasma density on the dynamic characteristics of a conductive rope through formulas (2), (3) and (4), and has stronger universality.
3. The invention discloses a high-precision high-efficiency dynamic modeling method for an electrodynamic force rope, which uses the gravity acceleration of a unit k mass center, the geomagnetic field magnetic induction intensity and the unit k mass center relative to each other when calculating the environmental external force
Figure BDA0003310739940000061
The relative speed replaces the corresponding quantity of the rest points of the unit k, the calculation of the environment external force through integration is avoided, the calculation complexity of the environment external force is reduced, the numerical simulation calculation quantity and time consumption are reduced, and the numerical simulation prediction efficiency is improved.
Drawings
FIG. 1 is a schematic diagram of an electric power roping system operating in track;
FIG. 2 depicts the attitude between adjacent cells of the conductive tether;
FIG. 3 direct method and recursive algorithm computation quantity comparison;
FIG. 4 shows simulation durations required by different units of a recurrence algorithm;
FIG. 5 is a schematic diagram of a discrete volume virtual bar model;
FIG. 6 comparison of virtual pole attitude angles at different voltages, where: FIG. 6(a) is VsInner angle at 5000V, fig. 6 (b)) Is a Vs8000V, V in FIG. 6(c)sOut-of-plane angle at 5000V, V in FIG. 6(d)s8000V out-of-plane angle;
FIG. 7 comparison of semi-major axis and eccentricity at different voltages, where: FIG. 7(a) shows VsSemi-major axis at 5000V, and V in FIG. 7(b)s8000V, and V in FIG. 7(c)sEccentricity of 5000V, V in FIG. 7(d)sEccentricity 8000V;
fig. 8 comparison of conductive tether currents at different voltages, where: FIG. 8(a) shows VsCurrent at 5000V, V in fig. 8(b)s8000V.
FIG. 9 is a flow chart of a high-precision and high-efficiency dynamic modeling method for an electrodynamic force rope disclosed by the invention.
Detailed Description
To better explain the objects and advantages of the present invention, the following detailed description of the embodiments and effects of the present invention will be made with reference to the examples and the accompanying drawings.
As shown in fig. 9, the high-precision and high-efficiency dynamics modeling method for the electrodynamic force rope disclosed in this embodiment includes the following specific implementation steps:
the method comprises the following steps: numbering each unit and hinge of the system;
the motion of the electrodynamic force rope system on the track is shown in figure 1, a main star and a sub-star in the electrodynamic force rope system are taken as mass points, the conductive rope is scattered into a plurality of rigid rod units, the main star is recorded as a unit 0, the conductive rope unit close to the main star is recorded as a unit 1, then the number of the rigid rod units is increased gradually to N, the sub-star is recorded as a unit N +1, all the units are connected through a spherical hinge, and as shown in figure 2, the relative attitude angle between all the units is taken as a generalized coordinate.
Step two: deducing a kinetic equation of a single unit by using a Kane equation;
the general form of the Kane equation is:
Fk,i I+Fk,i A＝0 (17)
wherein Fk,i IGeneralized inertial force representing ith generalized velocity of k-th order of object，Fk,i AThe generalized main power of the ith-order generalized velocity of the object k can be represented, and a point S on the conductive tether unit k can be obtained according to the algorithm of the vector array
Figure BDA00033107399400000712
Vector, velocity and acceleration components:
Figure BDA0003310739940000071
Figure BDA0003310739940000072
Figure BDA0003310739940000073
wherein r iskIs a hinge k at
Figure BDA00033107399400000713
Vector of middle, lSIs dotted at
Figure BDA00033107399400000714
The vector of (1) is selected as vkAnd ωkAs the generalized rate of unit k, then vSRelative vkAnd ωkYaw rate of
Figure BDA0003310739940000074
And
Figure BDA0003310739940000075
are respectively as
Figure BDA0003310739940000076
According to the definition of the generalized inertia force and the generalized main force of the Kane equation, the relative v of the unit k can be obtainedkHas a generalized inertial force of
Figure BDA0003310739940000077
Unit k vskHas a generalized inertial force of
Figure BDA0003310739940000078
Unit k vs vkHas a generalized main power of
Figure BDA0003310739940000079
Unit k vskHas a generalized main power of
Figure BDA00033107399400000710
Wherein
Figure BDA00033107399400000711
Substituting equation (20) into equation (17) for all the resultant forces acting on unit k can obtain the translational kinetic equation and the rotational kinetic equation of unit k, and the final result is equation (1).
The environmental external force borne by the electrodynamic force rope system in the low-orbit is mainly the earth gravitation, the Lorentz force and the atmospheric resistance.
For the gravity of the earth, the description is generally carried out by adopting a spherical harmonic series. The expression of the gravitational potential energy function of a point S on a unit k in the NEU coordinate system is as follows:
Figure BDA0003310739940000081
wherein G is a constant of universal gravitation, and θ is a point S
Figure BDA00033107399400000813
The latitude of middle (middle) is,λ is longitude, M is earth mass, R is distance of O to point S, R iseIs the radius of the earth, Cm,nAnd Sm,nIs a spherical harmonic series of pm,n(x) Is an m-order n-th order associated Legendre function, the gravity acceleration is
Figure BDA00033107399400000814
Is defined as follows:
Figure BDA0003310739940000082
wherein
Figure BDA0003310739940000083
In order to be the Nabla operator, the operation,
Figure BDA0003310739940000084
are respectively as
Figure BDA00033107399400000815
The three axes of rotation of the shaft(s),
Figure BDA0003310739940000085
specific expressions are as follows
Figure BDA0003310739940000086
Figure BDA0003310739940000087
Figure BDA0003310739940000088
The vector array is written in the form of
gs＝[gθ gλ gr]T
The gravity acting on the unit k can be obtained by substituting r at the centroid of the unit k into the above expression and then substituting the above expression into the expression (2).
The lorentz force of the unit k needs to calculate the geomagnetic field magnetic induction B at the centroid of the unit k and the current on the conductive tether respectively. The magnetic induction intensity is generally described by adopting a spherical harmonic series, and the magnetic potential energy of a point S on a unit k is
Figure BDA00033107399400000816
The expression in (1) is:
Figure BDA0003310739940000089
wherein
Figure BDA00033107399400000810
And
Figure BDA00033107399400000811
is a coefficient of a gaussian function, and is,
Figure BDA00033107399400000812
the latitude is the remaining latitude. The magnetic induction intensity of the geomagnetic field is in
Figure BDA00033107399400000817
Is defined as follows:
Figure BDA0003310739940000091
it is unfolded with
Figure BDA0003310739940000092
Figure BDA0003310739940000093
Figure BDA0003310739940000094
It is marked as
Bs＝[Bθ′ Bλ Br]T (24-b)
The electric field generated by the unit k centroid cutting magnetic induction line
Figure BDA0003310739940000095
In that
Figure BDA00033107399400000912
Component (A) of
Figure BDA0003310739940000096
Wherein
Figure BDA0003310739940000097
Is composed of
Figure BDA00033107399400000913
To
Figure BDA00033107399400000914
V.ofk,cIs a unit k with a center of mass of
Figure BDA00033107399400000915
Velocity component of (1), ωI'Of rotational angular velocity of the earth
Figure BDA00033107399400000916
Component rk,cThe vector from geocentric to k centroid of cell is
Figure BDA00033107399400000917
The induced electromotive force on the cell k can be obtained from the component of (b):
Figure BDA0003310739940000098
when an insulating cord is used as the conductive tether, the current on the conductive tether is:
Figure BDA0003310739940000099
wherein VsIs the supply voltage, VEMFInducing an electromotive force, R, to the entire conductive tethertIs a conductive tether resistance, e is an amount of electron charge, RacIn order to be the radius of the charge-exchange device,
Figure BDA00033107399400000910
is the average earth magnetic field at the centroid of each cell, meFor electron mass, AacIs the area of the charge exchange device, NeIs the electron density, kBIs the Bolmatz constant, TeIs the electron temperature, σ1And σ2The formula is a nonlinear equation which is an empirical coefficient, the conductive tether current can be obtained by solving through a Newton iteration method, and finally Lorentz force acting on the unit k can be obtained by substituting the formulas (24) and (27) into the formula (3).
According to
Figure BDA00033107399400000918
And
Figure BDA00033107399400000919
the velocity of a point S on a cell k relative to the atmosphere can be expressed as:
Figure BDA00033107399400000911
the atmospheric resistance of the action unit k can be obtained by substituting the above formula into (4).
Step three: and (4) further deducing the kinetic equation of the whole electrodynamic force rope system by utilizing a multi-body dynamics recursive algorithm based on the single rigid body kinetic equation obtained in the step one.
As shown in FIG. 2, the bit-vector relationship between the unit k and the unit k-1 is:
Figure BDA0003310739940000101
wherein
Figure BDA0003310739940000102
The bit vector from point O to hinge k,
Figure BDA0003310739940000103
the bit vector from point O to hinge k-1,
Figure BDA0003310739940000104
for the position vector from the hinge k to the hinge k-1, the first derivative and the second derivative of the above equation are respectively obtained with respect to time to obtain the relationship between the speed and the acceleration of the two units:
Figure BDA0003310739940000105
Figure BDA0003310739940000106
the kinematic recurrence relation of the expressions (8) and (9) can be obtained by rewriting the above two expressions into a vector array form, respectively.
After obtaining the kinematic recurrence relation of all units, the kinetic recurrence relation among the units is deduced as follows, and the pair star has the following dynamic equation of any unit given by the formula (7):
Figure BDA0003310739940000107
since the subsatellite has no constraint of the hinge N +2, FN+2,cN+1When it is equal to 0, then order again
Figure BDA0003310739940000108
The above formula is rewritten as:
Figure BDA0003310739940000109
f can be deduced according to the acceleration recursion relation of the units N +1 and NN+1,cN+1Expressed as the acceleration of unit N:
Figure BDA00033107399400001010
because the binding force of hinge N +1 pair child star and electrically conductive tether unit N are the reaction each other, consequently have:
Figure BDA00033107399400001011
wherein
Figure BDA00033107399400001012
For transforming the matrix of restraining forces, the restraining forces are transformed from a coordinate system
Figure BDA00033107399400001016
Conversion to a coordinate system
Figure BDA00033107399400001017
The kinetic equation of the unit N can be dynamically updated by using the formulas (30) and (31):
Figure BDA00033107399400001013
it is simplified to obtain:
Figure BDA00033107399400001014
wherein:
Figure BDA00033107399400001015
Figure BDA0003310739940000111
for the remaining units, the process of equations (29) - (30) can be repeated to perform dynamic recursion on each unit in turn until the main star:
Figure BDA0003310739940000112
since the main star is only constrained by the hinge 1, F0,c0Finally, combining the equations (9), (10), (11), (13) and (14) to derive the relative angular acceleration between all units
Figure BDA0003310739940000113
Step four: and analyzing the attitude and track dynamics characteristics of the electrodynamic force rope system under different initial conditions by using the dynamic model of the electrodynamic force rope system obtained by deduction in the first step and the second step, and analyzing the real calculation efficiency of a recursion algorithm.
Given system parameters and initial conditions of the orbit as shown in tables 1 and 2, the orbit eccentricity of the main satellite is 0.01, the initial attitude angle and the angular velocity of each unit are 0, the conductive tether is dispersed into 3 units, and a virtual rod of the electric power tether system is defined as a connecting line between the sub satellite and the main satellite as shown in fig. 5.
TABLE 1 electrodynamic roping system parameters
Figure BDA0003310739940000114
TABLE 2 Primary Star initial orbit conditions
Figure BDA0003310739940000115
When supplying a power supply voltage VsAt 5000V, the attitude angle of the virtual pole of the electrodynamic force rope system changes as shown in fig. 6(a) and 6(c), and the virtual pole is in an uncontrolled stateThe amplitude of the in-plane angular oscillation of the track is gradually increased, the amplitude of the out-of-plane angular oscillation is periodically changed, and the Lorentz force continuously applies negative work to the system, so that the semi-major axis and the eccentricity of the track are in a descending trend, as shown in FIG. 7(a) and FIG. 7 (c); when the supply voltage is increased to 8000V, the tether current increases, as shown in FIG. 8, which in turn causes an increase in Lorentz force, and hence increases in both the magnitude of the imaginary inside and outside angles, while the orbit semi-major axis and eccentricity are both greater than VsThere was also a drop at 5000V, indicating that the electrodynamic roping system is more efficient off-track at the same time.
Fig. 3 shows the calculation complexity of the direct solver and recursion algorithms as a function of the number of discrete elements, with the recursion algorithm being less complex when the number of discrete elements of the system exceeds 15. Fig. 4 shows the real simulation duration required by different discrete unit numbers, and the simulation duration curve in fig. 4 shows an approximately linear increase, which is close to the theoretical derivation result in fig. 3, and illustrates that the dynamic model of the electrodynamic force rope system derived based on the recursion algorithm has better calculation efficiency when the unit freedom is more.
The above detailed description is intended to illustrate the objects, aspects and advantages of the present invention, and it should be understood that the above detailed description is only exemplary of the present invention and is not intended to limit the scope of the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (4)

1. A high-precision and high-efficiency dynamic modeling method for an electrodynamic force rope is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
the method comprises the following steps: in an electrodynamic force tether system, a main satellite and a sub-satellite are taken as mass points, a conductive tether is discretized into a plurality of rigid rod units, namely a discretization model hinged by a plurality of rigid rods is adopted, the main satellite is recorded as a unit 0, a conductive tether unit close to the main satellite is recorded as a unit 1, then the number of the conductive tether units is sequentially increased to N, the sub-satellite is recorded as a unit N +1, the units are connected through a spherical hinge, and a relative attitude angle between the units is selected as a generalized coordinate;
step two: deducing a kinetic equation of a single unit by using a Kane equation, wherein the single unit refers to the main star, the conductive tether unit and the subsatellite defined in the step one;
the second step is realized by the method that,
in modeling of dynamics of an electrodynamic tether system, coordinate systems OXYZ, OX ' Y ' Z ', OX are definedoyozo，oxuyuzu，oxdydzd，okxkykzkSequentially representing the inertial system FIEarth-based, orbital coordinate system FoNorth-East-Up NEU coordinate system F for shortuNorth-East-Down is called NED coordinate system F for shortdUnit k body coordinate system Fk(ii) a Wherein the system of inertia FIThe origin O is located at the center of mass of the earth, the OX axis points to the spring equinox, the OZ axis is consistent with the rotation axis of the earth, and the OY axis and the other two axes form a right-hand system; earth following coordinate system FI′The origin is O, the OZ 'axis is consistent with the OZ axis, and the OX' axis and the OY 'axis rotate around the OZ' axis at the same time at the earth rotation angular speed; orbital coordinate system FoThe origin o is located at the centroid of the main star, oxoThe axis pointing along the earth's centre O to O, ozoThe axis is vertical to the orbital plane and is consistent with the orbital angular momentum direction of the main satellite, oyoThe shaft and the other two shafts form a right-hand system; NEU coordinate system FuThe origin is o, oxuThe axis pointing along point O, oyuThe axis is vertical to the meridian plane of the local area and points to the east, and the oz' axis and the other two axes form a right-handed system; NED coordinate system FdThe origin is o, ozdThe axis pointing along point O, oydThe axis is perpendicular to the meridian plane of the local region and points to east, oxdThe shaft and the other two shafts form a right-hand system; conductive tether unit k body coordinate system FkOrigin okWhen the whole conductive tether is consistent with the local vertical line of the main star, the body coordinate system of the unit k is superposed with the track system; selecting an attitude angle theta between the unitsk,k-1＝[θk,k-1z θk,k-1y θk,k-1x]TAs a generalized variable:
q＝[θ1,o θ2,1 … θN,N-1]T
and defines the attitude angle of the unit k with respect to the unit k-1 as thetak,k-1zAnd thetak,k-1yWherein thetak,k-1zIs xkAxis is at okxkykProjection of plane and xk-1Angle of axis thetak,k-1yIs xkAxis is at okxkykProjection of plane and xkAngle of axis thetak,k-1xFor the last attitude angle and specifying Fk-1To FkThe rotation sequence is 3-2-1, and the subsatellite is fixedly connected with the conductive tether unit N;
based on the algorithm of the vector array, dynamic modeling is carried out on the single unit by adopting a Kane equation, and the translation and rotation dynamic equations of the conductive tether unit k are obtained through derivation and arrangement:
Figure FDA0003577105820000021
in the formula mkMass of electrically conductive tether element k, JkIs unit k relative to okAt a moment of inertia of FkDescription of (1), SkIs unit k relative to okStatic moment ofk"x" represents a diagonally symmetric matrix of vectors, akIs unit k at FIMedium translational acceleration, ωkIs unit k relative to FIAngular velocity of FkThe component (b) of (a) is,
Figure FDA0003577105820000022
in order to be able to take into account the angular acceleration,
Figure FDA0003577105820000023
is FkTo FIOf the rotation matrix fk e、fk,ck、fk+1,ckRespectively, an environmental external force acting on the unit k, a constraint force of the hinge k to the unit k, and a constraint force of the hinge k +1 to the unit k, which are all represented at FIIn, Tk eRepresenting the ambient external moment, T, acting on unit kk+1,ckDisplay hingeThe constraint moments of chain k +1 on unit k, said moments all being represented at FkPerforming the following steps;
the external force acting on the electrodynamic force rope system in the low-orbit environment mainly comprises the earth attraction, the Lorentz force and the atmospheric resistance; using the gravity acceleration at the unit k centroid, the geomagnetic field magnetic induction and the unit k centroid relative FI′The relative speed of the unit k replaces the corresponding quantity of the rest points of the unit k, the calculation complexity of the environment external force is reduced, and the numerical simulation efficiency is improved;
high-order earth gravitational acceleration in FuThree components of (A) are
Figure FDA00035771058200000211
The gravitational force acting on unit k is at FIComponent (b):
Figure FDA0003577105820000024
wherein
Figure FDA0003577105820000025
Is Fu,kTo FIThe rotation matrix of (a);
the Lorentz force acting on each unit is obtained at F according to a calculation formula of the Lorentz forceIComponent (b):
Figure FDA0003577105820000026
wherein
Figure FDA0003577105820000027
Is Fd,kTo FkRotation matrix of, ItTo insulate the current on the conductive tether, Bk,cIs the magnetic induction at the centroid of cell k;
according to the calculation formula of the atmospheric resistance, the atmospheric resistance acting on the unit k is FIThe components in (a) are:
Figure FDA0003577105820000028
where ρ isaIs the atmospheric density at the centroid of cell k, CdIs a coefficient of resistance, SkThe area of the wind-facing surface is,
Figure FDA0003577105820000029
is FITo FI′V.ofS_I′The moving speed of the mass center of the unit k relative to the atmosphere is obtained;
the integration of (2), (3), (4) allows all environmental external forces acting on unit k:
fk e＝fk g+fk B+fk d (5)
from the relationship between force and moment, all the ambient external moments acting on unit k are derived:
Figure FDA00035771058200000210
wherein
Figure FDA0003577105820000031
Is FITo FkThe kinetic equation of the unit k is rewritten into a matrix form by combining the formulas (1), (5) and (6):
Figure FDA0003577105820000032
wherein:
Figure FDA0003577105820000033
a quality matrix for cell k;
Figure FDA0003577105820000034
is the acceleration vector of unit k;
Figure FDA0003577105820000035
is a non-linear term of the cell;
Figure FDA0003577105820000036
is the restraining force vector of the hinge k to the unit k;
Figure FDA0003577105820000037
is the constraint force vector of the hinge k +1 to the unit k;
step three: based on the dynamic equation of the single unit obtained by derivation in the second step, a discrete electrodynamic force rope model considering the coupling effect of multiple physical fields is established by combining a Kane equation and a multi-body dynamic recursive algorithm, and the characteristic that the calculation complexity of the recursive algorithm is in linear relation with the freedom degree of an electrodynamic force rope system is utilized, so that on the premise of ensuring the unchanged precision, the increase of discrete units in the discrete electrodynamic force rope model is reduced as much as possible, the high-precision and high-efficiency dynamic modeling of the electrodynamic force rope is realized, and the numerical simulation calculation amount and time consumption are reduced.
2. A high accuracy and high efficiency kinetic modeling method of electrodynamic ropes according to claim 1, characterized by: and step four, applying the dynamic model of the electrodynamic force rope system established in the step three to the field of dynamics and control of the orbit and the attitude of the spacecraft, and solving the technical problem of the related engineering of the orbit and the attitude of the spacecraft.
3. A high accuracy and high efficiency kinetic modeling method of electrodynamic ropes according to claim 2, characterized in that: the dynamics and control field of the spacecraft orbit and the attitude comprises the stability analysis of the attitude of an electric power rope system, the attitude angle control of the electric power rope system and the orbit transfer of the electric power rope system.
4. A high accuracy and high efficiency kinetic modeling method of electrodynamic ropes as claimed in claim 1, 2 or 3, characterized by: the third step is to realize the method as follows,
the dynamic equation of the electrodynamic force rope system is deduced by a multi-body dynamics recursion algorithm and mainly comprises the following two steps: a unit kinematics recursion and a unit dynamics recursion; firstly, performing kinematic recursion of a unit k, wherein the velocity and acceleration recursion relations of two adjacent units k and k-1 of the electrodynamic force rope system are respectively as follows:
Figure FDA0003577105820000038
Figure FDA0003577105820000039
wherein:
Figure FDA0003577105820000041
is the velocity vector of unit k, is the translation velocity, ωkIs the rotational speed;
Figure FDA0003577105820000042
for the transfer matrix of cell k-1 to cell k, I3×3Is a 3 × 3 unit matrix, 03×3Is a 3 x 3 zero matrix and is,
Figure FDA0003577105820000043
is Fk-1To FkThe rotation matrix of (a);
Ub,kprojecting a matrix for hinge k;
Figure FDA0003577105820000044
in order to rotate the non-linear term,
Figure FDA0003577105820000045
for the angular velocity of unit k relative to unit k-1 at FkThe component (b);
obtaining the speed and acceleration recurrence relation of all units from the main satellite to the sub satellite by using a speed recurrence relation (8) and an acceleration recurrence relation (9);
after the kinematics recursion of the electrodynamic force tether system is completed, the dynamics recursion of the unit k-1 is carried out, and the inertia matrix and nonlinear term recursion relations of two adjacent units k and k-1 of the conductive tether are respectively as follows:
Figure FDA0003577105820000046
Figure FDA0003577105820000047
wherein: mk-1Is the inertia matrix when cell k-1 is not updated,
Figure FDA0003577105820000048
updating the inertia matrix for the unit k-1;
Figure FDA0003577105820000049
Figure FDA00035771058200000410
an updated inertial matrix for unit k, "T" represents the matrix transpose;
Figure FDA00035771058200000411
Figure FDA00035771058200000412
is composed of
Figure FDA00035771058200000413
The inverse matrix of (d);
Figure FDA00035771058200000414
for the non-linear term when the cell k-1 is not updated,
Figure FDA00035771058200000415
the updated nonlinear term for the unit k-1;
Figure FDA00035771058200000416
updated non-linear terms for unit k;
the inertia matrix and nonlinear terms of all units from the subsatellite to the main star are derived by using the equations (10) and (11), and the kinetic equation (7) of the unit k is rewritten as:
Figure FDA00035771058200000417
and specifies that:
Figure FDA00035771058200000418
when k is 0, the acceleration of the main star is obtained
Figure FDA00035771058200000419
Expression (c):
Figure FDA00035771058200000420
the angular acceleration of the unit k relative to the unit k-1 is known according to a recursion algorithm derivation process
Figure FDA0003577105820000051
The relative angular acceleration between all the units is derived by combining the equations (9), (10), (11), (13), (14):
Figure FDA0003577105820000052
the obtained formula (15) is an electrodynamic force rope system dynamics model, the electrodynamic force rope system dynamics model utilizes the characteristic that the calculation complexity of a recursion algorithm is in a linear relation with the system degree of freedom, and on the premise that the accuracy is not changed, the problem that numerical simulation calculation amount is increased due to the fact that discrete units in the discrete electrodynamic force rope model are increased is reduced as much as possible, namely the high-accuracy and high-efficiency dynamics modeling of the electrodynamic force rope is achieved.
CN202111224605.5A 2021-10-19 2021-10-19 Efficient dynamics modeling method for electrodynamic force rope derailing device Active CN113935176B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111224605.5A CN113935176B (en) 2021-10-19 2021-10-19 Efficient dynamics modeling method for electrodynamic force rope derailing device

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111224605.5A CN113935176B (en) 2021-10-19 2021-10-19 Efficient dynamics modeling method for electrodynamic force rope derailing device

Publications (2)

Publication Number Publication Date
CN113935176A CN113935176A (en) 2022-01-14
CN113935176B true CN113935176B (en) 2022-05-10

Family

ID=79281016

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111224605.5A Active CN113935176B (en) 2021-10-19 2021-10-19 Efficient dynamics modeling method for electrodynamic force rope derailing device

Country Status (1)

Country Link
CN (1) CN113935176B (en)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007681A (en) * 2018-11-28 2019-07-12 北京理工大学 It is a kind of to realize that rope is formation spinning stability expansion optimization method using continuous propeller

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107122515B (en) * 2017-03-17 2020-06-23 北京航空航天大学 Dynamics analysis method of rope system transportation system based on absolute node coordinate method
CN107908855B (en) * 2017-11-10 2021-02-05 南京航空航天大学 Modeling method of rigid-flexible coupling space belt-shaped rope system
CN109002050B (en) * 2018-07-02 2021-04-27 南京航空航天大学 Modeling method for space three-body flexible tether satellite formation system under non-inertial reference system
CN110210047B (en) * 2019-03-19 2021-02-05 南京航空航天大学 Method for constructing release dynamics model of strap-shaped tethered satellite

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110007681A (en) * 2018-11-28 2019-07-12 北京理工大学 It is a kind of to realize that rope is formation spinning stability expansion optimization method using continuous propeller

Also Published As

Publication number Publication date
CN113935176A (en) 2022-01-14

Similar Documents

Publication Publication Date Title
Slavinskis et al. High spin rate magnetic controller for nanosatellites
Zhong et al. Libration dynamics and stability of electrodynamic tethers in satellite deorbit
Zappulla et al. Experiments on autonomous spacecraft rendezvous and docking using an adaptive artificial potential field approach
Liu et al. Collision-free trajectory design for long-distance hopping transfer on asteroid surface using convex optimization
Huang et al. Nonlinear dynamics and reconfiguration control of two-satellite Coulomb tether formation at libration points
Tsai et al. LQR motion control of a ball-riding robot
Oland et al. Reaction wheel design for cubesats
CN113935176B (en) Efficient dynamics modeling method for electrodynamic force rope derailing device
Ousaloo et al. Verification of Spin Magnetic Attitude Control System using air-bearing-based attitude control simulator
Biggs et al. Motion planning on a class of 6-D Lie groups via a covering map
Osterloh et al. A rigid body dynamics simulation framework for the analysis of attitude control systems of modular satellite systems
Shibli Unified modeling approach of kinematics, dynamics and control of a free-flying space robot interacting with a target satellite
Panosian et al. Tethered Coulomb structure applied to close proximity situational awareness
Quadrelli et al. Active electrostatic flight for airless bodies
CN112389679B (en) Moonlet constant thrust orbit recursion method considering multiple perturbation forces
Sakamoto et al. Verification of tether deployment system aboard CubeSat through dynamics simulations and tests
Simo et al. Solar sail trajectories at the Earth-Moon Lagrange points
Joe et al. Formation dynamics of coulomb satellites
Weis et al. Attitude control for chip satellites using multiple electrodynamic tethers
Zorita Dynamics of small satellites with gravity gradient attitude control
Kalita et al. Mobility and Science Operations on an Asteroid using a Hopping Small Spacecraft on Stilts
Saaj et al. Electrostatic forces for satellite swarm navigation and reconfiguration
Zhang et al. A high-fidelity high-efficiency model for electrodynamic tether system based on recursive algorithm
Sun et al. Attitude tracking control of gravity gradient microsatellite in maneuvering for space exploration
Feng et al. Dynamics and momentum equalization control of redundant space robot with control moment gyroscopes for joint actuation

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