Inﬂuence of Bending Stiffness on Snap Loads in Marine Cables: A Study Using a High-Order Discontinuous Galerkin Method

: Marine cables are primarily designed to support axial loads. The effect of bending stiffness on the cable response is therefore often neglected in numerical analysis. However, in low-tension applications such as umbilical modelling of ROVs or during slack events, the bending forces may affect the slack regime dynamics of the cable. In this paper, we present the implementation of bending stiffness as a rotation-free, nested local Discontinuous Galerkin (DG) method into an existing Lax–Friedrichs-type solver for cable dynamics based on an hp -adaptive DG method. Numerical veriﬁcation shows exponential convergence of order P and P + 1 for odd and even polynomial orders, respectively. Validation of a swinging cable shows good comparison with experimental data, and the importance of bending stiffness is demonstrated. Snap load events in a deep water tether are compared with ﬁeld-test data. The bending forces affect the low-tension response for shorter lengths of tether (200–500 m), which results in an increasing snap load magnitude for increasing bending stiffness. It is shown that the nested LDG method works well for computing bending effects in marine cables.


Introduction
Cables are used in many marine operations and installations. Examples include lifting operations, mooring lines, Remotely Operated Vehicles (ROV), umbilicals, power cables and towing operations. Snap loads may arise in all of these applications. In lifting operations and in the mooring of ships, the phenomena of snap-back when a cable accidentally breaks during operation is an important safety issue [1]. When mooring lines, power cables or umbilicals become slack (i.e., completely unloaded with zero tension), a snap load of potentially large amplitude is soon to follow when the cable enters back into tension. When such a snap occurs, it is characterised by a sharp rise in axial tension that effectively is a shock wave that propagates in the cable [2]. In the current guidelines, incidents of slack in mooring lines are to be avoided in design calculations [3]. However, in many developing applications such as wave energy converters or floating wind turbines, combinations of device motion, water depth and environmental conditions increase the probability of snap loads occurring [4][5][6][7]. Snap loads have been identified as a potential hazard to the structural integrity of the installation, but the uncertainties in their occurrence, identification and simulation is still high [8][9][10].
Snap loads in marine cables have been extensively studied experimentally [6,[11][12][13][14], analytically [15,16] and numerically [2,[17][18][19][20]. Snap load generation is due to the interaction of the axial tension and the transverse motion of the cable, which is a complex problem with significant nonlinearities [21]. Arguments regarding the influence of bending stiffness on the snap load however vary in the literature depending on the type of cable studied and how it is used or operated.
The most intuitive snap load appears when a mooring cable returns to tension after a period of slack. In this context, cable slack is associated with a loss of axial tension force at any point along the cable. The loss of stiffness in the equation of motion makes the system ill-posed [15] in the absence of bending stiffness, and simulations of a hanging chain were shown to be significantly improved by the addition of a small bending stiffness into a finite difference solver [22]. However, this problem is closely connected with the numerical solution procedure. Many implicit time stepping schemes rely on an estimate of the Jacobian matrix, which becomes singular during free-fall of the cable in the slack regime. Therefore, a small artificial bending stiffness is often recommended to model cable slack in chains, where there is no restoring bending stiffness [23,24]. In explicit simulations, the simulation typically remains stable also in the slack regime, although other restrictions such as the time step size may become prohibitive for very stiff cables [20]. However, with the exception of chains, mooring cables of rope or wire do have a small bending stiffness. The stiffening effect is typically negligible in relation to the working range of axial forces of the mooring in operational conditions, but during slack events, the bending forces may affect the response of the cable. The bending stiffness is also considered important in low-tension applications such as umbilical cables for ROVs, where the cable is designed to operate in slack condition. Cable model developers generally agree that models require bending stiffness to produce accurate and stable simulations of the low-tension cable position [21,25,26]. Although the axial stiffness is significantly higher than the bending stiffness, the forces in low-tension cables are by definition mostly small, and it has been shown by Buckham et al. [27] that transverse and axial forces are of the same magnitude. In the same study, the authors designed experiments specifically for ROV umbilical validation purposes with good results. Two experimental studies hosting both freely hanging cable slack-snap behaviour and bending stiffness effects (i.e., not made on chains) were made by Koh et al. [28,29] on a neoprene rubber cable. In [28], they made a swinging cable experiment, but in contrast to other similar tests (see e.g., [30]), the tension force at the pinned node was carefully recorded.
Numerical models for marine cable dynamics are many and frequently occur in the literature, both with and without bending stiffness included [31]. Examples with bending stiffness accounted for include solvers based on finite differences [2], lumped masses [32], linear finite elements [33,34], curved elements [35], cubic splines (finite elements) [26] and iso-geometric mapping [36,37]. In this paper, the bending forces of the cable are introduced with a new numerical approach: a rotation-free Local Discontinuous Galerkin (LDG) formulation within a high-order Lax-Friedrichs-type discontinuous Galerkin (DG) method.

Scope of the Paper
In our previous work (see Palm et al. [20]), we developed an hp-adaptive DG method for mooring dynamics, where both the mesh size (h) and the expansion order (p) are allowed to dynamically change during the simulation. The method was developed with the aim to capture snap loads with high accuracy. The present contribution aims to extend the capabilities of the solver (named Moody) to include cables with bending stiffness and to begin an investigation into the bending stiffness influence on snap load generation. The paper starts with describing the conservative form of the equations of motion in Section 2, followed by the DG formulation in Section 3. We then present two verification cases in Section 4, followed by a validation against the swinging cable experiment of Koh et al. [28] in Section 5. Section 6 discusses the bending influence on the deep water tether-cage system for ROV operations described by Driscoll et al. [13]. The paper ends with conclusions in Section 7.

Preliminaries
A structure can be approximated as a cable if its cross-section diameter d is significantly smaller than its length, L. The cable position r = [x(s), y(s), z(s)] is at each point expressed in terms of the curvilinear material coordinate s ∈ [0, L]. We denote the inertial (global) system base vectors by x, y and z and, following [34], introduce the Frenet frame coordinatesî, andk as the tangential, normal and binormal directions of the cable, respectively. The Frenet frame is defined from the mean cable line as:î The projection of any vector p onto the cable normal is written p ⊥ = p − ( p ·î)î. p ⊥ is thus in the osculating plane formed by andk.
We will throughout the paper use both standard notations of the time derivative of a quantity:ẋ = dx dt .

Governing Equations
The cable equation of motion is the balance between inertial, internal and external forces on the cable: in which ν = γ 0 v is the cable momentum per meter (mass per meter γ 0 times velocity v), T is the internal tension force and f e represents the external forces. We denote the axial cable strain (elongation) by with elongation factor l = 1 + . The cable tension force can be divided into axial and transversal forces as: where T is the axial force magnitude and T s is the shear force vector. The axial tangent vectorî and the strain can be written in terms of the cable position vector r as: We use ν and q as independent variables in the inertial frame and use the time derivative of Equation (6) to arrive at an expression for˙ q, Equations (4) and (9) can be written in conservative form as: which transforms to: in terms of a single state vector, u T = q ν , a flux vector, F, and a source term, G.

External Forces
The f ex in Equation (4) represents the total body force on the cable segment. The force can be divided into: where f a is the force from the added mass and the Froude-Krylov effect, f b is the buoyancy force, f c represents contact forces and f d stands for the drag force. Added mass f a : In the Morison equation [38], the added mass force on a slender body is assumed to be proportional to the relative acceleration between the body and the flow. Hence, the added mass acting on the cable cross-section is: where A is the cross-section area of the cable, ρ f is the fluid density, a * = a w −˙ v is the relative acceleration of the fluid with respect to the cable, and subscripts ı and ⊥ denote the tangential and perpendicular components of a vector quantity. C mı and C m⊥ represent the added mass coefficients in the tangential and perpendicular direction of the cable, respectively. Please note that there is no dependency of cable strain as we assume a volume-preserving material, in which the cable elongation factor l is cancelled by the same decreasing factor on the cross-section area. Buoyancy, f b : The buoyancy term includes the sum of the cable weight and the buoyant Archimedes force, written as: where g = 9.81m/s 2 is the gravitational constant. For the same reason as in the added mass force, the cable elongation factor does not affect the buoyancy force.
Contact forces, f c : Contact forces refer to the contact between the cable and the ground. This is modelled as a vertical bilinear spring-damper system with dynamic friction in the horizontal direction. The vertical force is computed from the cable penetration depth δ = z g − rẑ and the vertical cable velocity vẑ as: where K g is the bulk modulus of the ground, z g is the vertical coordinate of the ground and ξ is the damping coefficient (ξ = 1 indicates critical damping). The horizontal friction force is proportional to the horizontal velocity of the cable, vx y = v − vẑẑ, up to a threshold speed v µ beyond which the magnitude is constant. The dynamic friction force is computed from: where µ is the fully developed dynamic friction coefficient.
Drag forces, f d : The cable drag is also modelled as in the Morison equation [38]. Based on the relative velocity between the fluid and the cable, v * = v f − v, we compute it as: with d being the nominal diameter of the cable and C dı and C d⊥ being the drag coefficients in the tangential and perpendicular directions. √ l comes from assuming a volume preserving cable material as the cable stretches. The nominal diameter therefore decreases by a factor √ l with increasing , while the cable length factor is by definition l . In combination, the drag force therefore scales with √ l as the cable stretches.

Shear Force Modelling
We model the cable by adapting the local Lagrangian frame formulation described in [2] to the inertial frame of reference. The main assumptions are: (i) that there are no distributed moments acting on the cable; (ii) that the cable cross-section is axisymmetric; and (iii) that the inertial effects of rotating the cross-section can be neglected. The balance of moments equation for a linearly visco-elastic material with bending stiffness EI, bending damping ξ b , and torsional stiffness GI p thus reads: in which we define the curvature κ = ∂î ∂s and its time derivativeκ = ∂ ∂s ˙ q − q˙ l . The bending moment vector is thus acting in the binormal direction. Taking the cross-product from the left on Equation (18) allows us to solve for the shear force as: where M * = M l 2 represents the moment in the stretched domain used for boundary conditions. Finally, we note that the torsional stiffness is the onlyî component in Equation (18), which leads to a simplistic model for the torsional curvature Ω 1 : This is a consequence of all three assumptions (i)-(iii) stated above (see [2] for further information); however, the most important factor is the assumption of a circular cross-section. The torsional moment is constant along the cable, and Ω 1 can be globally computed from the boundary conditions of the cable system at each time step. Torsion has been left outside the modelling in the computational examples of this paper.

Tension-Strain Relations
To complete the set of equations, we need to define the scalar function T( ,˙ ) relating cable strain and its rate of change to the axial tension force. Several material models are implemented in the numerical method (see [39]), but the simulations in the present study are made with bi-linear visco-elastic material properties according to: Please note that the cable does not support compression loads in the simulations of this paper.

Finite Element Method
We use a Discontinuous Galerkin (DG) finite element formulation of Equation (11). The computational domain Ω is partitioned into N elements Ω e ∈ [s e L , s e R ] with size h e . A function x(s, t) is approximated to an arbitrary order P within Ω e as: where φ k (s) is the kth order trial function with expansion coefficientx e k . Legendre polynomials are used as test and trial functions in this study. Defining the inner product operator () Ω e as: the DG formulation expressed in strong form within Ω e reads: The DG method uses a numerical flux (denoted with · in Equation (23)) to express the value of a quantity on an element boundary. In this paper, we use: where n is the outward pointing unit normal, λ is the speed of sound in the cable and F is represented by the Lax-Friedrichs flux of F. Superscripts + and − indicate if values are taken from the interior domain (i.e., from Ω e in this case) or from the neighbouring element of the boundary (i.e., Ω e+1 or Ω e−1 , respectively). See [20] for details. Bending stiffness is introduced via the shear force using a local DG (LDG) [40] approach to compute the derivatives in Equations (19) and (20). We need to define three auxiliary variables to arrive at an LDG formulation for the modal shear forceT s : (i) κ, the cable curvature, (ii)κ, the curvature rate of change, and (iii) τ, the spatial derivative of the moment vector. The modal coefficients,·, of these auxiliary variables are for each element found from: remembering that M * is the stretched domain moment from Equation (20). T s is then recovered from The numerical fluxes are chosen as: where the β ∈ [−0.5, 0.5] parameter governs the amount of flux taken from the left and right side respectively for each equation. In this paper, we used a centred scheme (β = 0) to avoid one-sided bias on the auxiliary variable M * in the results.

Boundary Conditions
The boundary conditions are introduced via numerical flux values at the edges of the finite element domain. The original Lax-Friedrichs formulation without bending stiffness [20] required the numerical fluxes of velocity v or the tension force vector T to be given as the boundary condition on the domain boundaries. The implementation of the nested LDG method for bending stiffness additionally requires that fluxes for the cable tangential vector î or for the total moment vector M * be defined. Finally, if bending damping is included, also the time derivative of the cable direction ˙î is required. Boundary conditions for the auxiliary variablesî or M * associated with the bending stiffness may be specified independently of the conditions set for the cable velocity and tension force at each end point. Variables for which no boundary condition is specified are reactions and are simulated by collecting the flux from the interior domain. Table 1 describe the combinations of boundary conditions required to model some typical end point properties for marine cables. We use the term pinned to describe a point where no moments are transferred. When moments are transferred, we label the condition as a clamped one.
Rigid body connection at point P v

Time Integration
The third-order strong stability preserving explicit Runge-Kutta (RK) method [41] is used for the time evolution of Equation (23). Keeping with the traditions of the RKDG framework, we first let L n = L(u n , t n ) =u n at time index n and allow ∆ t = t n+1 − t n to be the time step size. Then: where u (1) and u (2) are intermediate solution states and u n+1 is used to update the solution at the next time step.

Ring-Shaped Beam
The bending stiffness implementation is verified in this example. We show how a clamped cantilever beam, subjected to a constant edge moment M y = 1 Nm at the free end conforms to a circle. See Table 1 for an explanation of the boundary conditions. The numerical example is taken from Raknes et al. [36], but we extend it here to a static convergence analysis. The beam is initially unstretched with material properties: L = 2π m, EA = 100 N, γ 0 = 10 kg/m and EI = 1 Nm. The moment is smoothly ramped from zero to one over the first 500 s of simulation. Snap shots of the transient beam response at six different times, taken 100 s apart, are shown to the lower left in Figure 1. The radius of the analytical circle is R = 1, centred at [x, z] = [0, −1]. The convergence rates are computed from the numerical error in the L 2 -norm: with r ana as the position of the analytical result.   Figure 1 shows the convergence rates of the simulations. The rates shown in the top right corner table are for each polynomial order P evaluated between the errors at three meshes (h 0 ,h 1 ,h 2 ), where the number of elements were doubled between each mesh (N h 2 = 2N h 1 = 4N h 0 ). The rate is then computed from: where n is the number of degrees of freedom in the simulation. The results in Figure 1 show optimal exponential convergence for P = 1 and for even orders of P with α = P + 1, while sub-optimal (α = P) convergence is achieved for odd orders of P > 1. The simulations were made using centred fluxes for computing the cable curvature and the shear force. Centred fluxes are known to produce sub-optimal convergence rates in LDG methods [42], but are also improving the stability of the simulation. Sensitivity tests with β = −0.5 and β = 0.5 showed an effect only on the absolute error of low resolutions, where β = −0.5 gave the smallest error, as expected due to the suitability for this particular case. The asymptotic converge rates shown in Figure 1 were all unaffected (to three digits) by the β-parameter, except the P = 2 convergence, which was contained in [3.09, 3.13]. From the results of Figure 1, we conclude that the embedded LDG method for bending stiffness is implemented correctly and that it is stable and accurate for very large curvatures and displacements.

Vibrating Cantilever
The nested LDG formulation for bending was also applied to a vibrating cantilever, to serve as a dynamic verification of the bending stiffness implementation. The analytic result of the free vibration of a cantilever Euler-Bernoulli beam is well known. The model equation of a uniform cantilever beam with no rotary inertia reads [43]: where m is the beam mass per m, EI is the bending stiffness, L is the beam length, y is the transverse displacement, s is the axial coordinate and t is time. Spatial derivatives are denoted with , and time derivatives are denoted with an overdot. The boundary conditions (39) and (40) represent a fixed clamped condition at the left side and a free cable end on the right; see Table 1. The initial condition AΦ i is the i th eigenmode (Φ i ) multiplied by a deflection amplitude at the right end (A). The modes of free vibration are written: in which Ω i is the frequency parameter found as the i th solution to the frequency equation: from which the frequency of vibration ω i is derived as: The dynamic solution to (37) is then for each mode of vibration Φ i found as [43]:

A Swinging Cable
The dynamics of a hanging chain or the swinging of a pendulum are frequently used for the verification and validation of new implementations of dynamic cables; see, e.g., [23,36,37,44]. The expected response is dependent on the bending stiffness of the cable. Raknes et al. [36] studied the rigid limit where the deformation of the cable is negligible and the response approaches the theoretical response of a swinging pendulum. On the other extreme, Gobat et al. [23] extensively studied the swinging chain where there is no bending stiffness, although a small artificial value was used to stabilize the finite difference simulations.
Here, we will focus on the experiments of Koh et al. [28] where the swinging motion of a circular rubber cable (25 mm diameter, solid neoprene of ρ c = 1430 kg/m 3 ) with total length L tot c = 2046 mm was investigated. It was suspended with pins a horizontal distance 1805 mm apart at z = 0. The pins were located δ = 12 mm from each cable end, and there was L c = 2.022 m cable between the pins. We will use L c as the cable length in this study. At t = 0, one pin was suddenly removed, and the cable pivoted around the remaining pin location. We apply free cable end conditions on the released pin end and pinned joint conditions at the remaining pin, according to Table 1. Young's modulus of the rubber was estimated to be E = 3.14 MPa, and a pull back test on the suspended cable gave a damping factor of ξ = 0.062 on the standing wave response. The cable properties used in our simulation are shown in Table 2. Table 2. Data used in the baseline simulation of Koh's experiment, computed from ρ = 1430 kg/m 3 and E = 3.14 MPa [28]. The simulated length (L c ) is the length between the two pins (the pivot pin and the release pin), thus neglecting the 12 mm extra cable at each end. Simulations were made with N = 10, P = 4, selected after a mesh convergence analysis. The root mean squared discretisation error in tension between N = 7 and N = 10 elements was less than 1%. The damping factor of ξ * = 0.062 * 2 √ EAγ 0 was in [28] used for a linear damping model proportional to cable velocity, but here, we use it as a material damping in the constitutive model of the rubber material. Consequently, only axial forces are damped in the simulation, and we require an additional parameter for realistic damping of the normal deformation modes. For this, we use a visco-elastic bending moment with damping factor ξ b . The value of ξ b is calibrated via a parameter sweep to provide a reasonably correct level of damping in the system. Figure 3 shows the results of the simulations with different sets of cable properties. We make the following observations:

•
The bending stiffness is important for the cable response. During the first swing in Figure 3a, mostly the end of the cable is affected, but on the return swing in Figure 3b, the EI = 0 cable deviates significantly from the simulations with bending stiffness. In particular, we highlight the "whipping" behaviour of the EI = 0 case, which results in a significant overestimation of the tension force at t = 2 s compared to the experiments; see Figure 3c. The whipping is effectively avoided by the addition of bending stiffness for all configurations studied.

•
Experimental results are reproduced well in the simulations. The peak frequencies match well in all cases studied (see Figure 3d), and the two simulations with bending damping provide a very good approximation to the experimental force reading in Figure 3c. The primary experimental loss factor is probably the frictional losses at the pivot pin, which were not taken into account in the numerical model. This can to some extent account for the overestimation of the primary snap load at t = 2 s.

•
The simulations with bending stiffness require a damping mechanism to reduce the transient bending modes. This is particularly important for the high-frequency response. In Figure 3d, both EI = 0 and ξ b = 0 simulations show significant forces around the return period of longitudinal waves in the cable T r , which are not present in the experimental readings. Please note that the high-frequency oscillation of, e.g., ξ b = 0 in Figure 3c is not numerical noise. It is the long-lived elastic response of the bending modes due to insufficient damping properties. The results are much improved by the addition of ξ b , and less sensitive to its value. Cable experiments in air are known to be sensitive to the damping properties of the setup and behave very differently from experiments in water [45]. We therefore expect the importance of parameter ξ b to be small to negligible in marine application of cables.

Snap Loading in a Deep Water ROV System
We consider also the case of a deep water operation with the Remotely Operated Platform for Oceanographic Science (ROPOS) device, which was studied in a series of publications by Driscoll et al. [13,16,19]. The ROPOS system consists of a heavy cage suspended from the operational ship with a tether. The cage is in turn connected to the ROPOS with a 30 m long slack umbilical cable. As the umbilical from the cage and ROPOS was slack at all times, said cable and the dynamics of the ROPOS can be neglected in the analysis of the tether-cage system response [19]. In [13], snap load in the tether was observed in several tests. A one-dimensional (vertical motion only) lumped mass model was developed and calibrated to analyse the system numerically in [19], showing very good agreement given the uncertainties of a field-test environment and the scale of the operation (operational depth L = 1730 m). In this section, we first make a new calibration to achieve similar snap loads for one of the documented events in [19], after which we will analyse the influence of bending stiffness at different operational conditions and depths.

Calibration and Validation
The cage is a 2.1 by 3.4 by 4.1 m steel structure, with the tether attachment at the top. Previous models in [16,19] studied the vertical response of the cage, leaving surge and pitch out of the modelling. Hence, the centre of gravity is not specified, and we assume it to be at the geometric centre of the cage. The cage was modelled as a submerged cylinder with the same lid area as the cage base. The numerical values and notation of the parameters are presented in Table 3. Please note that as the cage was open to the surrounding water, the added mass is higher than common values for solid bodies. The calibration procedure used by Driscoll et al. was well described in [19], where both the empirical/manufacturer's data and their calibrated values were presented. The calibration was made for operational depth L = 1730 m. We used the calibrated values for most parameters; however, in our model, we got a better match with the experiments when the empirical (higher) value of 9200 kg for the cage added mass was used. The calibrated lumped mass model in [19] did not produce snap loads as expected when the recorded ship motion was used. Therefore, the amplitude was increased by 10%, and the slack-snap condition was reproduced with satisfactory results. In our case, the heave motion of the fairlead was digitized from Fig. 5a in [19], pass band filtered in [0 − 0.8] Hz, and assumed to be at the still water level. We then needed an amplitude increase of 20% to get the appropriate extreme event behaviour. The final calibration was made with the vertical drag coefficient, which was changed from the recommended value of 2.3 in [19] to 2.7. The accuracy of the Morison drag approximation on the cage was also debated in [19] where the addition of a wake model substantially improved the simulations compared with the experimental data.
The results of the calibration are shown in Figure 4. Figure 4b shows the important difference between the tension at the cage (bottom line node) and at the ship (at the top node). The tether self-weight in water is 44 kN (at L=1730 m), which explains the offset between the cage and ship tension. Consequently, the tether is only slack at the bottom section, and otherwise in tension by its own mass. The snap load originates from the cage and takes L/c=0.445 s to propagate the tether length, where c = 3870 m/s is the speed of sound at which snap loads propagate in the cable [2,19]. This explains the time lag between the tensions at the cage and the ship. The snap load is gradually damped by the internal material damping factor ξ and by tangential drag forces, which become increasingly important for longer cables. The peak loads of Figure 4 are more long-lived in the numerical results than in the experiments, which suggests that higher damping factors could have been used. However, we consider the results to be acceptable given the uncertainties related to comparing numerical simulations to field-tests.  The simulations were made with N = 20 elements of order P = 4, selected based on a mesh sensitivity analysis presented in Appendix A. The mesh dependence was tested for the longest tether, as this constitutes the coarsest case. For shorter lengths of line, the resolution per meter of tether increases.

Bending Influence
The numerical and analytical models of Driscoll et al. [16,19] focused on the one-dimensional system including heave of the cage and the longitudinal dynamics of the tether. We now extend the investigation to include the transverse motion of the tether and cage system. Sixteen combinations of cable length (operational depth) and bending stiffness were studied, combining four lengths and four values of EI according to: where EI * = 2.62 kNm 2 is the bending stiffness of the field-tested tether computed from the cable area (0.25πd 2 ) and the axial stiffness EA = 45.5 MN [19] under the assumption of a solid circular cross-section. The fairlead motion is a circular orbit of radius a = 2 m at frequencies f = [0.10, 0.20, 0.25] Hz. The amplitude was chosen from the experimentally reported maximum of 4 m peak-to-peak, and the frequencies were chosen to encompass the typical wave band of the field-tests [13]. Figure 5 shows the tension at the fairlead for different combinations of line length, bending stiffness and excitation frequency. We make the following observations: • We confirm the predictions of the analytical model in [16] that slack-snap occurs at f = 0.2 and f = 0.25 Hz with this amplitude, while the f = 0.1 Hz case remains in tension.

•
There is no noticeable difference in system response due to bending stiffness for the largest operational depths L 975 and L 1730 . For L 500 and L 200 , a small but increasing difference can be seen in the snap load generation.

•
The small effect of bending stiffness can be observed in Figure 5d,f. The lower the bending stiffness value becomes, the earlier the generated snap load appears at the ship after each period of slack. The difference between the simulations increases with decreasing cable length. A more detailed analysis of the slack events is needed to explain the differences in snap load generation between the simulations of different bending stiffness. Figure 6 shows the tension force and vertical acceleration at the cage together with the position and tension along the whole cable around a snap instant of the L 200 simulation at f = 0.2 Hz. The α 0 case becomes taut slightly earlier than the stiffer cables (see Figure 6a, where the tension at the cage and at the ship is shown). This is consistent with the behaviour of the vertical cage acceleration in Figure 6b. In Figure 6b, we also note that the acceleration of the fairlead is positive and increasing during the first snap instant after the slack. The relative velocity between cage and fairlead is therefore increasing. The snap load magnitude is strongly affected by the relative velocity at the snap instant [12], which explains why the α 0 case has a snap load of smaller amplitude than the stiffer cables. The smaller snap load of α 0 in Figure 6a is almost fully dissipated at the peak of the load cycle, in contrast to the cables with bending stiffness where a noticeable snap load still propagates the cable. In Figure 6c, the tension distribution along the tether is shown for three times: t 0 (black), right before the first snap takes place; t 1 (red),short after the α 0 snap load has occurred; and t 2 (blue), when all simulations have become taut at the cage connection. The differences in snap load size and phase are evident from the results. Finally, the earlier tension response is explained by how the bending stiffness influences the cable position near the cage. As is seen in Figure 6d, the differences in position between the different α values consequently matches the difference in snap load appearance in Figure 6c. In summation, the inclusion of bending stiffness in the analysis does not change the overall tether response significantly for these load scenarios. However, it affects the slack regime dynamics of the tether to the point of changing the resulting snap load event. The point in time when the snap occurs changes with the exact position of the slack tether, which in turn affects the resulting snap magnitude. Our results show that the slack regime forces near the cage constitute the key parameter that governs the snap load. Driscoll et al. [13] described how the lowest section of the tether needed to be replaced after some time in operation. It had experienced bending fatigue, and the connection to the cage was lost. A bending stiffener at the connection to the cage was later introduced, which alleviated the issue. This failure is another indication that the end section of the tether is important. The effect of a transverse current (with only vertical fairlead displacement) was investigated numerically by Zhu et al. [46], where they found that a small current increases the severity of the slack-snap behaviour. However, as the current velocity was further increased, the fluid drag force began to dominate the slack tether behaviour so that slack periods became less frequent and the snap load size decreased. All results thus point to snap load generation depending on the intricate force balance between tensile stress, bending effects and fluid forces in the low-tension and slack region of the tether. We leave it for future work to more elaborately study the multitude of potential configurations and load scenarios that may or may not affect the low tension dynamics of cables, and consequently govern the ensuing snap load.

Conclusions
This paper presented a new implementation of the equations for cable dynamics. The Discontinuous Galerkin (DG) formulation of [20] was extended by the addition of a nested LDG method to include bending stiffness in the modelling. Verification studies on the method showed that the implementation works and is stable for both small and large deformation in dynamic and transient cases. The deformation of a beam into a ring showed exponential convergence of order P for even orders of the expansion base and P + 1 for odd orders. Comparison with a swinging rubber cable showed a good match between experimental and numerical results. By comparing results with and without bending stiffness included, significant changes in cable response were noted. The whipping of the flexible cable is reduced by the addition of bending stiffness, resulting in significant improvement of the match to the experimental data, as well as a reduction of the peak snap load. Because the experiments were made in air, there is little drag damping, and the results become more sensitive to the material damping properties of the simulation [45].
We also validated the code against field-test data of a vertical tether with a heavy cage in deep water. The dynamic behaviour of the tether tension during a snap load event in the experiments was well captured by the model; however, there are uncertainties in the amplitudes of motion used in the calibrated model. When forced to a circular motion at the fairlead, the tether behaviour with a varying amount of bending stiffness was investigated numerically. We found that the overall response is essentially unchanged over the range of depths tested L ∈ [200, 1730] m; however, the snap load generation was significantly affected by the difference in bending stiffness. The snap load of the perfectly flexible tether occurred earlier and with a smaller amplitude than in the cases with bending stiffness included. This is explained by how the bending stiffness changes the position of the tether section closest to the cage where the cable is in slack or low tension during a period of time. The bending forces and tensile forces there are of similar magnitude [27]. Our results indicate that the exact configuration of the tether prior to the snap affects the snap load generation process. Further work is needed to investigate in more detail during which conditions this difference is of significance to the fidelity of marine cable design calculations. Therefore, we strongly recommend that the sensitivity to bending stiffness influence be accounted for in design calculations where part of a tether or cable is in the slack or low-tension regime and snap loads may occur.
Author Contributions: J.P. and C.E. are the sole authors of the study. The following division of tasks was used in the paper: conceptualization, J.P. and C.E.; methodology, J.P. and C.E.; software, J.P.; validation, J.P. and C.E.; formal analysis, J.P.; investigation, J.P.; resources, J.P.; data curation, J.P.; writing-original draft preparation, J.P.; writing-review and editing, J.P and C.E.; visualization, J.P.; funding acquisition, J.P. and C.E. All authors read and agreed to the published version of the manuscript.
Funding: This research was funded by the Swedish Energy Agency under Grant Numbers P42246-1 and P47264-1.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Deep Water ROV Simulation: Mesh Resolution
This section describes the mesh resolution tests made for the ROPOS system tests described in Section 6. We selected an expansion base of order P = 4 for the mesh sensitivity tests and used a forced circular motion of a = 2 m, f = 0.2 Hz on the fairlead to test the system sensitivity to mesh resolution. The mesh dependence tests were done on the longest tether (L = 1730 m). Four resolutions were tested, with the number of elements N = [5,10,20,40]. The resulting tensions at the ship (surface) end point of the cable are shown in Figure A1. From the results, we concluded that the N = 20 mesh constitutes a high-resolution simulation. Little new information was gathered from the N = 40 mesh, and we therefore chose to stay with N = 20 to save computational time.