Modeling of a Tethered Testbed for a VTVL Vehicle

This paper proposes an algorithm for modeling a three-dimensional tethered environment for testing vertical-take off, vertical landing vehicles. The method is able to take several geometrical configurations into account and combines the classical catenary model with the elasticity theory to predict the forces acting on the lander in quasistatic conditions, i.e., in conditions of hovering, where the motion of the vehicle is reduced. Numerical results confirm that the method is potentially able to provide real-time solutions, which can be included as feedforward contributions in the design of tethered experiments.


Introduction
Recent and future missions involve a precise descent and landing in addition to the ascent phase to reach the target orbit. This can be on the one hand the powered descent and landing of a reusable first stage of a launch vehicle as it was demonstrated several times by SpaceX with its Falcon 9 launch system. On the other hand, precise descent and landing has been applied to planetary missions and is foreseen for many more future missions to Mars and to the Moon. The development of guidance, navigation, and control (GNC) techniques for these applications remains a challenging task although several missions have been already successfully completed. Several ideas to support and accelerate the GNC development using demonstrators have been conceived in the past, for example NASA's Morpheus lander [1,2] or the HOMER demonstrator of Airbus Defense and Space [3,4]. In these and other developments of space systems, a wide variety of tethered experimental setups has been created [1,2,[5][6][7][8].
For vertical-take off, vertical landing (VTVL) vehicles, the use of tethered solutions has pros and cons. One of the main drawbacks is the impact of the tethers on the flight characteristics of the vehicle. They introduce effects which do not exist in the final free-flight scenario. A further downside is that the allowed flight envelope is usually quite small due to limited tether length. Nevertheless, it is sufficient to test the hovering capability, one of the first milestones to be achieved towards the development of the full free-flight capability. On the other side, a tethered configuration provides a safe environment for the vehicle. It is mitigating the effect of failures until the technology under development is mature enough to allow a reliable free flight. Moreover, the use of tethered configurations allows for making the test facility a protective area for the team of engineers and researchers. Thus, they can safely and closely track the progress in the development of the vehicle. This paper addresses the problem of modeling the tethered testbed during the hovering experiments of the VTVL vehicle EAGLE (Environment for Autonomous GNC Landing Experiments) developed by the German Aerospace Center (DLR) [9][10][11]. The modeling focuses on the tethering effects on the vehicle in hovering conditions. It can potentially be extended to any vehicle where it is needed to include the forces generated by the tethers in the design of the controller.
The modeling starts with the use of the catenary, a wellknown concept in the field of structures and mechanics [12], which physical meaning is the ideal shape that a hanging rope or chain has when subject to its own weight while having its endpoints constrained in two points in the space. The catenary is slightly different from the parabolic profile, which is the shape that intuition would (erroneously) suggest. The catenary concept is widely used and can be extended to multibody net structures [13]. For instance, it can be used for the modeling of railway overheads or power lines as in [14] or [15]. In these applications, however, the shape is the most important aspect and is therefore the main factor to be studied. As a consequence, the forces at the suspension points are not so emphasized.
In [16], the catenary curve and the corresponding equilibrium of forces are discussed in the frame of the development of a controller for an unmanned vehicle. However, this work focuses on the design of a controller which minimizes the tension in the ropes. Richardson focused instead on the entanglement detection of swarms of robot in urban environments [17].
This paper focuses on the modeling and the analysis of complex tethered configurations including hanging ropes attached to the vehicle while hovering. In addition, elasticity is included in the model. The purpose is to know the forces which act on the VTVL vehicle and can be included as feed-forward contribution in the design of the EAGLE controller. We propose an algorithm implementing the aforementioned theories of catenary and elasticity that computes the forces acting on EAGLE in an iterative way.
The paper is structured as follows. In Section 2.1, the modeling of a hanging rope and the related concept of catenary are introduced. The corresponding forces at the suspension points with and without elasticity are computed. Numerical validations of the proposed modeling are carried out in Section 2.2. In Section 3.1, the specific geometrical test setup for EAGLE is described. The setup includes different rope materials and a configuration of three ropes which are linked to the VTVL vehicle in a triangular configuration.
In Section 3.2, the iterative solving procedure is illustrated. It is based on a MATLAB implementation, but it can be transferred to any other software. Section 4 shows some numerical results obtained with the proposed algorithm. In general, every single point within the test area can be tested. However, for a better characterization of the scenario involving the EAGLE motion, two specific hovering paths are analyzed, and the accuracy of the computed solutions is discussed. Finally, in Section 5, we draw some conclusions on the work done.

Rope Modeling
In this section, a brief introduction about the modeling of a static rope is given, and the computation procedure of the forces at the suspension points is explained.
The function describing the shape of a hanging rope subject to a constant gravity is called catenary [18]. The catenary curve is defined as where r is the radius in the vertex, x 0 is the x-coordinate of the vertex and y 0 is the y-offset of the vertex. A two-dimensional x-y reference frame is used for the description of the motion. The gravity is directed towards -y, with a horizontal x-coordinate and a vertical y-coordinate. An example is depicted in Figure 1.
We can derive three equations for three unknown parameters r, x 0 , and y 0 . We know the position of the two suspension points, defined as a = ðx a , y a Þ and b = ðx b , y b Þ. Moreover, the length of the rope l is known.
With this equations, we get the following system: For a simplification of the system of equations, the catenary curve is shifted. The suspension points are moved, such that the first suspension point is on the origin. Thus, a local coordinate system positioned at the first suspension point is used. Now, Equations (4) and (5) can be used to obtain a formula for x 0 : Here, y 0 is eliminated. Therefore, Equations (5) and (6) become a system of two equations for two unknowns x 0 We can rewrite the condition on the length of the rope represented by Equation (6), which gives This is a nonlinear equation, only depending on r, which cannot be solved analytically. At this point, a numerical method is needed.
Different methods can be employed. Classical Newton methods for this specific problem can experience numerical issues. This happens because the slope of the function rapidly tends to infinity. A better alternative is the Halley-method [19], which requires the second derivative to be continuous and the first derivative at the root different from 0.
The Halley method works with the following iterative scheme: Therefore, the derivatives of the functions f ′ ðrÞ and f ″ ðrÞ are needed. They can be computed as A reasonable initial guess for starting the algorithm is required. In this case, we have to look at the double dependency of the radius r on both the length l and the distance between the x components of the suspension points. The more the points are apart each other, the greater is the radius r. On the contrary, the greater l, the smaller r. Therefore, the initial guess ðx b − x a Þl is used. This causes problems because the curve radius r tends to 0. To solve this case numerically, the points are slightly shifted apart each other. The first suspension point is shifted 1 mm to the left and the second suspension point is shifted 1 mm to the right. Now, the x -distance between this two points is not equal to 0 anymore and a radius r can be calculated. This approach results in a small error in the radius. We tested it for suspension points which are above each other and with a distance of one to nine meters and a length of the rope which is two meters longer than the distance. These are the dimensions of the setup. The radius is in the range of 3 · 10 −6 instead of 0. Therefore, the error is small enough. Note that the case x a = x b is of no practical interest anyway and is analyzed for the sake of completeness of formulation.
To get the forces acting on the suspension points, the force equilibrium and the knowledge about the orientation of the forces are used. An example of configuration where we can benefit from this knowledge is depicted in Figure 2, and its exploitation leads to the following set of equations: where F G is the force caused by the weight of the rope. This is a system of four equations for four unknowns that are F a x , F a y , F b x , and F b y . To solve it, we can use the catenary function.
and the system can be rewritten as follows: So far, we considered two-dimensional representations of the problem. The general three-dimensional problem can always be reduced to the two-dimensional version as pictured Figure 2: Balance of forces acting on the rope. F a and F b are the forces at the suspension points. They act tangentially along the rope. F G is the weight force which acts in the negative y-direction. 3 Modelling and Simulation in Engineering in Figure 3. First, the first suspension point a = ðx a , y a , z a Þ T is shifted onto the origin 0, 0, 0 Therefore, the whole catenary is shifted by −x a ,−y a ,−z a ð Þ T . We perform then a rotation about the z -axis by the angle α, which is the angle between the plane containing the catenary and the yz-plane.
where the angle α is defined as After this transformation, all the points of the catenary lie in the yz-plane. In this new representation, the curve radius in the vertex r, the displacement of the vertex in z -direction (z 0 ), and the z component of the forces are already correct.
Once the forces are computed according to Equations (16)- (19), we can back-transform the catenary and the forces to obtain the three-dimensional solution to our problem. This is done with the transpose of the rotation matrix R z ðαÞ , which corresponds to R z ð−αÞ.
After the rotation, we have to back transform in terms of translation. This means that also the vertex of the catenary needs to be shifted.
The catenary is now described by The back-transformation for the forces is the same as with the displacements of the vertex.
In the following, some simple implementation results of forces at the suspension points of a simple hanging rope without elasticity are illustrated. For a 2.90 m rope with a mass of m = 1 kg/m, weight of F G = 28:449 N and suspension points a : ð0, 4Þ and b : ð2, 5Þ (as represented in Figure 4(a)), the forces at the suspension points are F a = −6:9996 It can be seen that the x-and the y components of the forces at the suspension points and the weight force cancel each other out. In addition, the magnitude of the force at the upper suspension point b is greater than that in point a.
If the left suspension point is changed to a : ð0, 3Þ as in Figure 4(b), the vertex is not between the two suspension points anymore and the forces change In this case, we observe that the sign of the y component in b changes, and both magnitudes are greater than in the previous case. This is caused by the changed slope of the catenary at the suspension points.
In both cases, the length of the rope is longer than the distance between the suspension points. But if the length is shorter, the rope has to be stretched. Thus, it is necessary to   Modelling and Simulation in Engineering consider the elasticity. Under the assumption that the rope is linear-elastic, the behavior of the rope can be modeled by using Hooke's law: where Δl is the change in the length of the rope, which is proportional through the spring stiffness k, to the force F. There are two cases which has to be described separately: (1) The length of the rope is at least the distance between the two suspension points (2) The length of the rope is shorter than the distance between the two suspension points In the first case, the rope has to be stretched by its own weight. To solve this case, we implement the following procedure: first the forces at the suspension points are computed for the non-stretched rope. These forces are used to calculate the corresponding stretch with Hooke's law through Equation (33). The stretch is then added to the length and the forces at the suspension points are calculated again. This iterative procedure is used until the variation of stretch in the rope reduces to a given threshold.
For the second case, the superposition principle is used. A massless rope, stretched exactly to the length of the distance between the suspension points, is overlaid with a rope which is stretched by its own weight. The force, which is needed to stretch the rope until it is as long as the difference between the suspension points, is calculated with Hooke's law through Equation (33). Its gradient is equal to the gradient of the linear connection between the suspension points. It is added to the force which results from the stretch due to the rope's own weight. These forces fulfill the equality system, which is described in Equations (11)- (14). The forces caused by weight stretch fulfill the equality system directly as explained before, and the forces caused by the stretch to the minimal possible length cancel each other out with respect to the sum of the horizontal forces and the sum of the vertical forces. Indeed, all these forces together fulfill the equilibrium of forces.
As example consider the case represented by Equation (25) and a stiffness of k = 1000 N/m, the stretch is dl = 1:1943 cm. The forces are now: The decrease in the horizontal components of the forces shows that the rope sags more than in the previous case. The vertical component of the force at the lower suspension point is larger, while the same component on the higher suspension point becomes smaller. In other words, the elasticity of the rope tends to slightly reduce the force gap at the suspension nodes. In the above second case as represented in Equation (29) with the same stiffness, the stretch is dl = 1:9263 cm.
The more important case is the one where the distance between the suspension points is larger than the length of the rope. For instance, if the length is changed to 2 m, the weight is F G = 19:62 N, and stretch becomes dl = 82:943 cm. The forces at the suspension points are in this case: F b j j= 1282:5037 N: As expected, these forces are large compared to the forces in the first two cases. This is because the stretch generates a significant additional force, which is significant larger than the forces discussed in (25), (29), and (34). In the application 5 Modelling and Simulation in Engineering of the proposed method to the VTVL note, we always have situations corresponding to case (1) described above, as we are interested to have small forces which do not modify the flight of the vehicle significantly.

Validation.
In this section, we want to check our results to see whether they fulfill the forces equilibrium, described by Equations (18) and (19), the length of the rope in Equation (5), and whether the suspension points are part of the esulting catenary function.
In the above first case, Equation (25) Thus, the catenary function with the given parameters connects the suspension points and the rope is as long as expected.
The same can be done for the above second case (29). As before the x components ±18.0433 N cancel out each other and the y components Thus, the catenary function with the given parameters connects the suspension points and the rope is as long as expected.
In the above third case (34) with elastic stretch by the ropes' own weight, the x components ±6.9157 N cancel each Thus, the catenary function with the given parameters connects the suspension points. As expected, the rope is longer than the initial length caused by stretch.
Finally for the above fourth case (38) including the stretch, which results from a too short rope, the Thus, the catenary function with the given parameters connects the suspension points. As expected, the rope is significantly longer than the initial length caused by stretch. It is even a bit longer than the direct connection between the suspension points, which is 2 ffiffi ffi 2 p = 2:8284. Therefore, the rope sags slightly and the parameters of the catenary can be specified.

Modeling Procedure for a Specific Experimental Facility
First, the given experimental set-up is introduced followed by the explanation of the modeling procedure for this set-up.
3.1. NEST. Tethered configurations have been used for testing several vehicles and spacecraft in a secured flight environment. In this section, DLR's experimental facility NEST (Nest Environment for Suspended flight Tests) for the vehicle EAGLE, depicted in Figure 5 is used as a test case for the proposed algorithm. The legs of the VTVL vehicle EAGLE are connected to three ropes, one for each of its three legs. These ropes are led over a turning wheel at the top of a traverse pole with height of 4:50 m through its middle and are fixed on the ground. Therefore, the configuration is now more complex as it includes a rope, and a change of direction at a turning wheel, assumed here to be frictionless. For each of the three ropes, one suspension point is the upper end of the leg of the VTVL vehicle and the other suspension point is the bottom of the traverse pole. Thus, the force balance cannot be applied as easy as explained in Section 2.1. The experimental setup is sketched in Figure 6.
The coordinate system has its origin at the base plate in the middle of the traverse pole. The z-coordinate point vertical. Therefore, the coordinates of the turning wheel are ð0, 0, 4:5Þ. The x-and y-coordinates represent the in-plane components, with x pointing towards the starting point of the lander (middle of the whole traverse system) and the y-coordinate appropriate normal, as showed in Figure 7. The coordinate of the middle of the whole traverse system on the ground is ð3:46,0, 0Þ, where L ground is for the NEST facility equal to 3:46 m.
We will limit the analysis to only one rope which is fixed at the lander and led over a turning wheel through a traverse pole. The procedure can be repeated for the other two ropes. Figure 8 shows the elements to be considered. The main part is a very stiff rope which ensures that the VTVL vehicle is restricted to stay in the allowed area. The last part is an expander rope. It is 1:50 m long and is less stiff. This part ensures that the VTVL vehicle is not subject to hard jerks when being captured by the tethers. Now that the scenario has been defined, we can see how to solve the problem. The proposed iterative method is the subject of the next section.

3.2.
Description of the Modeling Procedure. In this part, the solving procedure is explained. First, the cutting clear     technique is described and applied to the NEST facility. Two resulting systems are calculated individually considering the elasticity of the ropes.
It is assumed that the elastic behavior is linear, and Hooke's law can be applied. In addition, it is assumed that the expander rope does not get stretched above the turning wheel, which is geometrically represented only by a point (In the hypothetical case of the lander position outside of the triangular flight area the ropes are led over the turning wheels anyway, but this case is clearly excluded from this analysis, as we assume that the controller is able to keep EAGLE within the prescribed area).
A further assumption is that the rope is treated as elastostatic. This means that the catenary curve is not depending on time, and only the hanging of the rope with its forces at the suspension points is modeled. This assumption is justified by the small motion characterizing EAGLE during its hovering. Finally, as said, the rope is frictionless at the turning wheel.
To model the forces at the rope as it is seen at Figure 8, it is necessary to use the free-body principle and intersect the rope at the turning wheel. So, the system is divided into a traverse system and a lander system. This procedure is sketched in Figure 9.
F a is the force at the lander. F c is the force at the point where the rope is fixed on the bottom. F G l , F G t , and F G e are the weights of the stiff rope on the lander system and the traverse system and the weight of the expander rope in the traverse system. F t and F b are virtual forces which arise by reason of using the free-body principle. If the virtual forces are equal the whole system will be balanced.
If the lengths of the ropes in the single systems are known, it will be possible to calculate the forces in the single systems. The problem is to choose the rope lengths in the single systems such that the forces F t and F b are equal. The length of the expander rope and the total length of the stiff rope are known. The following algorithm make sure, that the stiff rope is split correctly into the two systems.
At the beginning, the stiff rope is divided into a landerside part and a traverse-side part. A specific offset is chosen as a starting shift.
Then, the forces in the single systems F b and F t are calculated and compared. If F b and F t are equal, the rope is split correctly and the algorithm is finished (exit condition). Otherwise, the split is changed such that the lander or the traverse system gets a longer section of rope depending on which of both forces is bigger. After that, the forces are calculated and compared and the split is changed again until the forces are equal.
To calculate the single systems with given rope lengths, it is necessary to consider the elasticity as described in Section 2.1. With those methods, we can treat the lander system.
(c) Free-body diagram lander system Figure 9: Free-body principle for a rope in NEST experimental setup. Intersect the rope at the turning wheel leads to virtual forces F t (at traverse) and F b (at virtual suspension point b) which have to be equal for satisfying equilibrium of forces at the whole system. The weight force F G is divided into its three parts, F G t (traverse side), F G l (lander side), and F G e (expander part).
Choose a specific starting shift; Calculate the forces Fb and Ft in the single systems; Change the split such that lander system gets a longer section of rope; else Change the split such that traverse system gets a longer section of rope; end end Algorithm 1: Algorithm to choose the correct shift at the roll. 8 Modelling and Simulation in Engineering The calculation of the traverse system is a special case and is discussed below. The traverse system consists of the expander rope and a part of the stiff rope. First, it is necessary to calculate the stretch in the rope due to its own weight. The formula for this stretch can be found by using After calculating the stretch, we can distinguish three different cases for the length of the rope in the traverse system: (1) The length is shorter than the height of the traverse (2) The length is equal to the height of the traverse (3) The length is longer than the height of the traverse In the first case, the rope has to be further stretched. The required force can be calculated with Hooke's law described in Equation (33) very easily. In the second case, the rope does not have to be more stretched. It has already the correct length. In the third case, a part of the rope is lying on the bottom. The rope does not have to be further stretched but it is necessary to consider that the part of the rope, which is lying on the bottom, does not stretch the rest of the rope by its own weight. The force F t on the top of the traverse is the sum of the force for the stretch caused by its own weight as given in Equation (38), and the force which is needed to stretch the rope even more.

Simulation Results
In this section, more complex examples than in Section 2.1 are shown. Moreover, the given experimental setup and the elasticity are taken into account. In general, every position in the allowed flight area can be analyzed but in this context two possible motions of EAGLE are considered as test cases. All the results have been computed with a desktop having the following specifications: Operating system: Windows 7 Enterprise 64 bits Processor: AMD FX(tm)-6300 Six-Core Processor; 3.50 GHz To simulate the proposed approach, the VTVL vehicle is assumed to be a point mass and is attached to three ropes. In the first example, we use a vertical trajectory resembling a lift-off followed by a landing maneuver. EAGLE is placed in the barycenter of the equilateral triangle formed by the supports. In this case, we model a 7 m-long rope, with the expander rope that has a length of 1.5 m when not stretched. The scenario is depicted in Figure 10. The vehicle moves upwards until a vertical displacement of 5 m is reached, from the initial height of 2 m to 7 m. For this scenario, this position corresponds to the maximum force exerted on EAGLE, which is constantly vertical due to the symmetry of the chosen configuration.
In Figure 11, the associated forces are plotted. Note that because of the symmetry of the scenario the norm of the three components are equal to each other and the horizontal component of each force is rotated by 60 deg with respect to each other. The maximum stretch is in the order of 27 cm, corresponding to a total force of about 180 N.
The second plot of Figure 11 compares the two fictitious forces. No peaks are observed and the virtual forces are equal; thus, the equilibrium of forces is fulfilled, confirming the validity of the numerical approach in computing the elastostatic equilibrium. The maximum residual observed is in the order of 9:5 · 10 −4 .

Example 2 Helix
Trajectory. The second trajectory is a helix-shaped profile. EAGLE changes its altitude slowly, while constantly staying within a given distance from the center of the NEST. This scenario is represented in Figure 12, where a three-dimensional view and three twodimensional views of the setup are depicted. Since this scenario shows a more complex behavior, to reduce the maximum force exerted on EAGLE, a length of 10 m was considered for the rope, with the expander rope again having an unstretched length of 1.5 m. The test trajectory of the VTVL vehicle is marked in yellow. The direction of the force caused by the tethering configuration is marked in red. The  magnitude is not in scale for visualization purposes. Only the direction of the force that the VTVL vehicle must compensate for to keep the position is plotted. During the trajectory simulation, a maximum length of 1.8 m is achieved, with a maximum stretch of about 30 cm observed in the expander rope attached to the third support. In Figure 13, the corresponding plot of forces can be seen. Note that at the beginning there is a big contribution coming from support 3 which exerts the larger force, while the other two contributions are basically equal to 0. When the lander moves towards support 3 the corresponding contribution decreases and the forces exerted by ropes 1 and 2 grow. Since the helix shows an harmonic motion in the plane X-Y the maximum force oscillates between rope 1 and 2 accordingly. Note that the residuals in the lower part of the figure always show a very good convergence of the algorithm.
Furthermore, there is an error plot in Figure 13, defined as the difference between the two fictitious forces at the turning wheel F t and F b . Also, in this case, the force residuals are very small, with a maximum value of 9:8 · 10 −4 . This indicates that also for a more complex motion of EAGLE the algorithm shows nice convergence properties.

Conclusions
In this paper, the modeling of a tethered configuration for testing a VTVL vehicle is described. After a short introduction of the catenary curve and the calculation of the forces at the suspension points of a single hanging rope, a specific geometry (NEST) including different values of length and rope stiffness is considered.
The model computes the forces of each rope at the lander and at the bottom of the traverse, where the ropes are fixed. The use of equilibrium of forces, together with Hooke's law provide a viable way to compute the forces acting on the vehicle in conditions of hovering flight.
The examples show that in the relevant area of the testbed the algorithm provides a high accuracy solution. So, it can be included in the design of the controller as feed-forward contribution to compensate for the effects of the tethers. Based on an implementation in MATLAB and Simulink, this contribution can be run 10000 times in 22:5 s, which means a frequency of 443:29 Hz. The method can represent a valuable help in reducing the gap between simulations and real testing of complex system such as the one represented by a tethered VTVL vehicle.
The model has been developed and can be applied also for other configurations of the tethering. For example, after further maturing of the control system of EAGLE, all tethers on NEST but one could be removed (see https://gnc.dlr.de/s/ yt/eagle-tethered-flight). For this setup, modeling needs to be applied to a single rope only.

Data Availability
The data used to support the findings of this study are included within the supplementary information files.