Numerical simulation and behavior prediction of a space net system throughout the capture process: Spread, contact, and close

In this paper, we develop an exhaustive numerical simulator for the dynamic visualization and behavior prediction of the tether‐net system during the whole space debris capture phases, including spread, contact, and close. First of all, to perform its geometrically nonlinear deformation, discrete different geometry theory is applied to model the mechanical response of a flexible net. Based on the discretization of the whole structure into multiple vertexes and lines, the internal force and associated Hession are derived in a closed form to solve a series of nonlinear dynamic equations of motion. The spread and deployment of a packaged net can be realized using this well‐established net solver. Next, a multidimensional incremental potential formulation is selected to achieve the intersection‐free boundary nonlinear contact and collision between the deformable net and rigid debris. Finally, for the closing mechanism analysis, a log‐like barrier functional is derived to achieve the nondeviation condition between the ring–rod linkage system. The C2 ${C}^{2}$ continuous log barrier functionals constructed for both the contact model and the linkage system are smooth and differentiable, and, therefore, the nonlinear net capture dynamic system can be efficiently solved through a fully implicit time integrator. Overall, as a demonstration, the whole capture process of a defunct satellite using a hexagon net is simulated through our well‐established numerical framework. We believe that our comprehensive numerical methods could provide new insight into the optimal design of active debris removal systems and promote further development of the online control of tether tugging systems.

Space debris is nonfunctional equipment in the atmosphere or the earth's orbit manufactured by humans. 1 In this century, the amount of space debris has increased dramatically with the development of launch technology; for example, based on the investigation by the US Space Surveillance Network in 2022, more than 22 000 space debris with sizes larger than 100 mm have been detected. 2These useless objects can be dangerous for both spacecraft and astronauts; for example, there may be physical contact and collision. 3As a result, there is an urgent need to address the space debris issues and motivate the demands for active debris removal.
Two types of systems are widely used for space debris capture and removal: (i) contact-based methods and (ii) noncontact methods. 4On the one hand, robotic arms and harpoons serve are two typical tools for the contact-based method, where the space debris needs to be linked to a parent satellite and contact is required during the debris capture. 5 the other hand, laser and ion beams are usually used for the noncontact method, during which physical contact is not needed. 6,7The tether-net system, different from the aforementioned technologies, has gained increasing popularity due to the following advantages: (i) it is able to travel a long distance, (ii) it can adapt to debris with different shapes, and (iii) it is cheap and lightweight. 8Due to these advantages, many studies have been carried out on tether-net systems, for example, dynamics, 9 control, 10 and experiments. 11Botta et al. 12 developed a commercial software named Vortex Studio to study the deployment dynamics of the tether-net system.Golebiowski et al. 13

used a Cosserat
Rod Theory to model the tether-net, during which the bending and twisting modes of the structure were included.Shan et al. 8 used the Absolute Nodal Coordinate Formulation to simulate the mechanics of a flexible net, and the effects of shot mass, shot speed, and shot angle were analyzed by developing a simplified order-reduced model. 1 From the viewpoint of computational modeling, the flexible tether-net dynamics used either a naive mass-spring-damper approach 12 or a sophisticated finite element method. 14However, the discrete differential geometry (DDG) formulation, initially developed for visual animation, has become popular in the engineering community in recent years [15][16][17] due to its computational efficiency and numerical robustness in handling nonlinear deformation, contact, and collision, which are important in net capture analysis.Previous DDG-based methods have shown surprisingly successful performance in simulating the nonlinear mechanics of slender structures, for example, rods, 18,19 ribbons, 20 plates, 21 shells, 22 gridshells, 23 and nets. 246][27][28] Besides the geometrically nonlinear deformation of flexible net structures, the boundary nonlinear contact between the deployment net and the target object is also a challenge for dynamic analysis.Different numerical frameworks were adopted to model the boundary nonlinear contact between a flexible net and a rigid satellite, for example, Hertzian contact theory, 29 quadratic penalty potential, 25 the modified mass method, 23 and the incremental potential (IP) formulation. 26Even though there are several numerical methods for the spread and contact analysis of the tether-net system, to the best of our knowledge, so far, the whole capture process that combines the deployment phase, contact phase, and close phase into a unit framework is missing.Nevertheless, Endo et al. achieved the closing mechanism by decreasing the natural length of a net mouth string; this formulation cannot be obtained physically and it is difficult to adopt directly in real engineering applications. 30re, in this paper, a hybrid numerical framework is established to simulate the nonlinear dynamics of a tether-net during all capture phases, which include spread, contact, and close.As shown in Figure 1A, a packaged net is deployed and moves to a defunct satellite by shooting six corner masses (shown as black dots), and the boundary rods (shown in red) for the closing mechanism are connected to a hexagon net by multiple rigid rings (shown as blue squares).To simulate such a complex system, three numerical frameworks are used together: (i) the DDG method for the The capture of a defunct satellite using a flexible net system.(B) Components of our net capture system.geometrically nonlinear deformation of the flexible net, (ii) the IP formulation for the boundary nonlinear contact between net and debris, (iii) and the reduced-order model for the complex coupling between the deformable net and the sliding ring.First of all, the tethernet system is discretized into multiple one-dimensional (1D) elements, for example, nodes and edges, and the total elastic energies are the sum of uniaxial stretching and misaligned bending, which is a geometrically exact formulation for flexible structures undergoing large deflection and rotation.The gradient vector and the Hessian matrix of the elastic energies are given in a closed form on the basis of DDG.Second, for the boundary nonlinear contact between the flexible net and the rigid satellite, the multidimensional IP contact barrier potential with C 2 continuity is used, such that its gradient and Hessian are also obtainable.Finally, to achieve nondeviation interaction between a rigid ring sliding on a flexible rod, a special log barrier functional similar to the IP formulation is constructed, which is also smooth and differentiable.Overall, a fully implicit time-marching scheme for the nonlinear elastodynamic system including net deformation, contact interaction, and ring-rod coupling is established.
As a demonstration, a comprehensive net capture process considering the spread phase, the contact phase, and the close phase is performed using our newly introduced numerical framework.
Our paper is organized as follows: The numerical formulation is detailed in Section 2. Next, Section 3 numerically explores the complex dynamic behavior during the tether-net capture process.
Finally, conclusive remarks are presented in Section 4.

| NUMERICAL METHOD
In this section, the numerical method used for the tether-net capture simulation is presented.As shown in Figure 1B, our net capture system includes the following components: (i) a locally refined hexagon net with length L = 15 m and mesh grid L Δ = 1 m, (ii) six bullet masses for deployment, (iii) six boundary rods, and (iv) multiple rings for the connection between the hexagon net and boundary rods.To simulate a complex mechanical system like this, three numerical methods are combined together.The geometrically nonlinear mechanism of a flexible net is modeled by the DDG theory, which can be used for the deployment investigation of a flying net.
Next, the boundary nonlinear contact between a flexible net and rigid debris is achieved using a multidimensional IP method.Then, the reduced-order model is introduced for the ring-to-rod coupling system, which is essential for the closing mechanics simulation of a tether-net system.Finally, the numerical simulation loop for the tether-net capture system is presented.

| Flexible net
First of all, we introduce a numerical framework for the dynamic simulation of tether-net systems using a geometrically exact formulation.We discretize the flexible net structure into N nodes to represent its deformable configuration, and, consequently, a degree of freedom (DOF) vector with the size of N 3 is produced, where 3 denotes the position of the ith node.The edge vector between the ith node and the jth node is expressed as and its tangential direction is defined as where ∘   indicates the 2 norm of a vector.The deformable potential of a net structure has two forms: (i) elastic stretching energy, E s , and (ii) elastic bending energy, E b .As shown in Figure 2A, the linear stretching energy is associated with the elongation of the ij ( )th edge, where ϵ ij is its uniaxial stretching strain, Here, a bar on the top represents the quantity in the undeformed configuration; for example,   e ¯ij is the undeformed length of the ij ( )th edge.The bending energy between the ij ( )th edge and the jk ( )th edge is given by their misalignment, where l Δ ¯ijk is its Voronoi length and can be computed as (7)   and κ ijk is an approximation of the discrete curvature | 267 where N s is the total number of stretching elements and N b is the total number of bending elements.The internal elastic force vector,

| Contact model
Next, the contact problems between the flexible tether-net and rigid space debris are described using a computationally efficient IP method. 31A nonlinear log-based barrier function with C 2 continuous is constructed to determine the potential energy between two approaching basic elements.For a complex contact interaction in real engineering applications, a combinatorial geometric computation is utilized for closed-form distance formulas, for example, edge-edge and point-triangle, as shown in Figure 3A,B, respectively.
The distance between two edge elements, x x { , } i i+1 and x x { , } j j+1 , can be formulated as a constrained optimization problem 31 and the distance between a point, x i , and a triangle mesh, First of all, when both constraints are imposed, the minimum distance between the two primitives is a point-point evaluation as follows: Second, if only one constraint is activated, the minimum distance between the two primitives is a point-edge scenario as follows: Finally, when there is no constraint, the distance computations are parallel-plane evaluations, that is, It should be noted that when two edges are nearly parallel to each other, the point-edge computation is used as a substitution to avoid numerical errors.
Thereafter, a smooth log barrier potential is constructed to simulate the penetration-free condition between ith approaching pairs as follows: where d i is its minimum distance between two primitives, K c is the barrier stiffness, and d ˜is a barrier distance.It is worth noting that the energy barrier is discontinuous if d ˜= 0.0, which is equivalent to the ideal contact in the real physical world.The normal contact force would be zero if the minimum distance between two elements is larger than d ˜and would gradually enlarge if the distance is decreased.Finally, the reaction force would be close to infinity at nearly zero distance.The contact barrier function is with C 2 continuous; thus, its gradient and Hessian matrix for implicit time marching can be derived analytically.
The overall potential for intersection-free contact is the sum of c interacting pairs constructed by continuous collision detection Similarly, the contact force, ∈  F N 3 , is the negative gradient of contact potential In the current study, we only focus on the configuration evaluation of the tether-net capture system, such that the tangential frictional force between the rigid satellite and the flying net is not included, which can be formulated based on the maximal dissipation principle. 27,31

| Ring-rod interaction
The previous IP formulation can deal with the contact interaction between a flexible ring and a rigid rod if both of them are discretized into multiple nodes and edges for contact detection and dynamic simulation.However, the redundant degrees of freedom would significantly increase the computational time.Thus, here, a reducedorder barrier energy is introduced to handle the coupling between a sliding ring on a deformable rod.
F I G U R E 3 (A) Rod-to-rod contact primitive and (B) point-totriangle contact primitive.
Here, an entangled ring-rod coupled system is considered, as shown in Figure 4, and the rod is represented by several vertexes, while the ring is represented by a single point.The contact primitive is built on the basis of the minimum distance between the ring and the rod, such that three nodes are used, ≔ ∈  x x x ( , , ) i i i r +1 9 , where x i and x i+1 are the rod nodes, and x r is the position of the moving ring.
The minimum distance r between the joint, x r , and the cable element, , is parallel to the point-edge calculation, that is, The constrained equation can be reduced either to the point-point distance given in Equation ( 13) or the point-edge evaluation formulated in Equation (14).Next, to maintain the nondeviation condition between the moving point and the elastic rod, we construct the following barrier functional: where K J is the linkage stiffness similar to K C and r ˜is a linkage parameter similar to the d ˜in Equation ( 16).Consequently, here, the ring radius is set to ideal zero, and the ring-to-rod measurement can be directly achieved by subtracting the ring radius.It is noteworthy that the rod-to-ring coupling pair i needs to be updated over time, simply because the minimum distance between the sliding ring and the elastic cable is changed.Again, the overall potential for nondeviation contact is the sum of the j coupling element ∑ E J r r r = ( ),with <˜.
Again, this newly introduced coupling functional for the ring-to-rod connection is also with C 2 continuous, and the normal contact force is close to zero when the distance is small, while the force is almost infinite when the distance reaches the threshold, r ˜.Also, the interaction force, ∈  F N 3 , is the negative gradient of the coupling potential With this barrier functional and the associated coupling force, the nondeviation condition between the sliding ring joint and the elastic cable is successfully achieved.

| Equations of motion
This section presents the time integration scheme for the nonlinear net capture system with nonpenetration contact and ring-to-rod interaction.The overall discrete system with N nodes is treated as a multibody dynamic system, with total degrees of freedom ∈  q N 3 .

Its continuous dynamic equation of motion is
where M is the mass matrix, is the Rayleigh damping matrix, and F ext is the external force vector.
Moreover, the IP energies need to be used to satisfy the nonpenetration condition between the approaching elements and the nondeviation condition between the ring-rod linkage.Continuous geometric computation is used at each time step to update the contact force vector F (for net-debris collision) as well as the coupling force vector F (for ring-to-rod sliding).
At the kth time step, t k , we use the implicit Euler method to numerically integrate the DOF vector and update it to where h is the time step size and F Cq = D is the damping force vector and is set to zero in the current numerical investigation.As the interaction potentials are based on a C 2 continuous function, such that the Newton-Raphson method can be adopted to solve the nonlinear equations, that is, the Jacobian matrix associated with Equation ( 24) can be derived, therefore, our numerical framework is fully implicit.

| RESULTS
In this section, the numerical results are generated by our discrete simulator.The proposed numerical framework is programmed within the C/C++ platform using Eigen; 33  | 269 step size for numerical integration is h = 1 ms.Even though the computational time would increase if an implicit numerical scheme is adopted, here, we still use the implicit Euler integrator due to its numerical robustness.It is worth noting that the real-time simulation can be achieved if the total degrees of freedom are less than 3000. [28]

| Spread
First of all, we present the spread behavior of a foldable net.Here, we assume that the shooting velocity of each bullet mass is v = 10.0 m∕s s and the shooting angle is ϕ = 45°.The net is initially packaged into a small net bag, which is located at the origin of coordinates; see Figure 5A.
Due to the velocity of the bullet mass, the net is stretched out of the original point and spread in space.When t < 0.8 s, the total spread area is almost zero, and only boundary rods are stretched out; however, when the physical time is larger than 1.0 s, the net can expand in space and the system is almost deployed at a time equal to 2.0 s.The overall shooting and spreading phases for illustration are shown in Figure 5A-F.

| Contact
Next, due to the existence of a defunct satellite, the flying net will come into contact with the object and, as a result, the contact and collision models need to be included.Here, to achieve the contact simulation, the contact parameters are selected as barrier distance d ˜= 0.2 m and barrier stiffness K = 1.0 N∕m c .A more accurate result can be achieved if a smaller barrier distance d ˜and a larger barrier stiffness K c are used; however, in that case, a smaller time step size would be required for the contact detection and numerical robustness.
In other words, on the one hand, for a larger clamped distance, the reaction force would become smoother, and, as a result, the simulator would perform more robustly; on the other hand, for a smaller clamped distance, the contact force is closer to the ideal contact model, that is, a sudden change, such that the numerical framework would be more precise.The contact stiffness parameter follows similar criteria.The parameter selection is a balance between robustness and accuracy and can be tuned manually based on specific needs.Specifically, the satellite is discretized into 40 triangular meshes for contact detection.
After reaching the maximal spreading area at approximately t = 2.3 s (Figure 6A), the contact between the fixed satellite and the flying net is detected and the net would relax due to the collision at t ≈ 2.5 s.
See Figure 6B for details.Interestingly, we found that the net would rebound to a lower position instead of always making contact with the satellite after the collision.The snapshots of the tether-net capture system during the contact and collision phase are shown in Figure 6A-F.Altogether, the space debris capture process using a tether-net system can be divided into three parts: (i) deployment during t 0.0s ≤ ≤ 2.2s, (ii) contact during t 2.2 s ≤ ≤ 3.6 s, and (iii) close during t 3.6 s ≤ ≤ 4.1 s.The overall dynamic rendering for the tether-net capture progress is presented in Supporting Information: Movie-S1.mov.

| CONCLUSION
In conclusion, in this paper, we have developed a comprehensive numerical simulator for the behavior prediction and dynamic analysis of the tether-net system during the whole debris capture process, which includes the spread phase, the contact phase, and the close phase.For this purpose, our well-established computational framework fully accounts for three components: elasticity of the tether net, contact between the net and the satellite, and the connection between the ring and the rod, which were achieved using the DDG theory, the IP formulation, and a reduced-order model, respectively.First of all, the whole tether-net structure was discretized into multiple nodes and F I G U R E 6 Snapshots of a tether-net capture system during the contact phase (A-F).
F I G U R E 7 Snapshots of a tether-net capture system during the close phase (A-F).
HUANG ET AL.
| 271 edges, and the deformation modes were the sum of stretching and bending.Second, the multidimensional IP formulation was selected to achieve boundary nonlinear contact and collision between a defunct satellite and a deformable net.Finally, a log barrier functional was developed to achieve the complex coupling interaction between a sliding ring and a deformable rod, to mimic the linkage used in the closing mechanism.Importantly, the C 2 continuous incremental energy for both nonpenetration contact and nondeviation link was differentiable and, therefore, the nonlinear elastodynamic net capture system could be updated robustly through an implicit Euler integrator.The whole space debris capture process using a hexagon net was simulated by our newly introduced numerical framework.In the current study, the closing mouth of the tether-net system is still relatively large; it would be interesting to numerically explore the optimal closing time, the closing speed, and the closing angle to minimize the closing mouth and increase the capture rate.In the future, the numerical setup can be combined with rigid-flexible coupling to study the detumbling behavior of a spinning satellite when twining with a flexible net.Also, it would be interesting to develop a model-based controller for tether tugging during the post-capture phase.][39]

F
for details.The total potentials are the sum of elastic stretching and bending throughout the net system, I G U R E 2 (A) Discrete stretching elements and (B) discrete bending element.HUANG ET AL.

1 , 2 , 1 ,
we use the PARDISO package to solve the sparse linear system.The overall simulation process can be divided into three phases: spread, contact, and close.The physical and geometric parameters used in our numerical investigation are as follows: net length L = 15.0 m, satellite wing width D = 1.0 m satellite wing length D = 4.0 m relative distance, H = 17.0 m, Young's modulus E = 1.0 GPa, radius of the net rod, r = 0.1 mm radius of the boundary rod, r = 0.5 mm 2 , material density ρ = 3000 kg∕m 3 , and bullet mass 1.0 kg.The grid size is selected as L Δ = 0.5 m and is refined at the inner area.The total degrees of freedom for our discrete system are N 3 = 28 371.The discrete time F I G U R E 4 (A) A rigid ring sliding on a flexible rod and (B) pointto-edge connection element.HUANG ET AL.

Finally, the flying
net needs to be closed to increase the success rate for space debris capture.The closing mechanism considered here is by stretching the boundary rod, and due to the ring connection between the boundary rods and the flexible net, the tether-net system can be closed.The connection parameters are as follows: barrier distance r ˜= 0.1 m and barrier stiffness K = 100 N∕m J .At time F I G U R E 5 Snapshots of a tether-net capture system during the spread phase (A-F).t = 3.6 s, the bullet mass would explode into two parts and fly back to back to achieve the close mechanism, and the shot velocity of each flying part is set as v = 10 m∕s c .Due to the velocity of the closing bullet mass, the boundary rod would be stretched out of the rings located at the corner, and, therefore, the open area of the tether-net system can be reduced to achieve the closing mechanism.Several snapshots for the closing process are shown in Figure 7A-F.