Ocean Engng, Vol. 10, No 6, pp.443--457,1983. Printed in Great Britain.
0029-8018/83 $3.00 + .00 Pergamon Press Ltd.
NUMERICAL SIMULATION OF UNDERSEA CABLE DYNAMICS C. M. ABLOWand S. SCVI~.CnTER SRI International, Menlo Park, CA 94025, U.S.A. Ah~rm~t~A fully three-dimensional code has been written to compute the motion of a towed cable. The code is based on a robust and stable finite difference approximation to the differential equations derived from basic dynamics. A 3500-ft (1.07 kin) cable pulled at 18.5 knots (34.3 km hr -j) through a circular turn of 700 yd (0.64 km) radius has been computed in about half of the real time of the maneuver. The computed displacements are close to the measured ones; the changes in depth are within 2%.
INSTRUMENTED cables are towed behind ships or submarines for sound detection in geophysical exploration, antisubmarine warfare, and similar pursuits. Accurate location oi the source of a sound requires accurate knowledge of the location of the receiver. When the location cannot be observed, it can be estimated from a model that computes cable motion under known ship maneuvers. A review by Casarella and Parsons (1970) pointed out a decade ago that predictions of even two-dimensional motion were unsatisfactory as to mathematical method and hydrodynamic loading function. A second survey by Choo and Casarella (1973) concluded that "there is a need for an easy, simple methoO that can solve any unsteady-state problem to a good accuracy, and yet require only a small amount of computation time." This report describes such a method. In this numerical simulation, the cable is treated as a long, thin, flexible cylinder in arbitrary motion under the action of gravity, a maneuvering towship and hydrodynamic loading. The dynamic m o d e l is a finite difference a p p r o x i m a t i o n to the three-dimensional differential equations for conservation of momentum. The loading function is taken to be the sum of independently operating normal and tangential drags, each given by a single coefficient. The numerical solution method is a completely stable, implicit, finite difference scheme. The present model is similar to models published by Paidoussis (1966), Choo and Casarella (1972) and Russell and Anderson (1977) each of which covers only a restricted class of motions. The model by Sanders (1982) applies to arbitrary maneuvers but differs from the present one by using a finite element approximation to the dynamics and a conditionally stable numerical method. Computation of a circular turn took approximately real time with the Sanders model but only half as long with the present one.
VECTOR DYNAMIC EQUATIONS The balance of forces at a point of the cable may be written 0 0--~T+W+F+B=0 443
444
C.M. ABLOWand S. SCHECHTER
where S is distance along the stretched cable, T is the tension, W the weight less buoyancy per unit length, F the force exerted on the cable by the fluid per unit length, and B the d'Alembert force of the cable motion per unit length. The bold face indicates a vector quantity. Let (i, j, k) be a fixed inertial Cartesian frame of orthonormal vectors with k pointing vertically downwards. Then W=
wk/S'
,
w = (m -
p A )g
,
where m is the mass per unit length, A is the cross-sectional area of the unstretched cable, p is the fluid density, and g is the acceleration of gravity. The prime denotes the partial derivative with respect to the unstretched cable length coordinate s, holding t fixed. Partial derivatives with respect to t holding s fixed will be denoted by a dot. Let r be the position vector for points on the cable r=xi+yj
+zk
The d'Alembert force of cable inertia is
The effect of the surrounding fluid is to add the mass ( p A / S ' ) per unit length to the cable inertia for sideways motions. Thus B =
[J + ( u . t ) t ] F s ' ] ' ,
= m + pA,
where t is a unit vector tangent to the cable, t
=
r'/S'
,
S' = (x '2 + y , 2 + z,2)½
,
J is the current, the velocity of the sea with respect to the fixed coordinate frame, and U =
i--j. Tension T has the direction of t and has amplitude T, a function of the strain ¢: T=Tt
,
~=S'-I.
A linearly elastic cable satisfies =
eT
, e
=
1/EA,
where E is Young's modulus. An inextensible cable satisfies the same relation with e = 0. For a cable with circular cross section, the drag force F per unit length can be written F = -(1/2)pd{'=C,(U-t)lu'tlt + C, IU - (U't)tl [U - (U.t)t]}
Motion of a towedcable
445
where d is the cable diameter, d = (4A/'rrS') ~
,
and Ct and C. are the tangential and normal drag coefficients. INTRINSIC BASE VECTORS Introduce unit vectors orthonormal frame with perpendicular to t will be The two frames can be
n and h at each point of the cable so that (t, n, b) is an the same orientation as (i, j, k). The direction of n or b specified later. related by rotations through the Euler angles (0, ~, 6).
(t, n, b) -- (i, j, k)
A =
B
C
=
=
cos0
sin0
0
-sin0
cos0
0
0
0
1
1
0
0
0
cos~
sinO
0
-sin~
cosO
cos6
-sin6
0
sin6
cos6
0
0
ABC
=
ABC
0
1
cos0cos6 + sin0cosqJsin6
cosesin6 + sinOcos~cos6
- sin0cos6 + cosOcos0sin6
- sinesin6 + cos0cosOcos6
- s i n ~ sin6
-sinO cos6 sinesin~
1
cosOsin~ /
cos~ ]. The formulas represent a rotation through angle ( - 0 ) about the k axis to bring the i axis into the plane of t End n, rotation about the new i axis through ( - ~ ) to bring k into coincidence with b, and rotation about b through 6 to bring i and j into coincidence with t and n. The signs have been chosen to permit agreement with the report by Rispin (1980a).
446
C.M. ABLOWand S. SCHECHTER
Differentiation with respect to s gives (t', n', b') = (t, n, b) (ABC) -1 ( A ' B C + A B ' C + A B C ' ) For these matrices the inverse equals the transpose so that the multiplications can be readily carried through to give
(t', n', b') = (t, n, b) (e'M0 + ~o'm, + ~b'm,),
Mo
=
0
cos~
sin~ cos~b
-cosd)
0
-sin@ sin~b
-sin,
M,
M,
cos,b
0
0
sin~
0
0
coscb
-sin~b
-cos~
0
--
--
0
sind~ sin~b
0
-1
0
1
0
0
0
0
0
The same relations are valid for derivatives with respect to t. The inertial coordinates (x, y, z) of a point are related to the Euler angles (e, ~0, ~b) and the strain t through r' = x ' i + y ' j
+z'k=(l+~)t
so that x' = (1 + ~) (cose cosd~ + sine cos~ sin~b) , y' = (1 + e) ( - s i n e cos~ + cose cos~ sin~b) , z' = (1 + ~) ( - s i n e sind~) Integrations in s determine (x, y, z) when ~ and the angles are known. If one puts i- = v , t + V,,n + V~,h
,
then
(V,, V,,, Vb) = (x, y¢, ~) A B C
Motion of a towed cable
447
or
(~, 5',~) = (v,, v,,, vb) (ABC) -I so that an integration in t recovers x, y, or z from known functions Vt, V,, and Vb. Let components of the relative velocity U in the two systems be defined by U -- Uli + U2j + U3k -- U# + Unlit + Ubb
Then
(u,, on, ub) = (u~, 02, 0 3 ) a ~ c and
(u1, u2, u3) = (v. v., ub) (aBC) -1 The same relations hold for the components of the current J.
SCALAR DYNAMIC EQUATIONS The vector dynamic equation can be written as a system of three scalar equations by taking components in three independent directions. Components in the directions (t, n, b) give T' - w sin0 sin, - (1/2)pd(1 + ¢)t'trCtUtl Url
= me', + (mlV~ - o A J . ) ( c o s * - mVM(I+¢) T(-cos00' + , ' ) - w
6-$)+(mlVb-pAJb)sin* cost 6 + sin,,, ~)
,
sin0 c o s , -(1/2)od(l+e)tC,,U,,(U~
+ U2) t
= mVt(-cosO 6 + 4 ) + m l f / , , - P A ) , , + ( m t V b + p A J b ) ( - s i n t ~ - (mlV,,-pAJ,,)id(l+¢)
sin,0+coM~ t~)
,
- T ( s i n ¢ cos,0' + s i n , 0 ' ) + w cos0 - (1/2)pd(1 + e)tC,,Ub(U2., + U~) t = mVt
( - s i n , c o s , {J - s i n , t~) +(mlV~-pAJ.)(sinO sin, 0 - cos¢ ~) + m l V b -- pA)b -(mlVb-pAJb)#./(1
+ e)
where U, = V, - J,, and similarly for the other components of U. From
(r')" -- (~)'
448
C.M. ABLOW and S. SCHECHTER
or
[(1 + e)t]" = (V,t + V,n + Vbb)' one obtains the compatibility relations = Vt' + Vn(cosd/O' - ~b') + Vb(cos6 sinqJ O' + sin6 to') , (1 + e ) ( - c o s 0 O + ~b)-- V~' + V/(-cos0 0' + ~b') + Vb(-sinto sin6 0' + cos d~ to') , (1 + ~ ) ( - c o s 6 sinto 6 - sinqb d/) = Vb' - V/(cos6 sinqJ 0' + sin6 tO') + V,(sinqJ sin6 0' - cos6 to') . BOUNDARY AND INITIAL CONDITIONS The previous section presents six scalar dynamic differential equations of the first order in space variable s and time variable t. Initial conditions for the system require specification of all dependent variables as functions of s. If the functions satisfy the equations, they represent an initial steady streaming. If not, the equations require some time derivatives to be different from zero initially so that the solution begins with an unsteady motion. Six boundary conditions are needed for all time. The time-varying position of the towpoint and therefore the three components of the velocity vector at that point will be known. Another condition is that the tension T is zero at the free end. Two remaining conditions can be obtained from the second and third of the dynamic equations under the assumption that the cable is smooth near its free end. However, simpler approximate end conditions described below were used in the calculations. INTRINSIC COORDINATE SPECIFICATION The cable dynamic equations are six equations in seven unknowns. In the derivation, the orientation of the local coordinate system was incompletely specified so that an additional condition is needed. Let that condition be f(y,s,t) = 0 where f is a scalar function of the time t, the distance s along the cable, and the 7-dimensional vector y, the components of which are the unknown variables. The cable dynamic equations are written My' = Ny + q where y = (T, Vt, Vn, Vb, 0,6,to) r An additional time-independent condition fits this format if the row vector of the partial derivatives of f with respect to the unknowns is added to M, a row of zeroes is added to N, and the one element row added to q is the partial derivative of ( - j 0 with respect to s. The given condition is recovered from this added row if proper conditions are given at s --~0.
Motion of a towed cable
449
m
The solution method involves inversion of M, the matrix M as augmented. For this, the determinant of M should be different from zero. According to the equations in the preceding section, 1 0 0 0 0 0 0
0
0
0
0
- TCO
T
0
0
0
0
0
- TSOC6
0
- TScb
0
1
o
o
v . c o + VbSOC*
- V . VbS6
0
0
1
0
- V t C O - VbSOS6
V,
VoC6
0
0
0
1
-SO(VtC~-
0
-V36-V.C,
fl
.t:2
f3
f4
f5
.[:6
.1:7
m
M
=
V.S~)
where S0 stands for sin(0), C0 for cos(0)~ and)~ is the partial derivative of f with respect to the i-th component of vector y. The determinant of M may be computed to be det M = T2 [SO(-f3Vb + f4V. - fs + f7C6) - Sd~COf6] m
One sees that det M vanishes where tension T is zero. This is probably an effect of the physics of the cable treated as a completely flexible string. Det M also vanishes when the bracketed factor is zero. No choice of f has been found that prevents this factor from vanishing. The simple choice f = 0 - ,d2, so that 0 = 'rr/2, reduces the factor to C6. This quantity vanishes only when t is parallel to k, i.e. when the cable is vertical. Solution for y' can be made if the cable is under positive tension and out of the vertical at all points. EQUATIONS IN MATRIX FORM The choice of 0 = ~r/2 .leaves 6 unknowns, the components of the new vector y:
y = (T, Vt, V.,Vb, O,cb)r The dynamic equations may be written y'--M -1(NiT+q) where the equations are reordered in a more convenient form that makes the new M upper triangular. Then
M
=
1
0
0
0
0
0
0
1
0
0
Vbco~
-V.
0
0
1
0
- Vbsin~
V,
0
0
0
1
V . s i n 6 - V,cos~
0
0
0
0
0
- Tcos6
0
0
0
0
0
0
T
450
C.M. ABLOWand S. SCHECHTER 0
0
0
0
0
0
1
0
0
Vb/T
V,,/T
0
0
1
0
-(Vb/T)tand~
-WT
0
0
0
1
(V~tand~-V,)/T
0
0
0
0
0
- l/Tcosd~
0
0
0
0
0
0
1/T
1
M-t
N =
Vb--pAJb)COStb
-meV/(l+eT)
m
0
0
(ml
e
0
0
0
0
0
0
0
0
0
0
l+eT
0
0
0
0
-(1 +eT)cos6
0
-e(mlVb-pAJb)l(l+eT) 0
0
m~ (ml V . - p A J . ) s i n 6 - m V t c o s 6 0
-e(mlV.-pAJ~)l(l+eT) 0
m1 0
- ( m t Vb--pAJb)sinrb
- ( m I Vn - [ a J n )
mV,
and wsind~ + (1/2 )pd( l + eT)iwCtUtt Utl 0 q
=
0 0 (ll2)f~l(l+eT)*C.Vb(U2~ + U2)i - pA)b
wcos~ + (1/2)pd(l+eT) ~ C.U.( U2 + U2)~ - pA).
Initial conditions in the form of given functions of s for each of the variables may be prescribed arbitrarily. Boundary conditions at the towpoint are obtained from its given velocity as a function of the time. The other boundary for the calculations has been taken at a point P, a short distance from the free end, to avoid the singularity in M at T = 0. Between P and that end, the cable is assumed to be straight. The derivatives with respect to s of the two angular variables are therefore zero at P. The tension T at P is taken to be what the first scalar dynamic equation gives for a straight rod by integration from T = 0 at the free end. THE NUMERICAL MODEL A mathematical model of towed cable dynamics provides a way of determining the location and tension of any material point of the cable as a function of the time. Points of the cable are designated by their distance s from the towpoint in the unstrained cable so that the equations of the model are partial differential equations with s and the time t as independent variables. There are three second-order, coupled, hyperbolic equations for the components (x, y, z) of the position vector in a fixed cartesian coordinate frame together with an equation relating the tension T to the strain. Because the force components normal and tangential to a section of the cable arise from different physical processes, use of local coordinate frames with an axis tangential to the cable tends to decouple the equations and facilitate their numerical solution. The
Motion of a towed cable
451
equations of this kind of model are conveniently written as a system of first-order partial differential equations for the Euler angles (0,~,~) describing the local orientation of the cable, the velocity components (Vt, V,,, Vb) in the local coordinate frame, and the tension T. One of these dependent variables is eliminated when the direction of one of the axes of the frame normal to the cable is fixed. The vector y of the dependent variables has therefore 6 components. The tint-order system to be solved for y has been written above in the form of a matrix differential equation: y' = M - I ( N y + q)
To avoid restrictions on the accuracy or applicability of the computation, we use a finite difference method in which the governing differential equations are replaced by algebraic difference approximations. An implicit method is used in the time. The calculation will then be stable for time steps much larger than those permitted by the Courant-Fredrichs--Levy wave speed condition (Richtmyer and Morton, 1967). At each time the equations present a boundary value problem in the spatial coordinate. Shooting methods for solving this problem are not likely to be successful because of the singularity at the free end of the cable. An appropriate computational scheme has therefore been set up to solve simultaneously for the whole set of values of the variables at the mesh modes along the cable for a given time. FINITE DIFFERENCE APPROXIMATION Assume the dependent variable to have known values Y0 at time to. Let At be the time step to later unknown values y, which are written without subscript, for brevity. A finite difference approximation to the matrix equation of second-order accuracy, centered in the time interval, is y'+y~-f=0 where f = (M-1N+M~lNo)(y-yo)/At
+ M-tq+M~'tqo The similarly centered discretization over the spatial step As from y~ to Yi+l gives (Yi+I - Yi + Y0,i+I - yo~)/As - (fi+l + [i) = 0 i = 1,2 ..... n
,
,
where n is the number of spatial intervals from one end of the cable to the other. There are 6n scalar equations for the six unknown components of y at the ( n + l ) nodes of the difference mesh. Six additional equations are provided by the boundary conditions. Let g be the vector of functions set to zero in the difference equations plus boundary conditions. Vector g has 6(n+l) components, each component being a function of the
452
C.M. ABLOWand S. SCHECHTER
values of the dependent variables at the mesh points at time level t~ + to + At. The equations to be solved then read g=0 where the components of g are functions of the Yij (i= 1,2 ..... n+ 1; j = 1,2 ..... 6). These are nonlinear equations to be solved by iteration. Denote a current approximate value for the solution by Yij and an improved value by yij + 8~/. The 8~j are found according to Newton's method by solving the linear equations JS+g=0 where J is the Jacobian matrix of the derivatives of g with respect to the yq. Analytic expressions were derived for the entries in g and J used in the computation. Let the components of g be arranged with those presenting the boundary conditions at one end of the cable tint, those for the first interval next, and so on to those for the last interval and then those for the boundary conditions at the other end of the cable. Matrix J is then of block tri-diagonal form because the tirst six scalar equations involve only values at points 1 and 2: the next six at points 1,2,3; the next six at 2,3,4; and so on with the laot-Six equations involving only values at points n and n + 1. The equation can be solved with inversions of matrices of only sixth order, the order of the submatrices of J, as shown by Schechter (1960). RESULTS The code TCABLE was written according to the above specifications. At each time step, the code produces the three-dimensional shape of the cable, the tensions along the cable, and the velocity components at each node, as well as the angles made by the cable with respect to coordinate directions. The code is self-starting from a straight-line run, but allows for a restart from a previously presented initial shape. The direction of the towing path and the towpoint speed may be changed at any time. Some automatic checking that the results fall within physically realistic bounds is included. In the variety of tests given the code, it has exhibited unusual robustness. A reasonable rate of convergence was maintained with large real time steps of as much as 30 sec and a sharp drop in towpoint speed from 18.5 to 10 knots (35 to 18 km hr-1). The number of the Newton iterations needed for convergence at each time step varied from 2 to 12, depending on the maneuver and the type of change in the maneuver. The iterations increased with the speed of the towpoint and the length of time step as well as with the change in curvature of the towing path. The specific tests performed were for realistic full-scale data on straight-line runs and turning maneuvers as supplied by Rispin (1980b) of the David W. Taylor Naval Ship Research and Development Center. The cable system used in these tests is illustrated in Fig. 1. Lengths of the towcable; array, drogue and positions of special points in the array section are given in Fig. 2. Other physical properties are given in Table 1. The towcable section in most of the test runs was subdivided into ten finite-difference subintervals. Some runs were made with 20 subintervals to check the accuracy and stability of the difference scheme, but the smaller number of intervals was found to be sufficient.
Motion of a towed cable
I DrogueI Array I
453
Towcable JA-Ia64-1
FIG. 1. Towed array system.
Towpoint dll
.._
E
8
-~30.i[~.7 T 156 T 71 ~ 8 . 2 ~ / /
F1G. 2. System subdivision lengths (m).
TABLE 1. CABLESYSTEMPHYSICALPROPERTIES Segment
Diameter (in.)
Length (ft)
Wt in H20 (lb)
C.
Ct
Rope drogue Array
1.0 3.125
100 900
0.039 0.0
1.8 0.02168 1.8 0.00898
Towcable
1.6
2372
0.16
2.0 0.0150
To avoid the material discontinuities at the junctures of the array with the towcable and drogue, dummy mesh points were inserted. For a straight run at 10 knots (18.5 km hr-1), the angle of the towcable with the horizontal at the towpoint was computed to be 2.6 ° as compared with 2.7 ° in data reconstructed by Rispin (1980b) from another cable program. Our primary test of the efficacy of the code was to apply it to a loop turn maneuver, since experimental data are available for a similar case. The code was able to model this motion in reasonable agreement with experiment. We believe this is the first dynamic computer program to do so in less than real time. A run was made with towpoint speed at 18.5 knots (34.3 km hr -1) with the towpoint running straight for one second, then into a circle of radius 700 yd (0.64 kin) for 440 sec,
454
C.M. ABLOW and S. SCHECHTER
and finally, coming out of the complete 375 ° turn, to run straight along a tangent line of 15° for 300 sec. A plot describing the x, y coordinate positions of points A and C is given in Fig. 3 and a plot of depth vs time for these points is given in Fig. 4. The total time of the maneuver was 741 see, whereas the VAX-11/780 computer time was 420 sec. Although the experimental data given by Rispin (1980c) is for a loop turn that does not contain a fixed radius of curvature, it is nevertheless almost circular with radius close to 700 yd (0.64 km). Our data compare favorably with his experimental data. In particular, the changes in depth of point A of the array as given by experiment and code differ by less than 2% as shown in Table 2. This run was made with a total of 34 time steps, with a time step of length 20 sec in the circular phase and a time step of 30 sec in the final straight phase. Figure 5 compares measured and computed locations of point C early in the turn, plotted with a large scale. The discrepancies appear to be within the bounds of precision of drag coefficients, experimental measurements and loop turn asymmetry. The code was run both for an inelastic cable and for an elastic steel cable. The results were substantially the same. The code was also run with an additional straight run added to the loop but with an abrupt change in speed down to 10 knots (18.5 km hr -1) from 18.5 knots (34.3 km hr-~). No instability occurred, although the number of Newton iterations needed for convergence rose at first from 3 to 12. Subsequent iterations converged in three steps. As expected, the depth of cable went down to its steady state level in about 300 sec.
311170
I
I
l
I
I
I
I .....
I
I
1
2500
2000
200
1500 u,J
Q
800
O 30OO
ol
500 -
800 800 650 A •
2OO
i~4~0~,.,
~O100
==n
0
0
O --
0
o100
-500 --
-1~
I -1500
-1000
I
I
-500
0
1
I 1000
I 1~
I 20OO
I
t
I
2~
3~
3~
DISTANCE(m) JA-lS54-3B
Fro. 3.
Loop turn, plan view, with times (sec) and location of tow~int, O, and of computed points A, 0 , and C,A.
Motion of a towedcable I
4
-
2
-
o
I
I
~
0
I
0
Straight -2
~-
I
Circular
I
0
455 I
I
1
0
0000
Turn
Straight Run
Run
"r Q. UJ
-4
r~
% A
% "8
A
%
--
~a
eA -10 -
A
eA
-12 -200
¢,
I
0
- 100
100
I 200
I 300 TIME
I 400
I 500
I 600
I 700
800
(sec) JA-1654-4A
Flo. 4. Depthsin loop turn of towpoint, O, and computedpoints A, 0, and C, A.
TABLE2.
COMPARISON OF MEASURED AND COMPUTED DEPTHS
Computed depths
Rispin (1980c)
Code TCABLE
Initial depth (m) Minimum depth Difference
10.04 2.51 7.52
10.95 3.54 7.41
Final depth
10.16
10.82
CONCLUSIONS Since the system of nonlinear hyperbolic equations governing the dynamics of an undersea cable was not solved analytically, we relied on experimental data and previous numerical solutions to evaluate the code. In view of the close agreement of the program output with the given experimental data for a loop turn maneuver, we feel confident in the efficiency and versatility of TCABLE. Moreover, the data are produced in less than real time, so we may assume that the method is economical and practical.
456
C.M. AaLow and S. SCHECHTER I
t
I
1
I
I
220Z~k q
R
"O
g 5
200 Z ~
al
I C3
-
18o ~
k-
Z IJJ
t&l
UJ t~ ¢O
rv-
160~ _
1%=
IJJ MJ
--
120
•
4~
•
100 --
a0
•
az~
A40
[
[
1
L
l
t
100 METERS BETWEEN TICKS m azimuth 0 degrees JA-364525-1B
FIG. 5. Loop turn, plan view of locations of point C as computed, L~, and observed, &, by Rispin (1980c), at times (sec) shown.
The method of solution llas proven to be quite robust in that it responded without fail to changes in input. This is exemplified by aorupt changes in towspeed from 18.5 to 10 knots (34.3-18.5 km hr-t), changes in time step from 1 to 20 sec, and changes from a straight run to a circular one. The robustness is also exemplified in the large variety of spacial steps along the cable, which vary in length from 513 to 27 ft (156-8 m). Any changes in cable diameter and material properties are likewise taken in stride. The code could therefore be used aboard the tow ship to predict the behavior of the cable system and so influence operational decisions. Acgnowledgemcnt--Support of this work by the David W. Taylor Naval Ship Research and Development Center under project number N00140-82-C-CE61, P. Rispin, technical monitor, is gratefully acknowledged.
REFERENCES CASAI~ELLA,M.J. and PARSONS,M. 1970. Cable systems under hydrodynamic loading, Mar. Technol. Soc. J. 4, 27-44.
Motion of a towed cable
457
Csoo, Y. and CASARELLA,M.J. 1972. Configuration of a towline attached to a vehicle moving in a circular path. J. Hydron. 6, 51-57. CHOO, Y. and CASSJ,ELLX, M.J. 1973. A survey of analytical methods for dynamic simulation of cable-body systems. J. Hydron. 7, 137-144. PAIDOUSSlS, M.P. 1966. Dynamics of flexible slender cylinders in axial flow. Part I. Theory. J. Fluid Mech. 26(4), 717-736. PdCMTM~a, R.D. and MoaTo~, K.W. 1967. Difference Methods for Initial-Value Problems. Interscience Publishers, New York. RJslqN, P. 1980a. Memorandum. David W. Taylor Naval Ship Research and Development Center. Bethesda, Maryland. R ] s ~ , P. 1980b. Memo. Hydrodynamic loading functions for a cylinder at small yaw angles. David W. Taylor Naval Ship Research and Development Center, Bethesda, Maryland. RIslqN, P. 1980c. Data package No. 1 for cable and array maneuvering. David W. Taylor Naval Ship Research and Development Center, l~ethesda, Maryland. RUSSELL, J.J. and ANDERSON,W.J. 1977. Equilibrium and stability of a circularly towed cable subject to aerodynamic drag. J. Aircraft. 14, 680-686. SANDERS,J.V. 1982. A three-dimensional dynamic analysis of a towed system. Ocean Engng 9(5), 483-499. SCHECHTER,S. 1900. Quasi-tridiagonal matrices and type-insensitive difference equations. Q. appl. Math. lg, 285-295.