{'fraction_non_alphanumeric': 0.042400307169107566, 'fraction_numerical': 0.01439855191706434, 'mean_word_length': 4.558384146341464, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Dynamic Modeling of a Tether-net System for Space Debris Capture

In this paper, a flexible tether-net system is applied to capture the space debris and a numerical framework is established to explore its nonlinear dynamic behaviors, which comprises four principal phases: folding, spreading, contacting, and closing. Based on the discretization of the whole structure into multiple nodes and connected edges, elastic force vectors and associated Jacobian matrix are derived analytically to solve a series of equations of motion. With a fully implicit method applied to analyze the nonlinear dynamics of a slender rod network, the involved mechanical responses are investigated numerically accounting for the interactions. Contact between the deformable net and a rigid body is handled implicitly through a cost-effective modified mass algorithm while the catenary theory is utilized to guide the folding process (from planar configuration to origami-like pattern). The dragging and spreading actions for the folded hexagon net could be realized by shooting six corner mass toward a specific direction; next, the six corners would be controlled to move along a prescribed path producing a closing gesture, when touch between the flying net and the target body is detected, so that for the space debris could be captured and removed successfully. We think the established discrete model could provide a novel insight in the design of active debris removal (ADR) techniques, and promote further development of the model-based control of tether tugging systems.


Introduction
Space debris, according to European Space Agency (ESA), is the non-functional, human-made object in Earth orbit or re-entering into Earth's atmosphere, and more than 20, 000 such objects of larger than 100mm diameter have been detected in 2020 [1], which threatens the safety of other orbiters and motivates the demands for active debris removal (ADR) systems. Traditionally, two types of ADR frameworks are widely used: (i) contact-based ADR method, and (ii) non-contact ADR method [2]. Robotic arms and harpoons serve as two typical systems of the former type, for which the space debris is connected to a parent satellite. As for the non-contact ADR method, laser and ion beams are usually utilized while physical contact is unnecessary. Different from these technologies, the tether-net capture system emerges in recent years and gains increasing attentions. It is thought as one of the most promising techniques to relieve the debris threats due to several advantages: (i) a tether-net system is able to travel a long distance between the parent satellite and targeted objects; (ii) it is such a flexible system that it can adapt to different shapes of the target objects; (iii) a tether-net is usually lightweight and cost-efficient relatively [3]. Because of these advantages, the tether-net system has already been largely investigated and used in ADR [4][5] [6].
Plenty of investigations have been conducted on tether-net capturing method. A commercial software, Vortex Studio, was established by Botta et al. for the dynamic simulation of tether-net systems [7], including net deployment process [8], net contact mechanics [9], and tension-based net closing mechanism [10].
The lumped-parameter method was employed by Benvenuto et al. to study the entire capture process, especially facilitating the navigation and control of flexible nets [11] [12]. Golkebiowski et al. developed a Cossrat rod theory and considered the bending and twisting forces for one-dimensional slender body [13]. A comprehensive parameter sweep was implemented by Shan et al. to quantitatively analyze the effects of bullet mass, shooting velocity, and shooting angle during the net spread phase [3]. A line-line collision detection theory was proposed by Si et al. to accurately simulate the self-contact behaviors of a tether-net structure during the capture process [14]. Endo et al. studied the acceptable distance between the net and the debris for a successful capture, based on which an optimized deployment framework was developed [5]. The Floquet theory was used to analyze the nonlinear dynamic behaviors of tethered satelite system by Zhu et al. [15]. In parallel with numerical modeling, experiments were also performed to verify the effectiveness of tether-net capture systems [13][16] [17].
From a computational viewpoint, the modeling and simulating of a slender object, e.g. rods and their networks, are of sufficiently general interest, due to the preponderance of nonlinear deformations involved. Finite Element Method (FEM) has been the most commonly used framework in structural mechanics over the past decades [18]. Some other numerical methods, such as Geometrically Exact Beam Formulation [19], Absolute Nodal Coordinate Formulation [20], Isogemetric Collocation [21], and Mixed Isogemetric Finite Element Method [22], were also developed to reveal the mechanics of slender structures experiencing finite deflection and rotation. In this decade, discrete differential geometry (DDG)-based methods are becoming increasingly popular in the computer graphics community as well as mechanical engineering community [23] [24], for their computational efficiency and robustness in handling geometric nonlinearity, collision, and contact, which play important roles in the space debris capture and removal systems. The DDG approach first discretizes a smooth structure into a mass-spring-type system, while the key geometric characteristics are preserved, e.g., bending curvatures [23]. Previous DDG-based methods have shown surprisingly successful performance in simulating thin elastic structures, e.g. rods [25] [26], ribbons [27], and plates/shells [28] [29]. Tether-net or gridshell, on the other hand, usually represents a curved surface comprised of multiple one dimensional elastic rods and differs from the traditional one-dimensional rods or two-dimensional shells. It leaves opportunity for a new numerical model realizing efficient and accurate simulation of the tether-net system. In 2018, Baek et al. developed an explicit numerical framework based on Discrete Elastic Rods (DER) algorithm to investigate the buckling and rigidity of a grid of rods connected by joints [30] [31]. This numerical framework was later improved by Huang et al. to explore the nonlinear mechanics of rod networks [32] [33] [34] [35]. Hou et al. applied the DER method to constructing an energyconserving integrator for the purpose of modeling the dynamics of a tether-net capture system [36]. As the discrete curvature in original DER formulation might be singular and thus encountered numerical issue, the bending energy at the folding points was not taken into account.
In this work, a DER-based framework is established to simulate the whole dynamic capture process of space debris by a tether-net system, which is tackled in a numerical way so that the geometrically nonlinear deformations can be captured automatically. A modified energy functional is proposed to evaluate the nonlinear bending deformations of folded structures. Instead of applying Automatic Differentiation [37] that was employed by Hou et al. [36], the gradient vector (associated with elastic force vector) and the Hessian matrix (associated with Jacobian matrix) of total potentials are derived analytically, as a result of which a rather simple mathematical format is obtained. Moreover, a computationally efficient modified mass method is implemented to perform the interaction between the flexible net and the rigid debris, different from the penalty energy-typed formulation for contact mechanics in [36]. In order to avoid the numerical issues experienced in the folding process (from a planar net to an origami-like pattern), a modified bending curvature for slender rod-like objects is constructed, inspired by the discrete bending energy functional formulated in Föppl-von Kármán plate equations [38][39] [29]. The elastic force vectors and associated Jacobian matrix required are derived analytically based on DDG to solve the dynamic governing equations of a net system [23] [40]. Combination of the implicit Euler method and Newton-Raphson optimization facilitates the numerical integration of the motion equations step by step, leading to a robust performance and unconditional convergence [41]. Contact between the flexible net and the rigid target is handled by a computationally efficient framework -modified mass method [28][42] [32]. After a catenary theory is utilized to guide the transportation, a folded hexagon is adopted and controlled via internal mechanical links to realize the capture process consisting of folding, spreading, contacting and closing. In the end, five representative scenarios are tried to show the effectiveness of our numerical models: (i) nonlinear dynamic vibration, (ii) the contacting process, (iii) the folding process, (iv) the shooting process, and (v) the closing process. Excellent performances show the success of our proposed discrete net framework in simulating the whole debris capture process, which not only provides a novel insight in space debris cleaning, but lays a basis for the interactions between tether-net-systems and rigid bodies in other fields.
The novelties of the current study are as follows. 1. A new bending energy functional is provided, and its first and second variations are given analytically based on DDG. This new energy functional can not only capture the nonlinear bending deformations of rod-like objects, but can also handle the numerical singularity at the folding area. 2. Our contact model is a fully implicit approach, such that a relatively large time step size can be used to speed up the simulations. 3. Catenary theory is considered to guide the folding process, and the difference between the bending-dominated beam theory and the stretching-dominated cable model is discussed in detail, resulting a general numerical tool for all 1D rod-like structures.
The remainder of the paper is organized as follows. In Section 2, the nonlinear dynamic simulation for tether-net system is formulated comprehensively based on DDG theory. Next, in Section 3, some numerical results for net capture process are presented in details. Finally, conclusive remarks and future research directions are discussed in Section 4.

Numerical model
In this section, we introduce a numerical framework for the dynamic simulation of tether-net system through a geometrically exact formulation. Referring to Fig. 1(a), a hexagon net woven by several threads is considered as the capture covering. Side length of the hexagon is L = 10.0m and the mesh density is characterized by an grid interval of dL = 1.0m. Each rod in the hexagon net is assumed to be manufactured by an isotropic, linearly elastic material with Young's modulus E = 1.0GPa, material density ρ = 1000.0kg/m 3 , and a circular cross section of r 0 = 1.0mm diameter (and, therefore, cross-sectional area A = πr 2 0 and moment of inertia I = πr 4 0 /4).

Discrete model
The continuous structure is discretized into N nodes in order to numerically represent the geometrically nonlinear deformation under mechanical loads. As a result, a 3N sized degree of freedom (DOF) vector is produced, where denotes the position of i-th node. The edge vector between the i − th node and j-th node is expressed as and a unit vector in its tangential direction is

= 20m
= 1m x i x j x k x l e ij e jk e kl x i ( ) ) x + i ( ),-)

Elastic force vector
The total potential of an elastic net system includes two manifests: (i) elastic stretching energy, E s , and (ii) elastic bending energy, E b . Referring to Fig. 1(b), the linearc stretching energy associated with the elongation of (ij)-th edge is where ij is its uniaxial stretching strain, Here, a bar on the top represents the quantity in undeformed configuration, e.g., |ē ij | is the undeformed length of the (ij)-th edge. The bending energy between the (ij)-th edge and (jk)-th edge is given by their misalignment, where ∆l ijk = (|ē ij | + |ē jk |)/2 is its Voronoi length and κ ijk is an approximation of discrete curvature, It should be noted that the real curvature is [25][26] [36] However, to avoid the singularity and numerical issues induced during the folding process, Eq. (7) is applied to formulate the bending energy in a flexible rod-like structures. The discrete curvature in Eq. (7) is demonstrated being able to well approximate the original formulation in Eq. (8), even in the geometrically nonlinear regime, referring to Appendix A for more details. The total potentials are the sum of elastic stretching and bending throughout the net system, 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, F int ∈ R 3N , is the negative gradient of total energies, To formulate a mathematical expression of energy gradient, the operator ∇ is introduced, Based on DDG theory, the variation of (ij)-th elastic stretching energy with respect to (ij)-th edge is computed as, Similarly, variation of the (ijk)-th elastic bending energy is where Here, I 3×3 is a 3 × 3 identical matrix and ⊗ represents a tensor product. Finally, the gradient of elastic potentials with respect to DOF vector q (or x i with i ∈ [0, N − 1]) can be easily derived following the chain rule in Eq.(2).

Discrete equations of motion
The equations of motion for a dynamic net system is [41] where M ∈ R 3N ×3N is the diagonal mass matrix, F d = −µMq is the damping force (µ is damping coefficient andq is the velocity of DOF), F g = diag(M)g (g is gravitational acceleration vector and its z-component is g) is the external gravity force, and() represents a derivative with respect to time. The Euler method is used to update the DOF vector and its velocity from time step t k to t k+1 = t k + h (h is the time step size): The Newton-Raphson method is applied to solve this set of nonlinear equations of motion. Here, the implicit Euler can simplify the calculation by using the prediction-correction formulation, and the higher order time integration method can be found in Appendix B [41]. At the time step t k+1 , a new solution is first guessed on the basis of the previous state, i.e., Then, the solutions are optimized by the gradient decent algorithm, such that the new solution at the (m + 1)-th step is where J ∈ R 3N ×3N is the Jacobian matrix associated with Eq. (16), and all the quantities are evaluated at the m-th iteration when t = t k+1 . The equations of motion are updated iteratively until the error at the current time step decreases within the preset tolerance. The main difficulty for the formulation of Jacobian matrix appears in producing the second variation of elastic potentials, known as Hessian matrix, We here provide an explicit formulation of Hessian matrix of elastic potentials. The second variation of (ij)-th stretching energy is given by The second variation of (ijk)-th bending potential is computed as Finally, the Hessian with respect to DOF vector q (or x i with i ∈ [0, N − 1]) can be transformed via the chain rule in Eq.(2).

Contact model
The contact mechanics between flexible net and rigid debris will be addressed in this section, with the modified mass method is employed [28][42] [32]. As shown in Fig. 1(c), the target body can be described by z = f (x, y). A node in a discrete net structure, x i (t k ), would approach to the target surface at time t = t k .
If the rigid surface is not considered, the nodal positio at next time step would bex During the time marching scheme of dynamic simulation, when the node falls under the target surface, and in this case correction is required to movex i (t k+1 ) along a prescribed path onto the surface, where p i represents the surface normal vector, as shown in Fig. 1(c).
To ensure the non-penetration condition and associated displacements, the equations of motion at the i-th node in Eq. (16) have to be modified slightly. For the three degrees of freedom at the i-th node, where m i is the lumped mass, F int i ∈ R 3 (as well as F d i and F g i ) is the 3-element internal force (as well as damping force and gravitational force) on the i-th node, and the modified mass matrix is, when free DOF of i-th node = 3, where p i (similar for p 1,i and p 2,i ) is the constraint direction. As for a free node, ∆x pre i = 0, and Eq.(26) can be simplified to Eq. (16). If the node is fully constrained i (t k+1 ) = 0 and the change in position (as well as the velocity) is enforced to take the prescribed value. In this paper, the node is constrained along surface normal, p i , and, therefore, the free DOF number of i-th node is 2. Due to the flexibility of elastic net, inelastic collision between the i-th node and the target 3D surface in our numerical implementation, i.e., the velocity of constrained node is manually set to be zeros after contact. The contact detection is employed at the end of each time step: if any node falls under the target surface, the optional second solution is used as a corrector step. Pseudo code of the computational flowchart can be found in Algorithm 1.

Results
In this section, the numerical results are presented from our discrete net simulation. The nonlinear dynamics of a net under gravity is first introduced, followed by the contact simulation between a flexible net and rigid hemispherical surface. The folding, shooting, and closing of a tether-net are later explored numerically. The physical and geometric parameters are illustrated in Table 1.

Nonlinear dynamic vibration
In this subsection, the nonlinear dynamic vibration of an elastic hexagon net under gravity is simulated. In addition to the geometric and physical parameters Section 2, the gravitational acceleration is chosen to be g = −1000.0m/s 2 to address the geometrically nonlinear deformation. The negative gravity means its direction is along the minus z-axis. The exaggerated gravity can make the structure deform into a nonlinear phase, which is the novelty of our discrete model. The damping coefficient is selected from µ ∈ {0.01, 0.10, 1.00}. 5 nodes are interpolated between each mesh length dL as shown in Fig. 1(a), producing a total nodal number N = 3631, total stretching elements N s = 3960, and total bending elements N b = 3843. The numerical tolerance is set to be 10 −4 . It should be noted that all these parameters can be changed conveniently in our auto mesh generating algorithm. With six corner nodes (marked as red dots in Fig. 2(a)) fixed through a Dirichlet boundary condition, an initially planar hexagon net (Fig. 2(a)) would drop and deform due to the existence of gravity. Two representative deformed patterns are illustrated in Fig. 2  coefficient, e.g., µ < 0.1. Equilibrium diagram can be derived directly when the damping is relatively large. It should be noted that with damping forces added into the dynamic system, the equilibrium configurations of elastic structures can be eventually derived through the well-known dynamic relaxation method, which has been widely adopted to solve nonlinear mechanics problems. On the other side, the dynamic behaviors could be maintained once the damping is removed from dynamic system and symplectic time integration scheme (e.g., Newmark-Beta method) is employed for the integration of equations of motion [41]. The computation time is discussed in Appendix C.
The stress distribution during the equilibrium configuration can be found in Fig. 4. Here, the stress of (ij)-th stretching element is and the stress of (ijk)-th bending element is Because both stretching stress and bending stress are along the rod longitude direction, the overall stress is the sum of the stretching and bending. Moreover, our numerical model is more general and can deal with nets with different grid size, e.g., Fig. 5 shows the equilibrium configurations of nets with different grid interval, L/dL ∈ {10, 15, 20}. Different simulations can be performed with different input file, e.g., initial nodal positions, bending elements, stretching elements, and constrained index (or initial conditions).

The contacting process
Next, a numerical setup shown in Fig. 6(a) is discussed to demonstrate the effectiveness of our developed contact model. The physical parameters stay the same to the previous scenario (except for gravity, which is g = −10.0m/s 2 ), and the damping coefficient here is chosen to be µ = 0.1. The geometry of the target surface is given by Γ (x, y, z) : where R h = 4.0m is radius of the hemispherical shell, and the net locates at z = 5.0m, slightly above the target. Due to the existence of gravitational force, the net would drop and touch the target object; Fig. 6 shows several representative snapshots of the contact process between the hexagon net and a hemispherical body. The contact area gradually increases from the beginning, e.g., t < 1.6s, and a vibratory motion of the net edge can be observed after t = 2.0s, because of the structural inertia. An equilibrium configuration is obtained after a sufficient long period of time, e.g., t > 8.0s; Associated dynamic rendering of the contact process can be found in Supplementary Material Contact.mov.

The folding process
The extremely large net needs to be folded before the transportation to space. We here propose a strategy to fold an initially planar net into an origami-like pattern. Due to its flexibility, it is suitable to guide the folding process with the catenary theory. All physical and geometric parameters are identical to the previous contact case. Referring to Fig. 7(a1)-(a4), we plot the topology of a suspended cable with a shrinking boundary condition; here, the cable length is L, compression is ∆L, and the maximum displacement is denoted by δ. If the two boundary points, {x 0 , x N −1 }, are fixed only and the structure is shrunk from L to L − ∆L, the maximum displacement of a suspended cable on the basis of catenary theory is, where the coefficient A c can be derived by solving the following transcendental equation, Norm compress, Δ / Norm deflection, On the other side, when more nodes are constrained and manually shrunk along a prescribed trajectory, referring to Fig. 7(a2)-(a4), the effective rod length L and the corresponding maximum deflection δ would linearly decrease as the enlargement of suspension points, N F , e.g., the maximum deflection for Fig. 7(a2) is only half of the deflection observed in Fig. 7(a1). Fig. 7(b) shows the dependence of normalized maximum deflection, δ/L, on normalized compression, ∆L/L, with different number of suspension points, N F ∈ {1, 2, 3, 4}. Here, the solid lines are derived from analytical solution in Eq. (31) and Eq. (32), and the markers are from our discrete simulation. Excellent agreements can be observed between the analytical solutions and numerical data. With a given size of net bag, the constrained nodal number and compressive distance can be selected based on the catenary theory. It is worth to mention that the bending curvature formulated in Eq.(8) would not be physically accurate and may induce numerical issue at the folding points, referring to Fig. 7(a2) for details. Difference between the modified bending curvature in Eq.(7) and the original bending curvature in Eq.(8) is detailed in Appendix A, and the difference between the stretching-dominated case and the bending-dominated case (which is essential for the folding phase) is discussed in Appendix D.
The folding process of a planar hexagon net is shown in Fig. 8. It should be noted that the curve between the suspended points (red dots) is a catenary. Here, 331 nodes (marked as red dots in Fig. 8(a)) are constrained and manually moved towards a specified position, such that the planar net could be folded into a small space with a smaller hexagon with the edge length being 1.0cm. Representative shapes for different compression distances are given in Fig. 8(b)-(f), separately. Associated dynamic rendering can be found in Supplementary Material Fold.mov. The constrained nodal number and compression distance can be tuned to adapt to different sizes of the net bag.

The shooting process
When a folded net is shipped to space, six corner mass would be shot out of the net bag such that its deployable form could be spread and enter the space to capture target debris. In this process, the gravity is no longer considered, and all other physical and geometric parameters stays the same with the previous folding case, Each corner mass is 5.0kg and their initial shooting speed is 20.0m/s. The shooting angle, θ, should be selected between the speed directional vector and z-axis. The initially configuration, q(0), is the fully folded state in Fig. 8(f) and the initial speed,q(0), is zero, except for six corner mass. The planar state of elastic net can be achieved and maintained by shooting and controlling the speeds and locations of six corners.  Fig. 9(a) shows the spreading area as a function of time, for different damping coefficients, µ ∈ {0.00, 0.01, 0.10}. The shooting angle is fixed as θ = 45 • . It is found that the effect of damping force is almost negligible as long as µ < 0.1. Dependence of the spreading area over time for different shooting angle, θ ∈ {30 • , 45 • , 60 • } is given in Fig. 9(b) where the damping is selected as µ = 0.1 to avoid unrealistic vibration. It reveals that the larger the shooting angle is, the faster the net spreads, which is consistent with the numerical data in [36]. Several representative renderings for θ = 45 • , µ = 0.1 are provided in both Fig. 10 and Supplementary Material Shot.mov. Actually, the bending energy was not included until the associated nodal position is above the net bag, e.g., z i > 0.0, which means the relevant bending element has already been dragged out of net bag and the structure could spread without any constraints. Otherwise, the net would expand without shooting.

The closing process
Finally, the closing mechanism of a tether-net system in post-contact phase is discussed. The numerical setup can be found in Fig. 11. All physical and geometric parameters are consistent with the previous shooting scenario, specially µ = 0.1 and θ = 45 • . Initial distance between the net bag and target spherical shell (with radius R h = 3.0m) is S H = 10.0m.
At the time step t ≈ 0.75s, the shooting net would touch the target, and high tension configuration appears as a result of the competition between contact and inertia, referring to Fig. 12(b). After t = 0.9s, the six net corners would be manually controlled and fly along a prescribed path until reach to the destination, [0.0, 0.0, 4.0]m, such that the target body could be packaged and tugged to an ideal environment. Speed of the six corner mass during the closing process is chosen as 10.0m/s. Several snapshots are provided in Fig. 12, and the associated dynamic renderings can be found in Supplementary Material Close.mov.

Conclusion
A DDG-based numerical simulation has been developed to reveal the nonlinear dynamics of a tethernet system during the capture of space debris. The whole process is realized with a tether-net system comprising four manifests: (i) folding, (ii) shooting, (iii) contacting, and (iv) closing. A continuous net structure is discretized into a mass-spring-damper system, with a lumped mass located as each vertex and the associated discrete stretching and bending energies. The nonlinear bending energy is modified to avoid numerical singularity during the folding phase. The unilateral contact between the deformable net and the rigid body is evaluated implicitly through a modified mass method, and the deployment process is guided by the classical catenary theory. Finally, several representative examples during the capture process is performed by the discrete algorithm, which validates the effectiveness of our proposed numerical framework. In the future, a systematic parameter sweep can be utilized to further optimize the design of a tether-net system that is more computationally efficient. It would also be of interest to develop a modelbased control theory for the tether tugging system to improve its performance during the post-capture phase of space debris. We hope our method can also be used for the modeling and design of other tether systems [43], e.g., tether satellite [44] [45], tether antennas [46], and tether space robot [47][48] [49].

A Nonlinear bending energy
In this appendix, we discuss the difference between the original bending curvature derived in Eq.(7) and the modified bending curvature formulated in Eq. (8). The discrete bending curvature in Eq.(8) would be singular and may experience numerical issue once the turning angle, φ, between two consecutive edges is relatively large, which is the case during the folding process, referring to Fig. 8(a2). Here, we demonstrate that the adapted discrete curvature given in Eq. (7) could not only get rid of this numerical discontinuity, but also still be accurate enough.
First of all, the variation of normalized curvatures, {κ, κ}, on the turning angle, φ, between two tangential directions are detailed in Fig. 13(a2). Quantitative agreement has been found when the turning angle is small, e.g., φ < 40 • ; however, κ would grow much faster than κ when φ > 60 • . This error could be eliminated with a dense mesh, because the turning angle φ could decrease as the increasing of discrete nodal number with a specified continuous curve. Notice that it would not be physically correct to use discrete curvature,κ, to formulate a 1D curve with C 1 continuity, referring to Fig. 7(a2). Also,κ = 2 tan(φ/2) and would be infinite when φ is close to 180 • , which is the case when a net is folded completely.
The accuracy of the bending curvature proposed in Eq. (7) is rendered by using a nonlinear cantilever beam, as shown in Fig. 14(a). For a clamped-free beam with length L, density ρ, cross section area A, Young's modulus E, moment of inertia I = πr 4 0 /4, and gravity g, its tip displacement in geometrically linear phase can be obtained classical Euler-Bernoulli beam theory, In this numerical setup, we choose L = 1m, ρ = 1000kg/m 3 , E = 1GPa, r 0 = 1mm, and g is varied for different loading conditions. The structure is discretized into N = 50 nodes, and the first two nodes, {x 0 , x 1 }, are manually fixed to achieve a clamped-like boundary condition. In Fig. 14(b), we plot the relations between the beam tip deflection and its own weight from two different bending models and linear analytical solution. Both numerical results match well with analytical solution provided in Eq.(33) when the normalized weight is moderate; also, the variation between the curvature formulated in Eq. (7) and in Eq. (8) are negligible, and the overlapped markers indicate the accuracy of our simplified bending formulation for folded tether-net system, even in a geometrically nonlinear range.

B Time marching scheme
In this appendix, we discuss the first order Euler method and the second order Newmark-beta method for the time integration. It is known that the first order Euler method would experience artificial damping when the time step size h is relatively large, and this issue can be overcome if the symplectic Newmark-beta method is used [41].
The Euler method is: The Newmark-beta method is: Here, the nonlinear vibration of a cantilever beam under gravity is considered. In this numerical setup, we choose L = 0.1m, ρ = 1000kg/m 3 , E = 10MPa, r 0 = 1mm, and g = 10.0m/s 2 . Time step size h is varied to show the accuracy of dynamic simulation. In Fig. 15, we plot the beam tip deflection as a function of time for both Euler method and Newmarkbeta method. When the environmental damping is ignored, artificial damping is observed when the time-step size h in Euler method is relatively large. However, the numerical error can be overcome if the symplectic Newmark-beta method is employed (here β = 0.5). On the other side, if the environmental damping is taken into account, the artificial damping from Euler integration would be trivial and the solution is acceptable. For the flexible system, the environmental damping and material damping are usually relatively large, such that it would not experience a long time vibratory motion, and, therefore, the naive Euler method is largely used. If the energy conservation is concerned, the time marching scheme can be easily switched based on our previous work [41].  |   Fig. 14 (a) Undeformed and deformed patterns of a cantilever under its own weight. (b) variation of normalized tip displacement, δ/L, on normalized weight, ρgAL 3 /EI, from different models. environment by Eigen [50]. We can see that real time simulation can be achieved when N ≤ 1000 if h = 0.01s. It should be noticed that the computational time would almost double if the contact is considered. PARDISO package is utilized when solving a sparse linear system formulated in Eq. (18) [51][52] [53]. The simulations are performed on a single thread of Intel Core i7 − 6600U Processor @3.4GHz. In the future, parallel computing could be implemented to speed up the solver.

D Comparison between beam model and cable theory
In this appendix, we discuss the difference between the bending-dominated beam model and the stretching-dominated cable theory, which is essential for the folding process. Beam is a structure that experiencing bending, while cable is an object that undergoing pure stretching. For a 1D rod-like structure with length L, stretching stiffness EA, bending stiffness EI, and experiencing external load F (or ρALg for gravity), the deformed mode is based on the ratio among EA, EI/L 2 , and F , i.e., the beam model would be effective as long as EA F ∼ EI/L 2 , and the cable model should be considered in the phase of EA ∼ F EI/L 2 . We here consider a 1D object with pin-pin boundary condition, referring to Fig. 17(a). The numerical parameters are: structure length L = 1.0m, cross section radius r 0 = 1cm, Young's modulus is E = 1GPa, and number of vertex N = 100. If the magnitude of midpoint load F is similar to the normalized bending stiffness, the deformed solution should be a curved line and can be predicted by the classical Euler-Bernoulli beam theory On the other side, if the external force is quite large, the deformed pattern should be a zig-zag line and can be computed based on equilibrium condition, F EA = 2 δ 2 + (L/2) 2 − L/2 L/2 δ δ 2 + (L/2) 2 .
Our rod-based simulation can predict both bending-dominated phase and stretching-dominated phase, referring to Fig. 17(b). For a 1D structure with compressive load, e.g., folding process, the predictions from beam model and cable model would be different: if its bending stiffness is large and is similar to the external gravity, its deformation can be obtained based on Euler buckling theory, known as Euler Elastica; however, if the bending stiffness is quite small and the stretching stiffness is similar to external gravity, it becomes a classical catenary model. The deformed patterns for bending-dominated modes and stretching-dominated modes are plotted in Fig. 18. Different modes are also related to the ratio among EA, EI/L 2 , and ρALg.
Usually, the net is soft and its bending rigidity would be small compared with EA and ρALg during the folding phase (which is on earth), such that the bending effect (e.g., discontinuous point) would be trivial to be included. However, when the gravity is missing (e.g., on-orbit phase), the bending would be more important and cannot be ignored any longer. Moreover, the folding size is also determined based on the ratio among EA, EI/L 2 , and ρALg. For example, if the distance between two folding points is large (which means EI/L 2 is small), the bending energy can be ignored (e.g., dashed red line in Fig. 18(c)); however, as the decreasing of folding distance L, EI/L 2 would gradually become large and the deformed pattern would switch from stretching to bending (e.g., solid green line in Fig. 18(c)), and the discontinuous bending would have large energy, such that catenary-based folding process would fail. Overall, the effective folding distance for catenary theory is determined by the ratio among EA, EI/L 2 , and ρALg.

E Video
We provide several videos corresponding to Fig. 2

Availability of Code
The code generated during the current study are available from the corresponding author on reasonable request.