The Coupling Orbit–Attitude–Structure Evolution of Rubble-Pile Asteroid with Earth Flyby in the Restricted Three-Body Problem

Abstract

Some asteroids flying close to Earth may pose a threat of impact. Among them, the structural and dynamical characteristics of rubble-pile asteroids can be changed because of the tidal force of the Earth in this process. This can provide key information for predicting the dynamical evolution of potentially hazardous asteroids. In this study, the long-term evolution of the coupling orbit–attitude–structure of these small bodies is presented numerically based on the integration of two models. One is the 3D discrete element method, which models the structure and irregular shape of the rubble-pile asteroid. The other is the dynamical model of the circular restricted three-body problem (CRTBP). This provides a more precise dynamical environment of the asteroid orbital deflection, morphological modification, and attitude angles analysis compared to the frequently adopted two-body problem. Parametric studies on the asteroid evolution were performed focusing on its flyby distance and the bulk porosity. Numerical results indicate that the Earth flyby can form different patterns of modification of asteroids, where the rubble-pile structure can be destructed by considering the bulk porosity. The asteroid orbital deflection and attitude variational trends are also summarized based on the simulations of multi-orbital revolutions.

1. Introduction

The impacting threat from potentially hazardous asteroids (PHAs) (https://minorplanetcenter.net//iau/Dangerous.html, acessed on 4 April 2022) has been widely recognized throughout the world in the last two decades. A list of Earth-passing asteroids has been identified [1,2]. Even though the capability of the ground-based observation has been dramatically increased, some small celestial bodies with lower albedos still cannot be recognized, such as the asteroid 2019 OK [3]. Therefore, asteroids have attracted more and more attention and become hot spots for space explorations [4,5,6,7,8,9,10,11]. If its orbit deviates slightly due to other perturbations such as the tidal effect, it is possible for it to collide with the Earth. In order to predict their future close approaches more accurately and impact possibilities, the orbital evolution of these PHAs coupled with their attitudes is indispensable. Furthermore, there is evidence that a group of asteroids are in rubble-pile structures, which are defined as unorganized collections of macroscopic particles held together by their self-gravity [12]. It is of concern that some rubble-pile asteroids may be disrupted when they cannot sustain the stress of the tidal effect from a major planet (such as Earth, Mars, or Jupiter). The fragments or the remained body with the orbit deviation due to the perturbation effect would cause new threats to Earth. As a result, it is necessary to consider the coupling orbit–attitude–structure evolution when investigating rubble-pile asteroids passing the Earth at a very near distance.

Previous works have revealed that the origin of rubble-pile asteroids comes from the catastrophic disruption of the parent bodies in the main belt [13]. This indicates that they form as a result of gravitational reaccumulation, which describes the unorganized collection of macroscopic particles held together by their self-gravity [14]. Walsh and Richardson [15] simulated the tidal disruption of rubble-pile asteroids by gravitational aggregates and the formation of binaries during close encounters with Earth. In the above work, the efficient and commonly used discrete element method (DEM) was established to study the origin and evolution of rubble-pile bodies, which were usually modeled as gravitational aggregates [16]. For instance, Cheng and Yu [17] investigated the common formation mechanism of a top-like shape, such as the asteroids (162173) Ryugu and (101955) Bennu, by performing a numerical simulation of YORP spin-up based on the soft-sphere discrete element method.

When small celestial bodies closely approach the Earth, their morphological modifications caused by the tidal effect highly depend on their initial conditions and structural characteristics. A critical distance was derived between small bodies and the flyby planet, which is called the Roche limit [18]. Holsapple and Michel [19,20] gave a typical theory about the Roche limit, which is applied to the event of 99942 Apophis approaching the Earth in 2029. Tóth et al. [21] investigated the frequency of the NEAs’ tidal disruption within the Earth Roche limit. Sometimes, small bodies may disintegrate into small objects due to the tidal effect of a major planet [22,23]. Such conditions have been previously investigated by focusing on the flyby period. Nevertheless, long-term coupling orbit–attitude propagations with multiple-orbit revolutions were rarely studied by taking into account the tidal effect.

To understand the process of tidal encounters, Sharma et al. [24] investigated the planetary flybys of asteroids by utilizing an approximate volume-averaged method. This method offers a flexible approach toward studying the rotational dynamics of ellipsoids and shape deformation. Yu et al. [25] investigated the surface evolution of Apophis by analyzing regolith material sandpile landslides during tidal encounters. The results indicate that the small-scale avalanches are distributed on the surface of the sandpile widely, and are only triggered by tidal external perturbations. Kim et al. [26] considered the elongation, spin speed, and orbital conditions in order to analyze the sensitivity surface of an asteroid under a tidal effect during a distant planetary encounter beyond the Roche radius. The asteroid is assumed to be rigid and its shape modification is ignored, which is reasonable when the flyby distance is relatively far. Zhang and Michel [27] studied the dynamical evolution of rubble-pile asteroids during the close approach to Earth. They meticulously classified the three typical scenarios into deformation, mass shedding, and disruption, where the effects of different parameters were investigated, such as encounter conditions and material strength. These aforementioned works only focused on a segment (or a short arc) of the asteroid orbit, which covers only the close-approach phase near the Earth. Thus, it cannot capture the relatively long-term evolution of an asteroid when considering its coupling orbit–attitude–structure variations.

In fact, regarding the multiple encounters with Earth at short distances, an accurate prediction of the asteroid orbit evolution should at least consider the solar gravity at the same time, along with the asteroid structure and attitude. In such a case, the orbital dynamics of the asteroid should be derived in the Sun–Earth system instead of the simple two-body problem. Therefore, the motivation of this study is to present a more accurate prediction model for the multiple-flyby cases of asteroids over the Earth. In this study, the asteroid orbit–attitude–structure long-term evolution is investigated with the tidal encounters. Some basic points are summarized at the beginning of this paper: (1) the coupling orbit–attitude dynamics of the asteroid are always valid throughout the whole-period propagation, where the Sun–Earth–asteroid system is modeled as the circular restricted three-body problem (CRTBP); (2) the possible shape deformation of the asteroid is assumed to only occur during the short-period of flyby, which will be specified in Section 3.1; (3) among millions of asteroids, the group of elongated ellipsoid asteroids was selected as the representative in this study, whose shape is similar to those of 216 Kleopatra, 1620 Geographos, and 2063 Bacchus. Other irregular shaped small bodied can be investigated in future works based on the current method. Although the discrete element model and the CRTBP model have been widely used in previous studies, their integration for the asteroid multiple-encounters with the Earth is still novel, since it provides a more accurate basic dynamical model. Additionally, some new insights on the asteroid porosity influencing the shape deformation can be obtained based on our numerical simulation.

This paper is structured as follows. Section 2 presents some basic models of this problem, including the Sun–Earth–asteroid CRTBP model and the attitude dynamical model of an asteroid in quaternion. The rubble-pile asteroid is simulated by using the discrete element method, where a MATLAB code, MatDEM [28] is used in this study. Section 3 summarizes our simulation settings and results. The investigated range of flyby to the Earth is determined by the Roche limit and Lagrangian equilibrium points. In particular, the bulk porosity of the rubble-pile asteroids is investigated when discussing their structure deformation of the asteroid. Some classical indicators, such as the semi-major axis, eccentricity, and inclination of the asteroid, are adopted to illustrate the variational trend of the asteroid state. Some discussions and the conclusion are, respectively, made in Section 4 and Section 5.

2. Methodology

The polyhedron gravitational model [29,30,31], and the finite element method [32] have been widely used in studying the surrounding particle dynamics around irregular asteroids. Nevertheless, when investigating the rubble-pile structure deformation of asteroids, the discrete element method (DEM) is usually adopted to deal with the interaction of the granular systems [33]. In the current DEM method, the assemblies of particles consisted of a series of spheres that were simulated by integrating the motion of each particle by using Newtonian equations of motion. Several codes have been developed for simulating the granular dynamics, where the well-known one in planetary science field is the parallel N-body code pkdgrav [34,35]. Hard sphere and soft sphere models are both featured in pkdgrav, where the hard sphere model is appropriate for high-speed impact and dilute particle flows, and soft sphere model is appropriate for dense particle flows. As an alternative DEM model, the asteroid in this paper was modeled as a self-gravitational aggregate by using the opensource MATLAB code MatDEM [28], in which, the interactions of particles were handled by using the soft sphere model. Given that most of the particles are in quasi-static states and long-range interactions during the tidal encounter, this model is reasonable for studying the evolution of the asteroid. More detailed comments on MatDEM can be found in Ref. [36].

2.1. Equations of Motion

The motion of an individual particle i in a particle swarm system is governed by Newton’s second law:

$$m_i\ddot{{\mathit{r}}_i}=\sum _{j=1}^{N_C}{\mathit{f}}_{ij}^{\left(c\right)}+\sum _{i=1,j\ne i}^N{\mathit{f}}_{ij}^{\left(g_1\right)}+\sum _{i=1}^N{\mathit{f}}_i^{\left(g_2\right)}$$

(1)

where $m_i$ is the mass of particle i, ${\mathit{r}}_i$ is the position vector of particle i, ${\mathit{f}}_{ij}^{\left(c\right)}$ is the contact force of particle j exerted on particle i, ${\mathit{f}}_{ij}^{\left(g_1\right)}$ is the gravitational pull of particle j on particle i, and ${\mathit{f}}_i^{\left(g_2\right)}$ is the gravitational pull from other major bodies, such as the Sun and the Earth. The symbol N is the total number of particles, and $N_C$ is the number of particles in contact with particle i, which is also referred to as coordination number.

Correspondingly, the attitude dynamical model of the asteroid is established as:

$$\mathit{I}·\dot{\mathbf{\Omega }}=-\mathbf{\Omega }\times (\mathit{I}·\mathbf{\Omega })+{\mathit{T}}_g$$

(2)

where $\mathit{I}$ is the principal moment of inertia, $\Omega$ is the angular velocity of the small celestial body rotating relative to the its principal axis of maximum inertia, and ${\mathit{T}}_g$ is the gravitational gradient torque. The torque ${\mathit{T}}_g$ can be obtained according to the following equation:

$${\mathit{T}}_{\mathrm{g}}=\mathbf{C}\left[\sum _{i=1,j\ne i}^N\left({\mathit{r}}_i-{\mathit{r}}_{\mathrm{c}}\right)\times \left({\mathit{f}}_{ij}^{\left(g_1\right)}+{\mathit{f}}_i^{\left(g_2\right)}\right)\right]$$

(3)

where ${\mathit{r}}_{\mathrm{c}}$ is the center of mass of the particle system, and $\mathbf{C}$ is the rotation matrix from the common center of mass of the Earth–Sun coordinate system to principal axis of inertia body-fixed coordinate frame. The attitude evolution of the asteroid is tracked by using the quaternion $\mathit{q}={\left[q_0,q_1,q_2,q_3\right]}^{\mathrm{T}}$. The quaternion matrix differential equation is given by:

$$\left[\begin{array}{c}{\dot{q}}_0\\ {\dot{q}}_1\\ {\dot{q}}_2\\ {\dot{q}}_3\end{array}\right]=\frac{1}{2}\left[\begin{array}{cccc}0& -{\Omega }_x& -{\Omega }_y& -{\Omega }_z\\ {\Omega }_x& 0& {\Omega }_z& -{\Omega }_y\\ {\Omega }_y& -{\Omega }_z& 0& {\Omega }_x\\ {\Omega }_z& {\Omega }_y& -{\Omega }_x& 0\end{array}\right]\left[\begin{array}{c}{q}_0\hfill \\ q_1\hfill \\ q_2\hfill \\ q_3\hfill \end{array}\right]$$

(4)

The attitude angles can be obtained by solving Equation (4). Then, the rotation matrix $\mathbf{C}$ can be obtained by using the quaternion $\mathit{q}$.

2.2. Contact Force

The contact force ${\mathit{f}}_{ij}^{\left(c\right)}$ between two particles comprises the normal and tangential components, which are denoted by ${\mathit{f}}_n$ and ${\mathit{f}}_s$. The normal force ${\mathit{f}}_n$ between the particles is modeled by linear spring and it is defined as:

$${\mathit{f}}_n=k_n{x}_n\widehat{\mathit{n}}$$

(5)

where $k_n$ is the normal stiffness, $x_n$ is the normal relative displacement, and $\widehat{\mathit{n}}$ is the normal unit vector at the contact point. When the displacement $x_n$ exceeds the breaking displacement $x_b$, the spring breaks and the connection will consequently disappear. Therefore, the maximum normal force shows:

$${\mathit{f}}_{n max}=k_n{x}_b\widehat{\mathit{n}}$$

(6)

where the shear force ${\mathit{f}}_s$ in the tangential direction was simulated by using the rupturable spring and can be expressed as:

$${\mathit{f}}_s=k_s{x}_s\widehat{\mathit{t}}$$

(7)

where $k_s$ is the tangential stiffness, $x_s$ is the tangential relative displacement, and $\widehat{\mathit{t}}$ is the tangential unit vector at the contact point. There is a maximum shear force ${\mathit{f}}_{s max}$ allowed by Mohr–Coulomb friction, which shows that:

$${\mathit{f}}_{s max}={\mathit{f}}_{s_0}-{\mu }_{\mathrm{P}}{\mathit{f}}_n$$

(8)

where ${\mathit{f}}_{s_0}$ is the initial shear resistance and ${\mu }_{\mathrm{P}}$ is interparticle friction coefficient.

2.3. Gravity Computation of Self-Gravitational Aggregate

As shown in Equation (1), the gravitational terms of the asteroid can be categorized into two components: the gravitational pull from internal particles and the gravitational pull from other major bodies. Since the rubble-pile asteroid is formed by the gravitational aggregates, there exists non-negligible self-gravity between internal particles, i.e., ${\mathit{f}}_{ij}^{\left(g_1\right)}$ term in Equation (1), which yields

$${\mathit{f}}_{ij}^{\left(g_1\right)}=\frac{G{m}_i{m}_j}{{∥{\mathit{r}}_{ij}∥}^3}{\mathit{r}}_{ij}$$

(9)

where $m_i$ and $m_j$ are the respective mass of the particles i and j, and ${\mathit{r}}_{ij}$ is the vector of particle i to particle j. The symbol G denotes the gravitational constant.

When the asteroid is in orbit, computations of gravity from other major planets ${\mathit{f}}_i^{\left(g_2\right)}$ on the asteroid can be divided into two cases according to its distance from the Earth. The influence sphere of the Earth is introduced here as partition condition, as illustrated in Figure 1. When the asteroid is out of the influence sphere of the Earth, only the attraction from the Sun is considered and the asteroid–Sun system can be modeled by two-body problem. However, when the asteroid is within the influence sphere, it means that the asteroid is in the vicinity of the Earth and the tidal effect from the Earth should be involved. Therefore, the asteroid–Sun–Earth system was modeled as circular restricted three-body problem (CRTBP). The radius of the influence sphere was specified by the collinear Lagrange equilibrium points $L_1$ and $L_2$ in circular restricted three-body problem (CRTBP). In particular, the distance from both $L_1$ and $L_2$ points to the Earth have the same values in CRTBP. Moreover, the minimum fly distance from the Earth is also constrained by the Roche limit. Once the flyby distance is less than the Roche limit, the rubble-pile asteroid could disrupt catastrophically. Such a case was excluded in the current study. In summary, when the asteroid is in the vicinity of the Earth, the fly distance is constrained between the Roche limit and the influence sphere, which is plotted in gray color in Figure 1.

Therefore, the calculation of gravity from major planets can be summarized as:

$${\mathit{f}}_i^{\left(g_2\right)}=\left\{\begin{array}{cc}-\nabla U_1\left(\mathit{r}\right)-2\mathit{\omega }\times \dot{\mathit{r}}-\mathit{\omega }\times (\mathit{\omega }\times \mathit{r}),& \parallel \mathit{r}\parallel \le ∥{\mathit{r}}_{\mathrm{L}}-{\mathit{r}}_{\mathrm{Roche}}∥\hfill \\ -\nabla U_2\left(\mathit{r}\right)& ,\parallel \mathit{r}\parallel >∥{\mathit{r}}_{\mathrm{L}}-{\mathit{r}}_{\mathrm{Roche}}∥\hfill \end{array}\right.$$

(10)

For the first case $\parallel \mathit{r}\parallel \le ∥{\mathit{r}}_{\mathrm{L}}-{\mathit{r}}_{\mathrm{Roche}}∥$, the asteroid–Sun–Earth system is solved by CRTBP. The Sun and the Earth are regarded as point masses and they are in a circular orbit rotating around their mutual center of mass with angular velocity $\mathit{\omega }$. Therefore, the gravitational force can be evaluated by:

$$-\nabla U_1\left(\mathit{r}\right)=-\frac{G{m}_{\mathrm{E}}{\mathit{r}}_{\mathrm{E}}}{{∥{\mathit{r}}_{\mathrm{E}}∥}^3}-\frac{G{m}_{\mathrm{S}}{\mathit{r}}_{\mathrm{S}}}{{∥{\mathit{r}}_{\mathrm{S}}∥}^3}$$

(11)

where ${\mathit{r}}_{\mathrm{E}}$ and ${\mathit{r}}_{\mathrm{S}}$ are the respective position vectors from the Earth and the Sun to the asteroid, and $m_{\mathrm{E}}$ and $m_{\mathrm{S}}$ are the respective mass of Earth and Sun. The angular velocity $\mathit{\omega }$ is given by:

$$\parallel \mathit{\omega }\parallel =\sqrt{\frac{G\left(m_{\mathrm{E}}+m_{\mathrm{S}}\right)}{d^3}}$$

(12)

where d is the distance between the Sun and the Earth.

When the asteroid leaves the influence sphere of the Earth, it corresponds to the second scenario, i.e., $\parallel \mathit{r}\parallel >∥{\mathit{r}}_{\mathrm{L}}-{\mathit{r}}_{\mathrm{Roche}}∥$. Only the gravity of the Sun was taken into account; thereby, the gravitational acceleration is formulated by:

$$-\nabla U_2\left({\mathit{r}}^{\mathrm{new}}\right)=-\frac{G{m}_{\mathrm{S}}{\mathit{r}}_{\mathrm{S}}^{\mathrm{new}}}{{∥{\mathit{r}}_{\mathrm{S}}^{\mathrm{new}}∥}^3}$$

(13)

Note that the coordinate systems in Equations (11) and (13) are different. The mass center of the coordinate system in Equation (11) is located at the common center of mass of the Earth–Sun system, whereas mass center of the coordinate system in Equation (13) is located at the mass center of the Sun. Thus, the coordinate transformation should be performed as:

$${\mathit{r}}^{\mathrm{new}}=\mathit{r}+{\mathit{r}}_0$$

(14)

where ${\mathit{r}}_0$ is the distance between the Sun and the common mass center of the Sun–Earth system.

3. Numerical Simulations

3.1. Tidal Encounter Simulation Setting

Numerical simulations are performed in Section 3, where some typical cases are regarded by considering the tidal effect. The material parameters of the asteroid and the initial conditions of the asteroid CRTBP are summarized in

Table 1. The physical and state parameters of asteroid.

Parameters	Symbol	Value
Asteroid size	$d_{\mathrm{L}}\times d_{\mathrm{W}}\times d_{\mathrm{H}}$	540 m × 270 m × 210 m
Young’s modulus	E	30 Mpa
Poisson’s ratio	$\epsilon$	0.15
Coefficient of friction	${\mu }_i$	0.6
Bulk density	${\rho }_a$	$2.0 \mathrm{g}/{\mathrm{cm}}^3$
Initial bulk porosity	p	0.4
Number of particles	N	5942
Particle radius	$r_{\mathrm{P}}$	25 m
Initial distance	$d_0$	$1.25\times {10}^{11}$ m

. In this simulation, the parameters of distance were described by multiples of the Earth radius $R_{\oplus }=6371\mathrm{km}$. The rubble-pile asteroid was modeled as an elongated ellipsoid granular aggregates consisting of 5942 spherical particles with a size of 540 m × 270 m × 210 m. The size of the asteroid refers to the 25143 Itokawa [37] and the value of mechanical properties was set by measuring the strength of meteorites [38]. In the flyby process, the rubble-pile asteroid could be elongated or disrupted due to the tidal force from the star (Sun). Zhang et al. [22] investigated the tidal disruption process and fragmentation at three minimum flyby distances. It showed that, when the distance is greater than 6 × ${10}^8$ m, the structure of the rubble-pile body can remain intact. Therefore, the initial distance $d_0$ between the asteroid and the common mass center of the Sun–Earth system in Table 1 was set as 1.25 × ${10}^{11}$ to prevent the shape modification of the asteroid caused by the tidal disturbance from the Sun. Note that the initial velocity of the asteroid at the beginning position is given according to the simulated flyby distance, which will be specified in Section 3.2.

It is also necessary to specify the flyby distance q of the asteroid; in which case, the tidal force from the Earth should be considered. Once the flyby distance is less than the Roche limit, the rubble-pile asteroid could disrupt catastrophically. Such a case was excluded from the current study. Therefore, the lower limit of the flyby distance was taken to be the Roche limit, whose formula yields [39]:

$${\mathit{r}}_{\mathrm{Roche}}=2^{\frac{1}{3}}{\left(\frac{{\rho }_m}{{\rho }_a}\right)}^{\frac{1}{3}}{R}_{\oplus }$$

(15)

where ${\rho }_a$ is the asteroid’s bulk density, and ${\rho }_m$ is the Earth’s bulk density. The parameters brought in Table 1 can calculate the value of the Roche limit, which is ${\mathit{r}}_{\mathrm{Roche}}\approx 1.76R_{\oplus }$. In the following simulations, the lower limit of the flyby distance was set by $0.8{\mathit{r}}_{\mathrm{Roche}}\approx 1.4R_{\oplus }$, and this value can avoid the catastrophic disruption of the rubble-pile asteroid. The upper limit of the flyby distance was set by Lagrange collinear equilibrium points $L_1$ and $L_2$ [40,41]. The $L_1$ and $L_2$ points are both approximately $1.4911551\times {10}^9\mathrm{m}$ away from the Earth by calculating the algebraic equation in CRTBP [42,43].

As observed from previous studies, particles from rubble-pile asteroids could be ejected due to the effect of the tidal force. In order to investigate the shape deformation of the main body of the asteroid, the remaining majority of the particles will be focused on. This means eliminating the ejected particles when calculating the gravitational gradient torque and the attitude of the asteroid. Hence, each element needs to detect its contact force during the numerical iteration. The handling principles for specific scenarios are divided into three kinds:

(1)

Each particle of all N elements is in contact with its neighboring particles without ejected ones. At this time, dynamical behaviors of all particles are integrated normally using the theories of Section 2.2.

(2)

If there exist $k(1\le k\le N)$ ejected particles that are completely free from contact with others, these particles will not be considered in the integration of the gravitational gradient torque exerted on the main asteroid body.

(3)

Another possible case is that $1\sim k(1\le k\le N)$ particles that make contact with each other do not make contact with other $(k+1)\sim N$ particles. In such a case, the rubble-pile asteroid can disrupt into two or more parts. The problem concerns the parent body being disrupted into a few child bodies, which is not considered in the current work.

3.2. Evolution of Rubble-Pile Asteroid’s Structure

A tidal encounter may cause the asteroid’s structure deformation, surface reshaping, and disruption. This subsection aims to investigate how the variations in the flyby distance and porosity of the asteroid affect its structure evolution during the tidal encounter. In the first set of the numerical simulation, three different flyby distances, namely, 1.4$R_{\oplus }$, 2$R_{\oplus }$, and 2.5$R_{\oplus }$, were considered by adjusting the initial velocity of the asteroid ${\mathit{v}}_0$ to govern. The magnitude of initial velocities corresponding to the above flyby distance were, respectively, 4.0 km/s, 5.2 km/s, and 6.3 km/s. The following presents the detailed tidal encounters outcomes for the morphological evolution of rubble-pile bodies with every condition.

In Figure 2, the snapshots of rubble-pile asteroid evolution during the tidal encounters at different flyby distances are shown. The first row presents the $q=2.5R_{\oplus }$ case with blue flowchart arrows in Figure 2. The result shows that there only exists a small amount of particle sputtering, and that the shape of the rubble-pile asteroid does not alter much compared to the origin. At the flyby distance $q=2R_{\oplus }$ in the second row with yellow arrows, the structure of the rubble-pile asteroid is significantly affected by tidal effects. The main body is elongated and the particles shed from both sides of the asteroid. The condition $q=1.4R_{\oplus }$ is shown in the final row with green arrows. This flyby distance is in the critical upper limit based on the previous discussions in Section 3.1. The amount of shedding particles from both sides of the rubble-pile asteroid increases. The main body is more elongated than in other conditions and the structure approaches the disruption when compared with the other two conditions. Comparing the three outcomes, the nearer the approach to the Earth, the stronger the effect of the tidal encounter for the rubble-pile asteroid. The tidal encounter may trigger three cases for the asteroid structural evolution: the surface particles sputtering, the main body elongated deformation, and the particles shedding. Generally, these three events are not independent, and they often occur at the same time. Therefore, it can be seen that the influence of the flyby distance on the structural evolution of rubble-pile asteroids is significant under the tidal force.

In addition to flyby distance, the physical parameters of the rubble-pile asteroid, such as porosity, also play an important role in structural evolution. For instance, in the Double Asteroid Redirection Test (DART) mission, it was reported that different material porosities may lead to different asteroid responses after the impact is implemented [44,45,46]. For the second set of the numerical simulation, the porosity of rubble-piles was varied in order to discuss the structural evolution during the tidal encounter.

The asteroids in the Solar System usually have a large scope of porosities, and even up to $50\%$ in the case of asteroid 253 Mathilde [47]. If ${\rho }_A$ defines the bulk density and ${\rho }_C$ defines the grain density, the porosity shows that:

$$p=1-\left(\frac{{\rho }_A}{{\rho }_C}\right)$$

(16)

Generally, the high porosity illustrates that the asteroid is packed tightly. In this simulation, the initial position and flyby distance were set the same as in Section 3.1 to guarantee that the other factors are unchanged. The minimum flyby distance was set to $q=2R_{\oplus }$ and the porosities were assigned three different values, namely 0.3, 0.4, and 0.5.

In Figure 3, the snapshots of different-porosity rubble-pile asteroids in the evolution of the tidal effect are shown. This simulation reveals the effect of the physical parameters of the asteroid in the evolution of tidal encounters. By comparing the three scenarios, the results show that a lower porosity leads to more sputtered particles distributed around the asteroid. As the porosity increases, the particles are basically no longer shed. For the porous material, the stress on the surface from the tidal force is used to compact the pores and more energy is scattered, making it less exploitable to set the particle material in motion and eject material. Therefore, the tidal stress can impose the motion of the surface of the asteroid, and the existence of porosity acts as a buffer during tidal encounters.

3.3. Effect of Flyby on the Asteroid Orbit Evolution

The previous section analyzes the structural evolution of the asteroid after its encounter with the Earth. In the process of a close approach to the Earth, the orbital dynamic properties of the rubble-pile asteroid may also be modified due to the tidal effect. This subsection summarizes the multi-orbital propagation of the asteroid after it is perturbed by the tidal force of the Earth. The initial condition and parameters can be seen in

Table A1. The initial conditions and parameters for the numerical simulations.

Model Parameters	Symbol	Value
Asteroid size	$d_{\mathrm{L}}\times d_{\mathrm{W}}\times d_{\mathrm{H}}$	540 m × 270 m × 210 m
Young’s modulus	E	30 Mpa
Poission’s ratio	$\epsilon$	0.15
Coefficient of friction	${\mu }_i$	0.6
Bulk density	${\rho }_a$	2.0 $\mathrm{g}/{\mathrm{cm}}^3$
Number of particles	N	5942
Particle radius	$r_{\mathrm{P}}$	25 m
Bulk porosity	p	0.3
		0.4
		0.5
State parameters	Symbol	Value
Initial distance	$d_0$	$1.25\times {10}^{11}$ m
Initial velocity	$v_0$	4.0 km/s
		5.2 km/s
		6.3 km/s
Flyby distance	q	1.4 $R_{\oplus }$
		2.0 $R_{\oplus }$
		2.5 $R_{\oplus }$

, Appendix A, where the flyby distance and bulk porosity of the asteroid are taken as the variables to investigate the orbital evolution. To clearly describe the orbital evolution during the tidal encounters, a set of comparative simulations were added. It is a two-body problem of an asteroid–Sun system, which ignores the influence of Earth. The orbital elements, which contain the semi-major axis a, eccentricity $\mathit{e}$, and orbital inclination i, were investigated as the indicators of orbital dynamic properties. The semi-major axis a and the eccentricity vector $\mathit{e}$ show that:

$$a={\left(\frac{2}{\parallel \mathit{r}\parallel }-\frac{{\parallel \mathit{v}\parallel }^2}{{\mu }_g}\right)}^{-1}$$

(17)

$$\mathit{e}=\frac{1}{{\mu }_g}[\mathit{v}\times (\mathit{r}\times \mathit{v})]-\frac{\mathit{r}}{\parallel \mathit{r}\parallel }$$

(18)

where $\mathit{r}$ and $\mathit{v}$ are the vectors of position and velocity. The ${\mu }_g$ is the gravitational constant of the center body and means the multiplication of the gravitational constant G and the mass of the Sun. The semi-major axis of an elliptical orbit is a finite positive number and the range of eccentricity $\parallel \mathit{e}\parallel$ is greater than 0 and less than 1. If the angular momentum and the unit vector in z-axis are known, the orbital inclination i is calculated by the next formula:

$$i=arccos\frac{\mathit{z}·\mathit{h}}{\parallel \mathit{h}\parallel }$$

(19)

where $\mathit{z}$ is the unit vector in the z-axis and $\mathit{h}$ is the angular momentum of the asteroid.

Here, the simulated results are presented in Figure 4, which contains the trajectory of an asteroid in the orbit and the curves of orbital elements over time. In Figure 4, the minimum flyby distance q is set as 1.4$R_{\oplus }$, 2$R_{\oplus }$, and 2.5$R_{\oplus }$, respectively, while the porosity of the asteroid remains 0.4 for all three scenarios. The red curves in Figure 4 presents the asteroid orbit evolution when the tidal perturbation from the Earth is ignored, and only the attraction from the Sun is considered. It serves as a reference to clarify how the Earth tidal force affects the asteroid orbit propagation. From Figure 4a–c, it can be concluded that the orbit evolutions are quite different with/without the Earth tidal force for all three scenarios. The graph of the orbital radius, which is the distance between the center of mass of the asteroid and the origin of coordinates at each moment, can describe the orbital deviation in Figure 4d. With the minimum flyby distance increasing, the Earth tidal effect becomes weaker. This leads to a decrease in the orbital radius curves, and the shape of the curves are close to the two-body problem.

Figure 5 presents the orbital elements, including the semi-major axis, eccentricity, and orbital inclination. For the strong encounters, where the minimum flyby distance is close to the $0.8r_{\mathrm{Roche}}\left(q=1.4R_{\oplus }\right)$, the curves of the semi-major axis show the upward and periodic trend. When the tidal encounters gradually weaken, such as the $q=2R_{\oplus }$ and $q=2.5R_{\oplus }$, the trend of curves tends to be stable and periodic. In the same way, it can be revealed that Earth’s tidal force can affect the orbital eccentricity. It shows a mutation in the middle of an orbital period, as indicated by the blue curve in Figure 5b. As shown in Figure 5c, for each run time, this simulation also gives four curves about the orbital inclination. The Earth tidal effect leads to orbital deflection, which is relative to the equatorial plane. The orbital inclination of a two-body is less than a three-body at the same time. These outcomes imply that, for the different intensities of tidal encounters, such as the flyby distances, the strong tidal effect usually triggers a greater orbital evolution and the structural modification is more obvious; refer to Section 3.2. Accordingly, the effects of structural modification and orbital transmission during the tidal encounter are coupled with each other.

Besides the minimum flyby distance, the porosity of the rubble-pile is another significant parameter for the orbital transmission under the Earth’s tidal force. This paper analyzes the evolution of the rubble-pile asteroid structure in Section 3.2. According to Equations (2)–(4), the motion of orbit for an asteroid is strongly affected by gravitational gradient torques. The porosity p in this simulation was set as 0.3, 0.4, and 0.5, respectively, while the minimum flyby distance of the asteroid remained $q=2R_{\oplus }$ for all three scenarios. The curves of orbital evolution are similar to Figure 5a–c, and the orbital radius indicating that different porosities can affect orbit deflection to some extent is shown in Figure 6. For high-porosity conditions, such as $p=0.5$, the curve of the orbital radius changing with time is shown in Figure 6a, which presents a complete periodicity, and the orbital radius does not change too much with the increase in the orbital period. This outcome is similar to the event of weak tidal encounters. For the lower-porosity conditions, such as $p=0.3$ and $p=0.4$, the deflection of the orbital radius shows an upward trend, which, eventually, stabilizes. The orbital elements, including the semi-major axis, eccentricity, and orbital inclination, are also discussed. The semi-major axis for different material porosities is shown in Figure 6b. It appears that the change in the semi-major axis of the asteroid with high-porosity material is stable. On the contrary, at the condition of $p=0.3$, there exist fluctuations of curves at the end of an orbital period due to the tidal effect. In addition, there exists a difference in porosity that can be explained by analyzing the eccentricity in Figure 6c. The intrinsic elliptical orbits are still similar to Figure 5b because the range of eccentricity is 0.025∼0.5. Comparing the conditions of three kinds of porosities, a lower porosity has a larger eccentricity. This is because shape deformation causes the deflection of the orbits. The orbital inclination has a rising trend with time increasing, as shown in Figure 6d. For $p=0.3$, it indicates that the value of degrees is more than other orbits of porosity. Since the porosity of rubble-pile material can disperse the stress for the tidal effect, a higher porosity could decrease the deformation of the asteroid. Therefore, it can obtain the relationship between the material porosity and orbital transmission. In the long-term evolution of an asteroid, the higher porosity can keep the structure completely, and the orbit elements have variable forms compared to the two-body orbits.

Consequently, the tidal encounter effects are dependent on the initial conditions such as physical and state parameters. This leads to an internal structure change; for example, deformation, shedding, and even disruption. These effects gradually influence the motion of asteroids so that they are reflected in orbit elements, which research the semi-major axis, the eccentricity, and the orbital inclination in this paper.

3.4. Effect of Flyby on the Asteroid Attitude Variation

According to Equations (2)–(4), it can be seen that the gravitational gradient torques are significant for the periodic orbits, as well as the asteroid body shape and attitude. The attitude angles, including the pitch angle ($\theta$), roll angle ($\varphi$), and yaw angle ($\psi$), are solved by a quaternion, which avoids the rotation of a singular matrix. The differential Equation (4) can be solved to obtain the value of the quaternion and the attitude angles. In this simulation, the attitude of the asteroid under the Earth’s tidal force was focused on because it is necessary for the evolution of the coupling orbits–attitude–structure. The parameters, such as the minimum flyby distances, are set by $q=1.4R_{\oplus }$, $q=2R_{\oplus }$, and $q=2.5R_{\oplus }$. The pitch angles are shown in Figure 7, in which, the red, green, and blue curves represent the flyby distance of $q=1.4R_{\oplus }$, $q=2R_{\oplus }$, and $q=2.5R_{\oplus }$, respectively. As the flyby distance is dramatically increased, the amplitude of the pitch angles is amplified, but the range is limited to $-{80}^{\circ }$∼${80}^{\circ }$. Figure 8 presents the relationship between roll angles and time, where the initial conditions are similar. The roll angles oscillate under $-{40}^{\circ }$∼${40}^{\circ }$. As shown in Figure 9, a graph of yaw angles versus time and the range of angles $-{50}^{\circ }$∼${50}^{\circ }$ are presented. Comparing the three curves can show that the variable of attitude angles has a relationship with the flyby distance. For the strong tidal encounters, i.e., $q=1.4R_{\oplus }$, presented by red curves in three figures, the fluctuating trend of angles is more gentle. This kind of phenomenon indicates that the spin state of the asteroid is gradually affected with the long-term evolution of orbit. On the contrary, with the tidal effect weakening, the spin state can be inadequately changed. In the outcomes before, the tidal encounter can reshape the rubble-pile body into an extremely elongated shape. The angular momentum due to the tidal torque in this process obstructs the initial increasing rotational speed and the spin period can be extended during the tidal encounter. This is reflected in the milder amplitude of the oscillation curves, which is the relationship between three attitude angles and time. Therefore, these results indicate that the evolutions of the attitude and structure are coupled with each other.

4. Discussion

In this work, the number of rubble-pile asteroid particles was always a concern. A trade-off had to be made when selecting the number of particles. The main principle for the selection is to ensure the accuracy of the numerical simulation while ensuring that the calculation efficiency is also acceptable during the evolution of tidal encounters. To clarify the accuracy of the above results, three different groups of particles were adopted to make a comparative study. Here, the investigated particle numbers were set to 1421, 5942, and 10,902, and these comparative numerical experiments define a shedding ratio $\sigma \%$, which presents the ratio of the shedding particle numbers and the total particle numbers. Referring to Section 3.2, the simulating parameters include the minimum flyby distance $q=2R_{\oplus }$ and the porosity $p=0.4$. The outcomes are presented in Figure 10, and the simulating times for different numbers of particles are summarized in

Table 2. The simulating time of different rubble-pile particles for five orbital periods.

Particles Number	$\mathit{N}=1421$	$\mathit{N}=5942$	N = 10,902
Simulating Time/h	8.25	14.5	22.5

. In all of the simulations, the GPU model of a computer used was the Tesla P100, where the memory was 16 GB.

From the shape of the asteroid model in Figure 10a, it can be seen that $N=5942$ and N = 10,902 are more similar to an elongated ellipsoid in appearance than the case where $N=1421$. By investigating the shedding ratios of rubble-pile asteroids with three particle numbers over time in Figure 10b, the results show that the shedding ratio of $N=1421$ is higher than the other two scenarios. This causes the mass of the main body to decrease so much that it disturbs the attitude–orbit propagation of the asteroid during the long-term evolution. Nevertheless, the cases $N=5942$ and $N=10902$ have a similar shedding ratio and the computation of $N=5942$ is more efficient in this simulation. Therefore, the particle number $N=5942$ is selected for the final particle numbers of the rubble-pile in this paper.

These works are different from the earlier studies investigating a segment of the asteroid orbit during the tidal encounters. The circular restricted three-body problem (CRTBP) model was applied in this study to expand the asteroid orbit from the Earth–asteroid to the Sun–Earth–asteroid. It is in favor of analyzing the asteroid orbital deflection in multiple flybys over the Earth. The results of the numerical simulation indicate that the orbit, attitude, and structure of the asteroid are coupled with each other with the accumulated effect. It was verified by simulation experiments that the parameters of the minimum flyby distance q and porosity p can simulate different strengths of tidal encounters.

Hence, various scenarios were set up to simulate the asteroid flyby in Section 3. With the strong tidal encounter, i.e., $q=1.4R_{\oplus }$ or $p=0.4$, the structure of the rubble-pile asteroid can be modified dramatically; refer to Section 3.2. As for the orbital simulation, the structural modification leads to a large orbital deflection. Following multiple-encounters, these orbital deflections are more obvious. For the opposite condition, a weak tidal encounter, the orbital shape is close to the two-body problem and the structure has little modification. Furthermore, the outcomes of the attitude show that, when the shape of the asteroid is elongated, the magnitude of the oscillation of attitude angles can decrease. This is possible because the angular momentum from the tidal torque of Earth slows down the rotational speed and explains the relationship between the attitude and structure.

This investigation could be used to track the long-term orbit–attitude evolution of the asteroid after it is perturbed during near-Earth flybys. After the discovery of (99942) Apophis in 2004, it was initially predicted to have the possibility of impacting Earth in 2029. Recently, the 2029 impact has been ruled out by scientists. NASA plans to explore the asteroid (99942) Apophis during its encounter with the Earth as an extended mission of OSIRIS-REx [48]. The method proposed by this paper facilitates the exploration of potentially hazardous asteroids during near-Earth flybys, such as (99942) Apophis, and provides more accurate predictions.

5. Conclusions

This paper investigates the coupling orbit–attitude–structure evolution of near-Earth rubble-pile asteroids under a tidal effect from Earth in the circular restricted three-body problem (CRTBP) using the MATLAB discrete element method (MatDEM) to mimic this process. Numerical simulations were performed regarding the multiple-flyby scenarios, considering two key parameters, i.e., the porosity of rubble-piles and the flyby distance. Compared to the Sun–asteroid two-body problem, the orbit, attitude, and structure of a rubble-pile asteroid can be affected by the Earth’s tidal effect. Specifically, the Earth’s tidal effect leads to the deflection of the asteroid orbit, showing a divergent trend on the curves of the asteroid orbital elements. As for the asteroid structure, the tidal effect causes the original shape to be more elongated. Along with the increasing tidal effect, three possible patterns of structural modification can form, including deformation, particles shedding from the surface, and disruption. On the contrary, the increase in the tidal effect could decrease the attitude variation magnitude, and the variational trend is smoother. The reason is that the morphological modification of the asteroid makes the tidal torques from Earth block the spin speed from increasing. Obviously, a closer flyby distance corresponds to a strong tidal effect. In such a case, the rubble-pile of a higher porosity can allow for part of the tidal effect by preventing the structure deformation of the rubble-pile asteroid. Therefore, asteroids with a more porous structure illustrate the weaker tidal effect within a certain range. Future studies can focus on more different initial shapes of asteroids. An asteroid model that is shaped as gyroscopic (such as Ryugu and Bennu) and irregular elongated (such as Itokawa and 4179 Toutatis) could be interesting for investigations.