Double asteroid system detector orbit design method based on finite element method

Abstract

The invention discloses a double asteroid system detector orbit design method based on a finite element method, and belongs to the technical field of spacecraft space detection. Comprising the following steps that double asteroid data obtained through calculation are imported, loading data are input into a loading model, and simulation is conducted in combination with movement data of a double-satellite system; modeling each vertex in the asteroid by adopting a finite element method, and calculating the coordinate position of each vertex in an inertial system; and importing the calculated double-star motion data, setting input position parameters and speed parameters of the detector, and calculating to obtain an orbit of the detector under the finite element gravitational field model. According to the method, the interaction force of the irregular double asteroids can be accurately represented, and the orbit of the detector is planned according to the interaction force. By representing the positions of the spacecraft orbit and the double asteroids of the system in real time, the characteristics of the satellite orbit can be obtained, and the method can be applied to program simulation in the aspects of probe task planning and spacecraft orbit optimization. The method is applied to the field of spacecraft space simulation.

Description

Technical Field

[0001] The present invention relates to the field of spacecraft space exploration technology, and in particular to a method for designing a probe orbit for a double asteroid system based on a finite element method. Background Art

[0002] With the development and progress of aerospace technology, humans have set their sights on the distant deep space. Researchers' research on asteroids in the solar system is one of the hottest topics in deep space exploration. Because asteroids were formed at the same time as the solar system, but did not participate in the evolution of planets in the solar system, they retain a lot of early solar system materials and contain relatively complete early information about the formation of the sun, which can provide important clues to the origin and evolution of the early solar system. At the same time, some asteroids may contain particularly precious rare metals and other resources, which have great potential mining value. This puts higher demands on the detection and research of asteroids.

[0003] The unique formation mechanism of binary systems helps to better study the evolution of asteroids, and detecting two asteroids at a time can reduce launch costs and enrich mission results. Therefore, binary asteroid systems have gradually become a hot target for asteroid exploration missions in recent years. In 1993, the Galileo probe discovered the first double asteroid system - Ida243 and its satellite Dactyl. More and more double asteroid systems have been discovered one after another. The Lucy probe launched by the United States in October 2021 will flyby and detect multiple Trojan asteroids in the next 12 years, including the two double asteroid systems of Eurybates 3548 and Patroclus 617; in November 2021, the Double Asteroid Redirection Test launched by NASA visited the near-Earth double asteroid system Didymos 65803, and changed the orbit of the asteroid through a kinetic impactor, thereby verifying the defense technology of near-Earth asteroids; in addition, NASA's Small Innovative Planet Exploration (SIMPLEx) project has also passed the Janus mission, which will use the near-Earth double asteroid system FG31996 and VH With the goal of 1991, two lightweight, low-cost small spacecraft were used to perform high-resolution imaging and high-precision modeling to analyze the formation and evolution of binary systems. According to observational data, there are many binary asteroid systems in the universe, accounting for about 16% of the total number of near-Earth asteroids.

[0004] However, the close-in exploration activities of binary asteroid probes cannot be separated from the gravitational effect of the binary asteroids. The gravitational field near the binary asteroid system is closely related to the orientation and position of the two asteroids. Therefore, studying the attitude motion and orbital motion of the two asteroids is the key to the gravitational field of the binary asteroid system. The attitude-orbit motion of the asteroid is determined by the dynamics of the system. Considering that the two asteroids have more complex shapes, smaller volumes and masses, uneven mass distribution, and faster spin speeds, the gravitational field distribution near the two asteroids is irregular, weaker, and more drastic compared to general spherical large planets. As a result, the dynamics of the binary asteroid system often presents the characteristics of attitude-orbit coupling, that is, the full two-body problem. From the kinematics of the binary asteroid system, it can be seen that the gravitational mutual potential of the binary asteroid system is the basis for studying the full two-body problem. However, in general, there is no closed-form solution for the mutual potential between two asteroids with irregular shapes and non-uniform mass distributions.

[0005] In the past few decades, different scholars have been committed to developing methods that can calculate the mutual gravitational potential of binary star systems. This work is necessary for the preliminary study of the gravitational field of binary asteroids during the mission planning of the probe. The primary consideration of the research is to find an effective method to describe this gravitational field, which is also the basis for subsequent research. The commonly used gravitational field modeling methods currently include: simple geometric model, spherical harmonic and ellipsoidal harmonic series expansion model, particle group model, homogeneous polyhedron method modeling, etc. Among these modeling methods, each has its own advantages and disadvantages in terms of versatility, calculation speed, calculation accuracy, scope of application and additional storage. For example, Werner and Scheeres used infinite series to express the mutual gravitational potential between a pair of uniform polyhedrons; Fahnestock and Scheeres further studied this work and derived a method for calculating the derivative of the mutual gravitational potential between two objects with arbitrary positions and attitudes. Based on the above research, Scheeres et al. studied the dynamics of the near-Earth binary asteroid (66391) 1999KW4 and conducted dynamic simulation and analysis; Compere and Using the symmetric free trace (STF) tensor form developed by previous researchers, the description of the mutual potential between two bodies was elaborated and completed; Shi et al. proposed a shape-based mutual gravitational potential model and used it to simulate the evolution of the asteroid (66391) 1999KW4, in which one binary asteroid was modeled as a uniform polyhedron and the other was modeled as an extended rigid body with arbitrary mass distribution. However, the simple geometric model and spherical harmonic series expansion method are not suitable for asteroids with irregular shapes. The particle group model is more accurate, but because the asteroid often needs to be discretized when describing the gravitational field of the asteroid, this will result in limited spatial resolution of the model, and the mean polyhedron rule cannot characterize the mass distribution inside the asteroid. To solve these problems, Yu et al. proposed a finite element method for calculating the gravitational interaction between two arbitrarily shaped celestial bodies. This method has no restrictions on the shape and internal mass distribution of the asteroid, and has acceptable accuracy in describing the mutual gravitational potential energy between two non-strained gravitational bodies, thus avoiding the non-convergence problem caused by traditional methods. In order to improve the calculation accuracy of parameters such as potential energy, gravity and torque, a large number of grids are required, which places high demands on computing power. Therefore, the scale of finite element simulation of such celestial systems is usually very large, requiring high-performance computers to achieve acceleration.

[0006] The orbital design of asteroid exploration spacecraft also requires real-time calculation of the gravitational field of the binary asteroids. If the motion of the binary stars and the probe is integrated at every moment, the amount of calculation will be extremely large. However, the weak gravitational field of the asteroid determines that any disturbance to the probe will interfere with its orbit. To optimize the orbit, it is necessary to improve the calculation efficiency. Summary of the invention

[0007] The purpose of the present invention is to provide a method for designing a probe orbit for a binary asteroid system based on the finite element method. The program for designing a probe orbit for a binary asteroid system based on the finite element method uses finite element modeling to simulate the gravitational field of the target asteroid. Different types of orbits that can meet different mission requirements are designed based on this. This technology is suitable for spacecraft orbit planning and simulation, and meets the requirements of high simulation accuracy and low computational cost.

[0008] To achieve the above object, the present invention provides a method for designing a probe orbit for a double asteroid system based on a finite element method, comprising the following steps:

[0009] S1. Import the calculated binary asteroid data, input the loading data into the loading model, and simulate it in combination with the binary system motion data;

[0010] S2. Use the finite element method to perform finite element modeling on each vertex in the asteroid and calculate its coordinate position in the inertial system;

[0011] S3. Based on step S2, the calculated binary star motion data is imported, the input detector position parameters and velocity parameters are set, and the orbit of the detector under the finite element gravitational field model is obtained by calculation.

[0012] Preferably, in step S1, the specific process design of inputting the loading data into the loading model and simulating in combination with the binary star system motion data is as follows:

[0013] In the model, the tetrahedron c k From object A, and tetrahedron c m from object B. Let ρ A =(x a ,y a ,z a ) represents the element c k For any vertex position, let ρ B =(x b ,y b ,z b ) represents the element c m Any vertex position, where ρ A and ρ B Calculated from A and B respectively, the following equation is used to determine c k and c m The interaction function h(ρ A ,ρ B ) is the volume-based field integral:

[0014]

[0015] By implementing natural coordinates and using bilinear interpolation to evaluate the k ,c m The function h(ρ A ,ρ B ), allowing the analysis of integral equations within the defined domain, such as equation (2):

[0016]

[0017] Formula (2), J k and J m is the Jacobian, the shape function of the tetrahedron, k and m are the indices of the elements in objects A and B respectively, and i and j are the number of vertices of the two tetrahedrons under consideration;

[0018] Let A and B denote the set of all elements in rigid bodies A and B respectively;

[0019] The integral of the interaction function in the entire two-body domain is expressed as:

[0020]

[0021] Preferably, in step S2, finite element modeling is performed on each vertex in the asteroid using the finite element method, and the specific process is as follows:

[0022] Since each node belongs to multiple different elements at the same time, there are many redundant calculations in formula (3);

[0023] By considering the mesh nodes of two objects, the sum in equation (3) is simplified as follows:

[0024]

[0025] In formula (4), α and β are the node indices of the integrand, H αβ represents the value of the interaction function at nodes α and β, w α and w β is a custom node weight calculated as follows:

[0026]

[0027] In formula (5), α is the set of all elements in the grid point set A, and β is the set of all elements containing nodes in object B. The most time-consuming part of the analysis is the calculation of the gravitational potential and the calculation of the gravitational force on object A and the moments on objects A and B during the integration process. As the grid size increases, the computational workload increases significantly.

[0028] Preferably, in step S2, the specific calculation procedure and formula for calculating the coordinate position and interaction torque of the calculation model in the inertial system are as follows:

[0029] For any grid point α in object A and any grid point β in object B, the gravitational force f acting on α is obtained from the gravitational force exerted by β on α. αβ and the moment p acting on the center of mass of object A αβ , the moment p acting on the center of mass of object B is obtained from the gravitational force exerted by α on β αβ , as shown below:

[0030]

[0031] According to formula (4), the gravitational force on objects A and B and the moment on the two objects are estimated using the bilinear tetrahedron:

[0032]

[0033]

[0034] Preferably, in step S3, the calculated binary star motion data is imported according to the coordinate position of the model obtained in step S2 in the inertial system:

[0035] Introducing about w α and w β New definition of:

[0036] W α =w α σ(ρ Aα ) (12)

[0037] W β =w β σ(ρ Bβ ) (13)

[0038] Equations (9), (10) and (11) can be rewritten as:

[0039]

[0040] The calculation process for the gravitational potential of the binary star system, the gravitational force on A, the torque on A and the torque on B is as follows:

[0041] S31, input double precision data r Ai 、r Bj , Ai , Bj , W Ai and W Bj , where W Ai Yes Ai The product of the sum, W Bj Yes Bj The product of the sum;

[0042] S32, calculation formulas (14), (15) and (16);

[0043] S33, output result F AB 、T A and T B ;

[0044] In the calculation process of a single-core CPU, the process is as follows:

[0045] First, iteratively read the parameters associated with the nodes in objects A and B, and use these parameters to calculate

[0046] Calculate f for each iteration ij ,p ij ,q ij ,in:

[0047]

[0048] q ij =ρBj ×f ij (19)

[0049] Then, all calculated f ij ,、p ij and q ij The values ​​are added and multiplied by G to get F AB , T A and T B ; Among them, rigid body A has N1 nodes, rigid body B has N2 nodes, and it takes N1×N2 iterations to get N1×N2 f ij ,、p ij and q ij Numeric value;

[0050] The most time-consuming part of the program is to calculate the gravitational force and torque between the two binary asteroids, namely F AB , T A and T B ;

[0051] w in equations (12) and (13) α and w β The term is a constant that is computed once for each node during preprocessing;

[0052] Then, calculate equation (17), for p ij Calculation, only calculate ρ Ai and f ij The cross product of is as shown in formula (18).

[0053] Therefore, the present invention adopts a double asteroid system probe orbit design method based on the finite element method with the above structure, which has the following beneficial effects:

[0054] The finite element method used in the present invention can more accurately characterize the interaction force of the irregular binary asteroids, and use it to plan the orbit of the probe. By characterizing the spacecraft orbit and the position of the binary asteroids in the system in real time, the characteristics of the satellite orbit can be obtained, which can be applied to program simulation in probe mission planning and spacecraft orbit optimization. The present invention is used in the field of spacecraft space simulation.

[0055] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] Figure 1 A schematic diagram of the vertical orbit of the 1999KW4 double asteroid system probe according to a double asteroid system probe orbit design method based on the finite element method of the present invention;

[0057] Figure 2A schematic diagram of the horizontal orbit of the 1999KW4 double asteroid system probe according to a double asteroid system probe orbit design method based on the finite element method of the present invention;

[0058] Figure 3 This is a method for designing a probe orbit for a binary asteroid system based on the finite element method. k and c m The interaction function h(ρ A ,ρ B ) is a schematic diagram of the volume-based field integration. DETAILED DESCRIPTION

[0059] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0060] Unless otherwise defined, the technical terms or scientific terms used in the present invention should be understood by people with ordinary skills in the field to which the present invention belongs. The words "first", "second" and similar words used in the present invention do not indicate any order, quantity or importance, but are only used to distinguish different components. "Include" or "comprise" and similar words mean that the elements or objects appearing before the word include the elements or objects listed after the word and their equivalents, without excluding other elements or objects. "Connect" or "connected" and similar words are not limited to physical or mechanical connections, but may include electrical connections, whether direct or indirect. "Up", "down", "left", "right" and the like are only used to indicate relative positional relationships. When the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0061] Example

[0062] like Figure 1-2 As shown, the present invention provides a method for designing a probe orbit for a double asteroid system based on a finite element method, comprising the following steps:

[0063] S1. Import the calculated binary asteroid data, input the loading data into the loading model, and simulate it in combination with the binary system motion data;

[0064] S2. Use the finite element method to perform finite element modeling on each vertex in the asteroid and calculate its coordinate position in the inertial system;

[0065] S3. Based on step S2, the calculated binary star motion data is imported, the input detector position parameters and velocity parameters are set, and the orbit of the detector under the finite element gravitational field model is obtained by calculation.

[0066] In step S1, the specific process design of inputting the loading data into the loading model and combining the binary system motion data for simulation is as follows:

[0067] like Figure 3 As shown, in the model, the tetrahedron c k From object A, and tetrahedron c m from object B. Let ρ A =(x a ,y a ,z a ) represents the element c k For any vertex position, let ρ B =(x b ,y b ,z b ) represents the element c m Any vertex position, where ρ A and ρ B Calculated from A and B respectively. The following equation is used to determine the c k and c m The interaction function h(ρ A ,ρ B ) is the volume-based field integral:

[0068]

[0069] Next, we implement natural coordinates and use bilinear interpolation to evaluate the value of k ,c m The function h(ρ A ,ρ B ). This allows an exact analysis of the integral equation (1) within the defined domain.

[0070]

[0071] J k and J m is the Jacobian determinant and the shape function of the tetrahedron. k and m are the indices of the elements in objects A and B, respectively, and i and j are the vertex numbers of the two tetrahedrons under consideration. Let A and B denote the set of all elements in rigid bodies A and B, respectively. Then, the interaction function h(ρ) in the entire two-body domain can be A ,ρ B ) is expressed as:

[0072]

[0073] In step S2, the finite element method is used to perform finite element modeling on each vertex in the asteroid. The specific process is as follows:

[0074] Since each node can belong to multiple different elements at the same time, there are many redundant calculations in equation (3). By considering the mesh nodes of two objects, the summation in equation (3) can be simplified as follows:

[0075]

[0076] Here α and β are the node indices of the integrand, H αβ represents the value of the interaction function at nodes α and β, w α and w β is a custom node weight calculated as follows:

[0077]

[0078] α is the set of all elements in the grid point set A, and β is the set of all elements that contain nodes in object B. The most time-consuming part of the analysis is the calculation of the gravitational potential during integration and the calculation of the gravitational force on object A and the moments on objects A and B. The computational effort increases significantly as the size of the grid increases.

[0079] The specific calculation procedure and formula for calculating the coordinate position and interaction torque of the calculation model in the inertial system are as follows:

[0080] For any arbitrary grid point α in object A and any arbitrary grid point β in object B, the gravitational force f acting on α can be derived from the gravitational force exerted by β on α αβ and the moment p acting on the center of mass of object A αβ , we can also derive the moment p acting on the center of mass of object B from the gravitational force exerted by α on β αβ .

[0081]

[0082] According to formula (4), the gravitational force on objects A and B and the moments on the two objects can be estimated using the bilinear tetrahedron.

[0083]

[0084] According to the coordinate position of the model in the inertial system obtained in step S2, the calculated binary star motion data is imported:

[0085] In order to facilitate the derivation, the following about w is introduced α and w β New definition of:

[0086] W α =w α σ(ρ Aα ) (12)

[0087] Wβ =w β σ(ρ Bβ ) (13)

[0088] Therefore, equations (9), (10) and (11) can be rewritten as:

[0089]

[0090] The calculation process for the gravitational potential of the binary star system, the gravitational force on A, the torque on A and the torque on B is as follows:

[0091] S31, input double precision data r Ai 、r Bj , Ai , Bj , W Ai and W Bj , where W Ai Yes Ai The product of the sum, W Bj Yes Bj The product of the sum.

[0092] S32. Calculate equations (14) to (16).

[0093] S33, output result F AB 、T A and T B .

[0094] During the calculation process on a single-core CPU, the process is as follows. First, the parameters associated with the nodes in objects A and B are read iteratively. These parameters are then used to calculate f for each iteration. ij ,p ij ,q ij ,in:

[0095]

[0096] p ij =ρ Ai ×f ij (18)

[0097] q ij =ρ Bj ×f ij (19)

[0098] Then, all calculated f ij ,、p ij and q ij The values ​​are added and multiplied by G to get F AB 、T A and T BIt is worth noting that rigid body A has N1 nodes and rigid body B has N2 nodes, and it takes N1×N2 iterations to get N1×N2 f ij ,、p ij and q ij Numeric value.

[0099] The most time-consuming part of the program is to calculate the gravitational force and torque between the two binary asteroids, namely F AB 、T A and T B . w in equations (12) and (13) α and w β The terms are constants; they only need to be computed once for each node during preprocessing. ij When r Ai The term is the sum of the center of mass position vectors of object A. Calculate ρ Ai and r Bj The calculation method is similar. This step involves two addition operations of the 3D vector. Then, equation (17) is calculated, which involves one subtraction operation of the 3D vector, one modulus operation, four double-precision floating-point multiplications, and one scaling operation of the 3D vector. For p ij The calculation of Ai and f ij The process involves a cross product of the 3D vectors to calculate q i In general, each time the interaction between two objects is calculated, the summation process involves the addition of 3×(N1×N2-1) 3D vectors.

[0100] The code to implement the above steps is as follows:

[0101] First, import the binary asteroid data calculated previously:

[0102] %% Load data

[0103] % Load model

[0104] AlphaFe=loadFe("didymosa.fe");

[0105] BetaFe=loadFe("didymosb.fe");

[0106] % Load binary star system motion data

[0107] FEM=loadBinaryFEM("BinaryOutput.bt");

[0108] Then, the finite element method is used to model each vertex in the asteroid and calculate its coordinate position in the inertial system:

[0109] % Calculate the position and attitude of the binary star in the inertial system () Use the Quaternion function to create and operate quaternions and create the rotation reference point Orien array

[0110] Alphaq=quaternion(FEM.AlphaOrien);()

[0111] Alphaq_1=quaternion([FEM.AlphaOrien(:,1),-FEM.AlphaOrien(:,2:4)]);

[0112] Betaq=quaternion(FEM.BetaOrien);

[0113] Betaq_1=quaternion([FEM.BetaOrien(:,1),-FEM.BetaOrien(:,2:4)]);

[0114] %% Video demonstration results

[0115] % Coordinate system data

[0116] lA_axis=1e3;

[0117] OA=FEM.AlphaPos(1,:); % Initial position of Alpha center of mass in inertial system

[0118] XA=[lA_axis,0,0];%Alpha rigid body axis unit direction vector under fixed connection

[0119] YA=[0,lA_axis,0];

[0120] ZA=[0,0,lA_axis];

[0121] lB_axis=0.5e3;

[0122] OB=FEM.BetaPos(1,:); % Initial position of Beta center of mass in inertial system

[0123] XB=[lB_axis,0,0];%BetaUnit direction vector of the coordinate axis under the fixed connection

[0124] YB=[0,lB_axis,0];

[0125] ZB=[0,0,1B_axis];

[0126] OAt=FEM.AlphaPos;%XAt=pointRotationQ(XA,Alphaq,Alphaq_1); %Unit direction vector of Alpha rigid body axis at time t in inertial system

[0127] YAt=pointRotationQ(YA,Alphaq,Alphaq_1);

[0128] ZAt=pointRotationQ(ZA,Alphaq,Alphaq_1);

[0129] OBt=pointRotationQ(OB,Betaq,Betaq_1); %Beta center of mass position at time t in inertial system XBt=pointRotationQ(XB,Betaq,Betaq_1); %Beta rigid body axis unit direction vector at time t in inertial system

[0130] YBt=pointRotationQ(YB,Betaq,Betaq_1);

[0131] ZBt=pointRotationQ(ZB,Betaq,Betaq_1);

[0132] %Vertex Data

[0133] for i=1:AlphaFe.VertNum%Import all asteroid vertices

[0134] VertTI=pointRotationQ(AlphaFe.Vertices(i,:),Alphaq,Alphaq_1);

[0135] AlphaVerticeTI=mat2cell(VertTI,ones(length(FEM.t),1),3);

[0136] ifi==1

[0137] AlphaVertice0=cell(length(FEM.t),1);

[0138] AlphaVerticeT=cellfun(@(x,y)[x;y],AlphaVerticeTI,AlphaVertice0,'UniformOutput',false);

[0139] else

[0140] % Alpha polyhedron vertex vector at time t in inertial system

[0141] AlphaVerticeT=cellfun(@(x,y)[x;y],AlphaVerticeT,AlphaVerticeTI,'UniformOutput',false);

[0142] end%

[0143] end

[0144] for i=1:BetaFe.VertNum

[0145] VertTI=pointRotationQ(BetaFe.Vertices(i,:),Betaq,Betaq_1);

[0146] BetaVerticeTI=mat2cell(VertTI,ones(length(FEM.t),1),3);

[0147] ifi==1

[0148] BetaVertice0=cell(length(FEM.t),1);

[0149] BetaVerticeT=cellfun(@(x,y)[x;y],BetaVerticeTI,BetaVertice0,'UniformOutput',false);

[0150] else

[0151] %Vertex vector of Beta polyhedron at time t in inertial system

[0152] BetaVerticeT=cellfun(@(x,y)[x;y],BetaVerticeT,BetaVerticeTI,'UniformOutput',false);

[0153] end

[0154] End

[0155] On this basis, the calculated binary star motion data is imported:

[0156] load('AlphaNodeT.mat');

[0157] load('BetaNodeT.mat');

[0158] Finally, set the input detector position parameters and speed parameters, and calculate the detector's orbit through the following code:

[0159] trajectory_satellite=zeros(num_steps,3);

[0160] for i=1:num_steps

[0161] % Calculate the force on the satellite and traverse the nodes of the binary star to calculate the attraction

[0162] forj=1:AlphaFe.NodeNum

[0163] newNode=AlphaNodeT{i,1};

[0164] r1=OAt(i,:)+newNode(j,:);

[0165] m1=dm1*m1nodeweight(j,1)*dm1*m1density(j,1);

[0166] force_satellite_a=G*m_satellite*(m1*(r1-r_satellite) / norm(r1-r_satellite)^3);

[0167] force_satellite=force_satellite_a+force_satellite;

[0168] end

[0169] forj=1:BetaFe.NodeNum

[0170] newNode=BetaNodeT{i,1};

[0171] r2=OBt(i,:)+newNode(j,:);

[0172] m2=dm2*m2nodeweight(j,1)*dm2*m2density(j,1);

[0173] force_satellite_b=G*m_satellite*(m2*(r2-r_satellite) / norm(r2-r_satellite)^3);

[0174] force_satellite=force_satellite_b+force_satellite;

[0175] end

[0176] % Update satellite speed

[0177] v_satellite=v_satellite+force_satellite*dt / m_satellite;

[0178] % Update satellite position

[0179] r_satellite=r_satellite+v_satellite*dt;

[0180] %Save the current satellite position to the trajectory array

[0181] trajectory_satellite(i,:)=r_satellite;% trajectory position

[0182] speed_satellite(i,:)=v_satellite;%satellite speed

[0183] trajectory_force(i,:)=force_satellite;% Gravitational force on the satellite

[0184] force_satellite=0;% reset gravity

[0185] End.

[0186] Therefore, the present invention adopts the above-mentioned method for designing the probe orbit of a binary asteroid system based on the finite element method, which can more accurately characterize the interaction force of the irregular binary asteroids and plan the orbit of the probe. By characterizing the position of the spacecraft orbit and the binary asteroids of the system in real time, the characteristics of the satellite orbit can be obtained, which can be applied to program simulation in probe mission planning and spacecraft orbit optimization. The present invention is used in the field of spacecraft space simulation.

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

Claims

1. A method for designing a probe orbit for a binary asteroid system based on the finite element method, characterized in that: The following steps are involved: S1. Import the calculated binary asteroid data, input the loading data into the loading model, and simulate it in combination with the binary system motion data; S2. Use the finite element method to perform finite element modeling on each vertex in the asteroid and calculate its coordinate position in the inertial system; S3. Based on step S2, the calculated binary star motion data is imported, the input detector position parameters and velocity parameters are set, and the orbit of the detector under the finite element gravitational field model is obtained by calculation.

2. The method for designing a probe orbit for a double asteroid system based on the finite element method according to claim 1, characterized in that: In step S1, the specific process design of inputting the loading data into the loading model and combining the binary system motion data for simulation is as follows: In the model, the tetrahedron c k From object A, and tetrahedron c m from object B. Let ρ A =(x a ,y a ,z a ) represents the element c k For any vertex position, let ρ B =(x b ,y b ,z b ) represents the element c m Any vertex position, where ρ A and ρ B Calculated from A and B respectively, the following equation is used to determine c k and c m The interaction function h(ρ A ,ρ B ) is the volume-based field integral: By implementing natural coordinates and using bilinear interpolation to evaluate the k ,c m The function h(ρ A ,ρ B ), allowing the analysis of integral equations within the defined domain, such as equation (2): Formula (2), J k and J m is the Jacobian, the shape function of the tetrahedron, k and m are the indices of the elements in objects A and B, respectively, and i and j are the number of vertices of the two tetrahedrons under consideration; let and Represent all element sets in rigid bodies A and B respectively; The integral of the interaction function in the entire two-body domain is expressed as:

3. The method for designing a probe orbit for a double asteroid system based on the finite element method according to claim 2, characterized in that: In step S2, the finite element method is used to perform finite element modeling on each vertex in the asteroid. The specific process is as follows: Since each node belongs to multiple different elements at the same time, there are many redundant calculations in formula (3); By considering the mesh nodes of two objects, the sum in equation (3) is simplified as follows: In formula (4), α and β are the node indices of the integrand, H αβ represents the value of the interaction function at nodes α and β, w α and w β is a custom node weight calculated as follows: In formula (5), α is the set of all elements in the grid point set A, and β is the set of all elements containing nodes in object B. The most time-consuming part of the analysis is the calculation of the gravitational potential and the calculation of the gravitational force on object A and the moments on objects A and B during the integration process. As the grid size increases, the computational workload increases significantly.

4. The method for designing a probe orbit for a double asteroid system based on the finite element method according to claim 3, characterized in that: In step S2, the specific calculation procedure and formula for calculating the coordinate position and interaction torque of the model in the inertial system are as follows: For any grid point α in object A and any grid point β in object B, the gravitational force f acting on α is obtained from the gravitational force exerted by β on α. αβ and the moment p acting on the center of mass of object A αβ , the moment p acting on the center of mass of object B is obtained from the gravitational force exerted by α on β αβ , as shown below: According to formula (4), the gravitational force on objects A and B and the moment on the two objects are estimated using the bilinear tetrahedron: According to the coordinate position of the model in the inertial system obtained in step S2, the calculated binary star motion data is imported: Introducing about w α and w β New definition of: W α =w α s(r Aα ) (12) W β =w β s(r Bβ ) (13) Equations (9), (10) and (11) can be rewritten as: The calculation process for the gravitational potential of the binary star system, the gravitational force on A, the torque on A and the torque on B is as follows: S31, input double precision data r Ai 、r Bj , Ai , Bj , W Ai and W Bj , where W Ai Yes Ai The product of the sum, W Bj Yes Bj The product of the sum; S32, calculation formulas (14), (15) and (16); S33, output result F AB 、T A and T B ; In the calculation process of a single-core CPU, the process is as follows: First, iteratively reads the parameters associated with the nodes in objects A and B, using these parameters to calculate f at each iteration. ij ,p ij ,q ij ,in: q ij =ρ Bj ×f ij (19) Then, all calculated f ij ,、p ij and q ij The values ​​are added and multiplied by G to get F AB , T A and T B ; Among them, rigid body A has N1 nodes, rigid body B has N2 nodes, and it takes N1×N2 iterations to get N1×N2 f ij ,、p ij and q ij Numeric value; The most time-consuming part of the program is to calculate the gravitational force and torque between the two binary asteroids, namely F AB , T A and T B ; w in equations (12) and (13) α and w β The term is a constant that is computed once for each node during preprocessing; Then, calculate equation (17), for p ij Calculation, only calculate ρ Ai and f ij The cross product of is as shown in formula (18).

Images