Small perturbation analysis of oscillatory tow-cable motion

DEPARTMENT OF THE NAVY

N4VAL SHIP RESEARCH AND DEVELOPMENT CENTER BETHESDA, MD 20034

SMALL-PERTURBATION ANALYSIS OF OSCILLATORY TOW-CABLE MOTION

by

Keith P. Kerney

APPROVED FOR PUBLIC RELEASE; DISTRIBUTION UNLIMITED

TABLE OF CONTENTS

Page

ABSTRACT...

. . . .. .

. ...

ADMINISTRATIVE INFORMATION

... .

i . . i

INODUCTION . .

...

1

EQUATIONS AND BOUNDARY c0NDmONS FOR WO.DIMENSIONAL

CABLE MOTION ... ...2

SMALL-PERTURBATION FREQUENCY DOMAIN EQUA IONS AND

THEIR SOLUTIONS

COMPUTER PROGRAMS 18

ACKNOWLEDGMENTS

...

21

APPENDIX A - DERIVATION OF EQUATIONS OF MOTiON 23

APPENDIX B - LISTING AND SAMPLE OUTPUT OF OMWAY 25

APPENDIX C - LISTING AND SAMPLE OUTPUT OF.:ÓMFLO 35

REFERENCES . . 46

LIST OF FIGURES

Figure 1 - Quadrant I Towing Configuration 22

ABSTRACT. .

The equations of two-dimensional motion" of a' flexible inextensible cable are linearized by â small-perturbation approximation and sinusoidal

time dependence is assumed. The simplified. equations are integiatéd

nu-merically by the Kutta-Merson method. Separate computer programs,.

OMWAY and OMFLO, have been written for the Quadrant I and Quadrant

II cablé-towed-body' problems. and are listed iw the appendixes. .

ADMINISTRATIVE INFORMATION . .

-This research Was supported by the Naval 'Air Systems Cômmând and b.y the

General Hydromechanics Research Program under SR 009. 0101,. Task 0102 sponsored

by the Naval Ship Systems Command COde 034 l2B... . ..

INTRODUCTION .' . ..

AccUrate prédiction of the motion of cable.towed-body systems is essential for

their rational design and- for simulation of their operation. For the. simplest of these

systems, a single body towed. by a cable attached to a tow g ship, adequate.

represen-tations of the' towPointmotion, the dynamic behavior of the cable as it is affected

by the tensión at its two énds and by.. hydrodynamic ànd. gravitational f rces,'.

and the motion of the towed body are needed. If only steady-state motion is

considéred Whéré the éntire system is in steady rectilihear tratislation, analysis is possible. by numerical integration of the differential equations that-. describe the cable configura ation and tension. The case of steady-state. mOtion of ah inextensible toW.. cable in thé

plane described by

ait

and th dfrectiöÌi öf translation has been 'solved numerically

by Cuthill.' The tow speed and 'the tension. and angle at. the towed-body end Of the.,

cable are prescribed and theviscous force on the cablé 'is represented by functions of

cable angle based on experimental data. Integration of the two simúltànecMls first-ordér'

nonlinear ordinary differential equations fOr the cable tension and angle as functions of

distance along the cable is performéd by the Kutta-Mrson method.

The unsteady two-dimensional motion of the same system was considered by

Whicker.2 _{Instead of prescribing tension and angle at the towed-body end of the cable,.}

the equations of motion of the body become a bundary condition for the equations

differential equations for the tension, angle, and twO velocity components as fUnctions

of time and distance along the cable. Unless the restrictiOn tO an inextensible cable is

relaxed, the mathematical classification of these equations is parabolIc since two are

hyperbolic and two aré pàrabóliç.: Consequently they cannot be solved by the method

of characteristics.

The ÌIeXt section of this report describes the system of iñextensible cable equations given in Reference 2 and a particular way of representing the viscous force

on the cable. Then two importantj simplifying approximations are introduced: the

unsteady dependent variables are assumed to be small perturbations of the steady-state

variables, and their time dependence is assumed to b sinusoidal. The first of these

permits the lInearization of the unsteady equations; .the. second changes them from partial to ordinary differential equations with distance along the cable as independent

variable. The digital computer programs used to integrate the equations are described

in the next Section .and are listed in Appendixes 'B arid C.

hus the accuracy of the analysis reported here rests on the validity of four rnan. 'assumptions: (1) that the motion of a real cable can be tepresented by the solution' of the equations of motion of an inextensible cable, (2) that it remaiis in a vertical plane, (3) that it consists of small excursions from a steady-state configuration, and (4) that thése excutSions are sinusoidal in time.

EQUATIONS AND BOUNDARY CONDITIONS FOR TWO DIMENSIONAL CABLE MOTION

The two-dimensional motion of an inextensible cable is described by four

nonlinear first-order partial differential equations in which independent variables are distance along the cable s and time t and dependent variables are the cable angle Ø, tension T, and normal and tangential velocity components U and V. For the case where the motion is in a plane containing gravity and Ø is measured from the hori-:

zontal thé equations are those derived in Appendix A whiç are

au

_acb

-as as

-a

av

u +-

ös as

and

4av

._+

(1+À)U__]=____

a aT - w sin Ø

(4)

where ,t is the cable mass per unit length, pA iS its added mass (in water) per unit

length, w is its weight in water per unit length, and FL and FG are the nOrmal and

tangential components of the viscous fOrce per unit length acting on the cable. These

equations were also derived by Whicker2 fòr the case of X=O.

Two cable-towed systems are considered in this. report: in one a surface ship

in a seaway is towing a deeply submerged body, and in the other,a deeply submerged

submarine is towing a float which is slightly beneath the free surface and subject _to

disturbances from its seaway. The two configurations are, shown in Figures 1 and

2. _{In either configuration the (mean) towing direction ¡s to the right and the origin}

of coordinates is at the center of mass of the towed body with x in the (mean)

towing direction and y up. _{For the two configurations Ø lies between Q and ir/2 and}

between ir/2 and ir so they are known as the Quadrant I and Quadrant II configurations

respectively. _{The cable is of length L and the quadrant choice of Ø causes. the towpòint}

to be locateß at sL in Quadraht.

I' and sL in QuadrantI!.

More detailed discus-sion of these quadrant conventions is provided by Springston.3

In the Quadrant I problem the cable motion is. excited by the: effect of the.

seaway on the. towing ship so the boundary conditions at s=L are the, motiOn of the towpoint; the toundary conditions at the bottom are, provided by the equations of.

motion of the body. tethered by the cable. . In the Quadrant II problem the towpoint.

is on a. deeply submerged SUbmarine arid the towed body. is. subject to.

forces due to the. seaway Therefore the boundary conditions. at .s0 are . the equations

of motion of the . tethered body .subject to these forces; the other boundary. conditions

are that thç point sL be in steady rectilinear translation. _{Thus each problem has} kinematic boundary cOnditions at one end and dynamic conditions at the other, añd in

either problem the excitation is due to the action of the _{seaway; this' excitatiOn is}

through the kinematic boundary conditions in Quadrant I and through. _{the dynamic}

boundary conditions in Quadrant II.

In both problems the steady towing Speed is c. _{In the Quadrant I problem the}

(t) and Ch (t) and are assumed to be known. Therefore the kinematic boundary

con-dition for Quadrant I is

U(L,t) sin 4(L,t) + V(L t) cos Ø (L,t)=c + c (t)

and

-_U(=L,t)cos(L,t)4-V(L,t)sifl'ø L,t)=O

'The towed-body equations of motiOn used here are those bf a very sinip1

representation of the. body dynamics. The body is assumed to be a point of mass m

with added masses in surge and heave m and mh, .ard weight (in water) WB on

whiöh known drag and lift fòrces DB and L, act. Nô body 'pitching is represeñted,'

any body motion which is excitéd is undamped and the mass center of the bódy,

center of buoyancy, and point Of c'able atthchment are coincident. 'In the Quadrant I

problem: DB. and L are takeñ as constants but the Quadrant II 'problem 'they

represent the exciting force from the seaway and are timie-dependent. If Uß and

VB are the horizontal and vertical componénts of thebôdy Velocit', giveñ by

Uß(t)= U(O,t) sin Ø (O,t) + V(O,t) cos (O,t)

- U(L ,t) .cos (L,t) V(L,t) sin Ø (L1t)= Ch (t)

In the Quadrant II probleiù the töwpôint i in steady, horizontal translation at speed c

so the kinematic boundary condition is

-U( L,t) sin Ø (- Lit) +

V(L,t)cos(L,t)C

(7)

differentiation with respect to t and substitution into the body equations of motion leads tO the dynamic boundary conditions.

au

av

(in+ m) (0,t) Sin Ø (Q,t) + - (O,t) COS (O,t) '

at

+ [U(O,t) cos (0,t)-= V(O,t) sin(O,t)] --- (O,t)} bB ± T(0,t) cos'Ø (O,t)

aU 8V

(in +in17 ) .---- (O,t) cos Ø (O,t) + -i-- (0,t) sin Ø (O,t)

E U(0,t) sin Ø (O,t) + V(O,t) cs Ø (0,t)] -±. (O,t)}=LB - W ± T (O,t) sin Ø (Q,t) (10) with the signs in front of the T(o,t) chosen as plus for the Quadrant I problem and minus fOr the Quadrant II problem.

In steady-state towing-cable 'theory' the viscous' force components F1, and FG are assumed to be equal to D times known functions of Ø, where D is the drag per

unit length of the. câble when is ir/2. ':FL'/D and FG/D 'are known as the nóri1ial

and tangential loading functions. In this report. it is assumed that the viscous fOrce 'On

the cable in unsteady flow can be represented by an unsteady generalization of the steady-state force. D is expressed 'in' terms of CD, the coefficient of steady-state drag"

per unit length when the cable is normal to 'the flow, 'as D=

Ç

hCr. where is the.

water. density and h is' the cable thickness, or dimensión in the direction normal to thé'

plane in 'which the motion lies. The unsteady generalization of this expréssion for D

consists of assuming that D is given by ' . .

U2 +V2

D.p

.'

hCD

where CD remains the drag coefficient for the cable normal to a steady flow but U

'and V are instantaneous values of the velocity components. The unsteady

generali-zation of the loading functions consists of replacing the functional dependence on Ø

with the same functional dependence on Ô, defined by O=tarí' U/V since', as can be

seen in Figures 1 and 2, this is the angle between the tangent to the cable and the

steady-state values of U c sin Ø and V c cos .) Comparison of Figure 2 with

Figure 1 shows that the positive viscous force components in Quadrant II, FL and -FG, shoúld exhibit the same functional dependence on r-0, the acute angle between

the cable tangent and the flow direction, that FL and FG exhibit on O in Quadrant 1. Since the generalized loading functiOns proposed in Reference 3 are used in this report,. this requirement will be satisfied if the loading functions are given by

where plus signs are used for Quadrant I and minus signs for Quadrant II. A and B are coñstants to be determined.

An example of how A añd B may be related to experimental data is provided by Eames,4 who suggests that the viscous frce on the cable can be assumed tO consIst of a pressure drag, D( l-f) sin2 0, which acts normal to the cable tangent and a

frictional drag, Df, which acts parallel to th flow direction or horizontal for the

steady-state case. Consequently FL D(lf) sin2 O + Df sin O and FG Df cos 0, and A

and B are given by AO(A)

1-f

A2(")L

li-f

,

B'

= f, A1 f and the

other six. equal to zero. This representatio9 results in reasonable values at the nds of

the range of 0; for O = ir/2, FL is D an4 FG is zero, forO O, FL is zero and FG

is Df. Furthermore only two quantities, C1, and f, need to be determined from

ex-periments. They should be expected tO be dependent on Reynolds number.

The parameter X, the ratio of the added mass per unit length m water to mass per unit length for motion in the xy piane Will be represented in the form

ph27r

p

4i

(14)

=A"

±A1" cosO +A2 ₍₁₂₎

and

FG (I') (r') (F)

where p Vis a constant to. .be determined. Ideal-fluid theory predits v1 for bare:

round cables.

In, either quadrant the coordinates x and y of a point on the cable. are related to s and

by k cos Ø and k sin Ø so they are given by

as òs

x(s,t) cosØ (s',t) ds' 15]

and

y(s,t)

=J

sin Ø (s',t) ds' _[16]

To summarize, U, V, Ø and T are given as fünctions of s and t by solutions of Equations (1)-(4), with loading functions given by Equations [12] and [13] (wIth appropriate signs), D given by Euation E 11), X given by Equation [14], kinematic boundary conditions provided by Equations [5] V,and [6] or [7] and [8], and dynamic

boundary conditions by [9] and [la] (with appropriate signs). Then x and y can be

found using Equations E 15] and E 16].

SMALL-PERTURBATION FREQUENCY-DOMAIN EQUATIONS AND THEIR SOLUTIONS

Approximate solutions to the system of equations given in the previous sectiOn are obtained

V

for the case when the exciting disturbance - the towpoint motion in, the

Quadrant I problem and the time-varying force in the Quadrant II problem _{- is small}

and sinusoidal in time by assuming that Ø, T, U, and V are equal to the steady-state

values plus a small perturbation term which is proportional to the cosine of wt minus a

phase angle. is constant and the magnitude and phase angle Of the perturbation terms are functions of s found as solutions of eight linear ordinary differential equations

obtained from the four nonlinear partial differential equations, E 11 -[4]. Since the

solutions are the magnitudes and phases of a sinusóidal oscillation they are known as

frequency-domain solutions. V

The steady-state solutions are denoted by terms with subscript zero. For

ex-ample q(s,t) is given by- -

V'

V

çi(s,t)= Øo(s) + ØM(s) cos kot - 6(s)] +...

øo(s)₊GM(s) Ecos c,it cos (s) + sin cot sin 6(s)] +. . [17] = (5) + 1R (s) cos cot Øj(s) sin cot+. . .

where Re means that thé real part is to bé taken, i is thé irnaginat' unit, and +

indicates that OM (s) cs [wt 60(s)] is the leading term in a perturbation series

for-Ø(s,t) - Øo (s); Ø and Øj aré rélated to thé magnitude and phase OM and by

OR = OM CO5 and 0! =0M Sifl Ø so 0M

=J0R2OJ2

- Oi J R 0 is given by 01 ØR iOj. Similarly T(s,t)= To(s) TM(s) cbs [wt-5r(s)i +; =

t0(s)+Re[T1(s)e'°'t] +.

[18]

with cörresponding relatiôns between TM, r, TR,

ti, and T1.

As wäs seen in the

previous séction, U0 and V0 are given by

U0 (s) =c sin 0o(s)

and

VQ(s) = cos Ø0(s)

and satiSfy Equations [i] and [2]. It is convenient to represent U in the fotm

U(s,t)= _{c sinlØo(s) + ¿Re[0i(s)e"ft] } +}tRe[Ui(s)eit1 +:

= c sin Oo(S) C cós Oo (s) tRe[01 (s)e1"] +

e[Ui(s)et] +...

sO the éxpansión fòr U is

[19]

In the same way

V(s,t)=c cos 0(s)+ 6e [cØ1 (s) sin Ø0(s) +V1

(s)1et }

_+.

Then if U1 = UR + ¡Uj and UM and & are given by

UM =V'(cØR COS øo + UR )2 + (cj cos + Uj)2 U(s,t)= c sin_Øo(s) 6e.1 [cØ1(s) cosØ0(s) + U1

(s)]eUt

} . [21]

and

and

to give

V(s,t)= V0(s) + VM(S) COS [wt iSv(s)] +.

In the Quadrant I problem the excitation is due to speeds in surge and heave at the towpoint, c(t) and ch (t), yet it is desirable to have prescribed displacements in surge and heave, a(t) and ah (t), as inputs.

These are given by

a(t)

=aSM cos (oat _tSe)

cçbj COS øo + U1

u=tan

-CØR OSøo + UR

Equation [211 reduces to.

U(s,t) =U0(s) UM(S) COS [wt 6u(S)] +.

Similarly, V1 = VR + ¡V1 and VM and v are given by

VM

='/(

CØR Sfl øo + VR)2 + (-cØj sin øo + Vj)2

=

tan'

cØj sin øo VI CÇ5R Sifl0 + VR

and

so c is given by

where Cs!

Similarly, ch is given by

Ch (t) = tRe _(chie'

where c ch R ich!, ch R (.oah M sin 5, , and _h! =

h

M cos '5h

The excitation for the Quadrant II problem is due to time-varying drag añd lift

forces on the towed body Therefore the forces are given by

DB =DBO+DMCOS(Wt_D)

=DBO +tRe(DietWt) and

LB =LßO Wß +LM cos(cot&L)

=LBO WB +tRe(Liet) [26]

where D1 and L1 are zero for the Quadrant I problem and given by D1 DR + iDj

and L1 LR + iL7 in the Quadrant ÏI problem. tvhee DR _{DM cos} ,

D1 _{DM sin 6D LR} = LM CO5 5L and L LM Slfl L

-Substitution of Equations [211 and [221 into Equation [Il] gives P/iCD

D

= 2

c.2 _{+2c sin øo tRe [(cØ1 cosØ0 +} Ui)el(t]

+ 2c coS0 tRe [(cØ sin

+ V1)e"t]+

_} 127]

=.0

1

+ _{tRe [C1 (s)e'}

J }

+..

where D0 = (p/2)c2hC' and c1 U sin o + V1 cos Equatiòns [21] and [22] give

ah(i) a,1M COS(c)tSh)

[23]

e5(t)

= - 'sM

SIfl (wt

-wa5M (siñ o.,t cos & cos(sit sin

= CsR CÔS cot + Cs! Sifl cot

where

= tan Ø

=_{tan 00 {}

TherefOre

U_ csinQ +e [(cO1 cosØ0 +U1)ei)tJ+...

COS Oo + ¿Re [( cçb1 sin Oo + V1 )e)t] +.'..

i+ó:e(C01C0S00+Ul

ftf

-

\

csinØ,

/

i+e(00"1

_ccosØ

/

+geI(1

cosØ0+U1 cØ1 sinØ0V1

L

\

csinØ0 - ccosØ0

0tan{tanØo+tanØo4eÍ('

cosØ0+U1

cØisiñØo_Vi\.,

L CsiflØ0

änd use of the Taylor expansion for the inverse tangent gives

o = tan Ø, ¿Re

1('

COSO + Ui + ç

l+tan2O0

_L\

csinØ0

= Øo

¿Re [0 ()eCt] +...

01

tanØ0 _cO1 cosØ0+U,1

(

+ cosØ

-

sinØ0

Substitution of EquatiOns [27] and [28] into Equations [12] and [131 shows that the expansions for the viscous-force components are

FL=FL (D0,Ø0) 2FL (D0,Ø0) ¿Re

[cì(s)

e1)t] ±L.

(DO) ¿Re[O(s)e

and

FG FG (D0,Ø0) + 2 FG (D0,00) ¿Re

[cl(s)

e

The tWo Steady-state cable equations are obtained by substituting the

steady-state terms in the expansions for 0 T, U, _{V, FL, and FG given by Equations [171}

ccosØ0

r

-sinO0 V1

_\

e' ccosØ0 _J

cØ1 srnoV1)

Ccos4 dFG

-ti+

døo (D0,00) ¿Re[O1(s)e?t]

t

[29]

[30] +

through [22], [29], and [30].into Eqtïatioñs [3] and

k].

Thus

-To-1--FL (D0,0)+.wcosØ0=0

døó [31] dT ds-

-F (D0,Ø0)-wsin00

[32]

Boundary conditions at s O are fOund by substitution of theSe steady-state terms into

Equations [9] and [10]. The result is

-DßO± T0

(0)cçs0 (0)=0

LBO - WB±TO(0) sinØ0(0) =0

(o) = tan

L80 -

WB

DBO.

[34]

DB o is positive so Equation [341 shows that WB must be greater then LB for the

Quadrant problem and less than L80 for the Quadrant II problemi

Equations [31] and [32], which are Simultaneous nonliñear first-order ordinary

differential equations for Ø0 and T0, are integrated numerically on O <s < L in

Quadrant I and 0> s -L in Quadrant II, with initial conditions provided by

Equations [33] and [341. The integration is performed by the saine Kutta-Meton

subroutine, KU TMER, that is used in Reference 1. Thus the results of Reference i are

duplicated although it is necessary to write separate proams for Quadrants I and II and the loading functions are restricted tö those that can be described by Equations [12] and [131.

Substitution of Equations [17] [221, [29], and É30] into Équations [1]

-T0 (0)=v/DB02 +(LBO - WB)

[41, use of small-angle approximations to the multiple-angle formtilas, subtraction of

Eqúations [31J and [32], and cancellation of a mutual factor e..t)t result in

dU1 d0

-

_{Vl .=iØ}

ds ds U1

+-

0 _[36] ds ds and dT1 2 . dFG 138

iwji(V1 +XcØ1

sinØ0)

F(Do,o)c1

. (D0.,Ø0)01wØ1 COSØOL

ds c ...

Since Ø and T0 are given by xrimericál integration of Equations [311 nd [321 and

appear as coefficients in Equations [35] through [381, these. four linear ordinary

dif-ferential equations must be integrated numerically. This is also done by the subroutine

KUTMER.

Boundary conditions at s = O are found by substitution of Equations [17] [22]

into [9] and [10]. Thus

- iw(m + m)[Uj(0)sin øo (0) + V1 (0) CO5 o (0)] = - ± [T1(0) cosØ0 (0)

-

T0 (0)0i:(0)sin0o(0)]

and

- :w(m +mh)[ - U1 (0) cosOo (0) V1 (0) sin Øo (0)] = L1 ± [T1 (0) sinØ0 (0)

+T0(0)Ø1 (0)çsØ0(Q)]

or

T1'(0)= ± U1(0)(m mh) sin Ø(0) coØ0(0)+ V1(0)[.m + m _{cOs2Oo (0)}

[39]

+ m sin2 Ø(0)] } +D cos Oo (0) - L1 Sin Ø (0))

añd

1WM (XcØ1 cos øo + (I + X) U1 J

= -

T0 -

T1

(Do,Øo)c

ds ds c

(D,)O

wØ1

sinØ0

dØ0

[35]

and

and (0)

= T0(0)

with plus signs for QuadÌ-ant I and minus signs for Quadrant II. Recall that D.1 and

L1 are zero for the Quadrant I problem.

Boundary conditibns at s = L for the Quadrant I problem are obtained by substitution of Equations [17] and [21] - [24] into Equations [5] and [6]. Thus

U1 (L) Sfl øo (L) Chi COS Øo (L) [41]

V1(L)=ci cosO0 (L)+chl .Sin Oo (L) [42]

Substitution of Equations [171, [21], and [22] into Equations [7] and [8] shows that the boundary conditions at s = L in the Quadrant H problem, are

[43,]

V1(L)0

[44]

Satisfaction Of the two-point boundary conditions is achieved by assigning differ-ent sets of initial values tö the dependdiffer-ent variables at the lower, or passive, end of the cable and integrating Equations [35] [38] up the cable twice, once with each set of

initial values, so that two different solutions are Obtained along the cable. then the

linearity and homogeneity of the lower boundary conditions and the equations permit linear superposition of the two solutions such that the inhornögeneous boundary

condi-tions at the upper, or excited, nd are sat sfied. This is known as a "shooting" method

of solving a two-point boundary-value problem.

The two sets of solutions are called the A-mode and Bmode solutions and the solution to the full problem is given by

- U1 (0)[m .+ ,n sin20(0) + mh cOs2Ø0(0)] V1 (0)(in 'nh) [4Q]

and

CA ViA (L) + CB V1B(L) = C1 cosØ0(L)+ Chi sinØ0(L)

are satisfied. Therefore

[cSlsinØO(L)chl cos00(L)] V1B(L)[csi cosØ0(L)+ Chi sinØ0(L)] U1B(L) _49]

U1A(L) V1B(L)U1B(L) V1A(L)

[c

cosØ0(L) + chisinØo(L)] UIA (L) - E ci SinØ0(L) - chi cos00(L).] VLA(L) [50] B U1A(L)V1B(L)U1B(L)VIA(LÏ Ul=CAU1A +CBU1B [45] Vl=CA VIA +CB V18 [46] 01

EI

T1 =CAT1A +CBTIB [481

where the subscripts on U1, V1, , and T1 denote the solution mode and CA and CB

are constants to be determined below;

In the Quadrant I problem, values of U1 A (0), U1 (0), viA (0) and V1B (0) are

assigned and T1A (0), T1B(0), iA (0) and iB(0) are found from Equations [391 and.

[401 with plus signs and D1 = L1 = O Next Equations [351[38] are integrated on

O < s L so that both modal solutions are known on O s L. Then Equations [45], [46], [41], and [42] -shoW that C'A and CB must be given values such that

Finally U1, V1, Ø, and T1 are found on O s L by substitutiñg back into Equations [45][48].

The Quadrant II problem is solved in much the same way. Values are assied

to iA (-L), lB(-L), T1A(-L), and TIB(-L). Then Equations [43) and [44] show that U1A (-L),

UIB(-L), ViA (-L) and VIB(-L) are all zero (Otherwise Equations [45) and [46), with

left-hand sides equal to zero, show that C'A and GB would be determined by conditions

at s = -L). Equadtions [35] [38] are iñtegrated on -L < s O so that both modal

solutions are known on -L s 0. Equations [45] [48] and [39] and [401 with

fniius signs show that, if P and Q are defined by

P= T(0) iw U1(0)(m - mh) cosØ0(0) sin0(0)+ _V1(0)[m+m cos2Ø0(0)

+ _{mh sin2 Ø(0)] }}

and

Q=T0(O»1(0)iw1U1(0)Lm+»l sinØ0(0)+m cos20(0)]

- V1(0)(m=mh)

cos0(0) sinØo(0)}

with PA PB QA, and QB defined with appropnate modal values of T1 , U1, and

V1 on the right-hand side, the boundary conditions at s O become

CAPA. +CBPB

= -

cosØ0(0)+ L1 sinØ(Q) and

CAQA +CBQB =D1 sinØ0(0)+L1 cosØ0(0)

Therefore CA and CB are given by

[D1cosØ0(0) L1 sinØ(Ô)1 QB [D1 sinøo(0)+L1 CÒSØO(0)]PB _[51]

A --. _r4n

B -

F9

[D1 sinó(0) +L1 cosØ0(0)]

A -

[D1 cos0(0)+L1 sinØo(0)].QA

_[51

CB = =

PAQB PBQA

Finally U1, V1, , and T1 are föufld on -L s O by substituting back into Equation

The analysis is. completed by computmg the cable configuration. This is denoted

by x and y and is found from Equations [15] and [16] Equation [17] suggests that

x and y are given by expansions of the form

x(s,t) = x0(s) + Çe [xi($)et] +. [53]

and

and

y (s,t)=y0(s)-f-tRe [y1(s)A)t +. . [541

Then if X1 = XR + ix1, xç = V/xR2 + x12, and 6

tazf'--- it

is found that

XR.

X X0 + XM cos (wt

&) +..

.

Similarly Yi = YR + iYi, YM

+ yj2, and 6,, = tan'

22.

YR

Subthtution of Equations [17], [53],. and [541 into Equations [15]. and [16] and use of the smaliangle approximatións to the multiple-angle formulas give

x0

=J

cos40(s') ds' [55]

x1 (s') sin0(s') ds' [56]

Yo sin0(s')ds' [57]

fSØ() cosØ0(s'.)ds' [58]

These integrals are to be evaluated for O < s L in the Quadrant I problem and for

O > s ' -L in the Quadrant II problem. Because of. the choice of origin, the

con-figuration given is that seen by an observer moving with the, towed body.

After the complex-valued functions , Ti,' cØ1 cos øo + U1, -cØ sin Ø +

V1, x1, and y1 have been found the amplitudes and phases of the perturbation

quantities are easily found. '

the two-dimensional cable problem sinçe there is no way of assigning a meaning to the

tüne t = 0.

Consequently the phase angle of either the horizontal or vertical ex-citation is meaningless and may. be set equal to zero.

COMPUTER PROGRAMS

The digital computer programs which solve the Quadrant I and Quadrant H problems are OMWAY and OMFLO respectively, and are listed in Appendixes B and C,'

together with sample calculations. They are written in FORTRAN IV for a CDC 6700

computer.

Although Equations [35][38]., [531, [541, and, the boundary conditions are written for complex-valued functions of s, they must be separated into real and

imaginary parts because the integrating subroutine KUThIER is not available in complex

arithmetic. However, the existence of a complex-conjugate relationship between the

real and im aginary' parts makes solution for two instead of 'fôur modes sufficient

Whenever possible the notâtiòn of the previous .section is retained; some important ex-.

ceptions are the angle which is replaced by PH, and the frequency w which is

repláced by 0M. In generai, quantities whiCh pertain tó the towed body start with aB.

Thus m, the body added mass in surge, is BMS. The cable is of thickness CH and is

N-1 segrn ents long; each segment is of length DS. Its mass per unit length .i is

replacéd by ULM and its added mass parameter y by AMP. Its weight per unit length

in water is' WUL and the density of water is DN. Further explanation of the prOgram

is provided by comment cards' in the listing.

Each program includes four subroutines. A TA prevents the inverse-tangent

function from dividing by a prohibitively small number if an angle is very close to

± ir/2. K UTMER performs integrations with a fourth-order Kutta-Merson method. This

is a generalization of the Runge-Kutta method that provides an automatic reduction in

integration-step size When an errOr ctiterion is not' met. DAUX and DARN provide

integrands för KUTMER; DAUX for , T0, x0, and Yo and DARN for U1, V1, ,

and T1. KUTMER is not used fór x1 and Yi; they are' calculated by a = simple

first-order integration Wthin the main program.

In the prOgram OMWAY, J is an integer which increâses 'from J i at the towed

body (s=0) 'to J=N at the towpoint (s=L). Successive 'data cards are used to 'Specify;

I. Name of program, OMWAY. . ' '

Body lift BL, drag BD, and weight BW.

. '' «

'

and frequency OM . :

4.Bódy. mass and added masses ¡n surge añd heave, BM, BMS, and BMH.

WUL, ULM, DS, and N. .

Tow speed G DN, Gli, cable drag coefficient 'GD, and AMP.: Constants for Equation [12] to give FL/DO.

Constants for Equation E 13] to give FG/DO.

These data are written out. ,D0 is computed and Equation [14] is used to compute X, which is represented by AMC/ULM. Equations [33,] and [34] are used to find T0

and. Øo at J=1. 'and U0, Vo, o. 7'0, x0, and Yo are computed by using Equations [1.2],

[13], [19], [20], [31,], [32], [55], and [.57] and afe written out as functions off. Values are assigned at J=1 as UJA real. and VIA zero, ØJA 'and TJA are found from Equations [39] and. [40], and the A-mode solution is computed by integrating

Equations [35] [38] on

1 < J

N. Then UIB is set zero and V1B is set real at

J1 and the B-mode solution is found in the same way.. Next c51 and Chi are

com-puted and' the required values .of U and V at JN are found 'from Equâtions [41]

and [42] and cA and B are found from Equations 149] and [50]. ' U1,' V1, Ø and

T1 are found along the cable from Equations [.45]'[48] and their magnitudes and

phases are' computed and written out. Finally. Equations [56] and [58] arC used to'

find x and Yi and their magnitudes and phases are computed and written out

OMFLO uses the nomenclature of OMWAY wherever possible.. J is an .integer

which decreases from J-1 at the .towed body (s0) to J-N. at the towpoint (s=-L).

For performing integrations down the cable, which is done when . Te,. X0, Yo, xj,

and yj are computed, an integer K is used which is defmed . by K'-J and.thus

in-creases from K1 at . the towed body to KN at the towpoint. Since these inte-grations are 'made. in the direction of decreasing s, the integrands computed ' in DAUX

are negatives of those in OMWAY. U1, V1,. .Ø, and T1 thC computed' by integration

in the direction of increasing s so an ' integer I is defined by' l=N+i-K which increases

from 11 at the towpoint to 1N at .t,he towed body. The first three data cards

differ fröm those of OMWAY and specify: ' ' .

Name .of program, QMFLO. ' ' '

Steady!state body lift and drag BLO and BDO, and body' weight BW. 3, Magnitude' nd phase' of oscillating lift and drag'forces on the body, BLM,

DBL', DBM DBD, and frequency OM '

To and

o at

K1 and Uo, V0, _o, T0, x0, and Yo are computed by using Equations

[12], [13], [19], [20], [31],

[32],.[551, and [571 and are written out as functions

of J. '1.o and To are relisted as functions of I, ØJA is Set real and TIA UJA and VIA

are set equal to zero at 11, and the A-mode solution. is corn puted by integrating Equations [351-438] on I <IN. Then TIB is set real and GIB' jJB, and VIB are are set equal to zero at 11 and the B-mode solution is computed in the same way.

Next Øi, T0, UJA,. ViA, 1Ai 7'JA» VJB, and TJB are re-listed as

functions of K instead of I. The required values, of P and Q at K1 are found from

Equations [391 and [40] and CA and CB are found from Equations [51] and [52].

Uj, V1,

j,

and T1 are found along the cable from Equations [45] [48] and their

magnitudes and phases are computed and written out. Finally Equations [56] and

[58] are used to find x1 and yj, and their magiitudes and phases are computed and written, out.

Representative times for a run involving a cable nine integration steps long (of ten feet each) are 28 seconds compilation time and 14 seconds computation time on a CDC 6700 computer.

Certain modifications to the programs are easily made. For example, the

programs listed require angles in radians and any consistent mass-length-time units for

dimensional quantities. They have been adapted to treat, as input and printout

quantities, angles in degrees, masses in pounds, lengths iii feet, and speeds in knots.

They have also been adapted tO. examine behavior over a range of frequencies, by repeating the unsteady part of the program in a DO loop, with lowest, highest, and incremental frequencies specified on an additional data card.

A computer experiment was performed On the Quadrant I program, OMWAY, by examining the frequency range 0.01 to 0.80 hertz in steps of 0.01 hertz fôr

cable lengths of 200, 400, 600, 800, 1000, and 1200 feet at tow speeds of 6, 10, and

14 knots. As criteria for successfúl. performance it was required that the perturbation

quantities must remain below the steady-state. quantities. It was found. that VM, the

magnitude of the tangential velocity, at (or very near) the towpoint was the quantIty which failed at the lowest frequency. The only cable length where its behavior was satisfactory throughout the frequency and speed. ranges was 200 feet; all the longer

ones failed but showed better behavior with increasing speed. The 600-foot cable

failed at 0.26, Ó.34, and 0.48 hêrtz at 6, lO, an4 14 knots While the 1200-foot

cable failed at 0.13, 0.20, and 0.30 hertz at the same speeds. Values for the 800-and

ERRATA for

Report 3430

Page 21 The paragraph under Acknówledgments reads:

The assistance of Dr. Elizabeth Cuthill, Dr. J. Nicholas Newman, Mr. William

E. Smith, Mr. George B. Springsn, Jr., and Dr. Henry T. Wang in conducting the research

reported here and in writing the computer programs are gratefully acknowledged. Dr. Newman deserves special. thanks for suggesting the small-perturbation frequency-approach.

It should read:

The assistance of Dr. Elizabeth Cuthill, Dr. J. Nicholas Newman, Mr. William E. Smith, Mr. George B. Springston, Jr., and Dr. Henry T. Wang in conducting the research reported here and in writing the computer programs is gratefully acknowledged. Dr. Newman deserves special thanks for suggesting the small-perturbation frequency-domain approach.

speed. The 400-foot ca1e failed at 0.49 and 0.57 hertz at 6 and 10 knots and

behaved satisfactorily at .14 knots.

The next most sensitive quantity was the perturbation tension at or very near

the ship. It behaved well throughout the frequency and speed ranges for cable

lengths of 800 feet and less. For the 1200 foot cable it failed at 0.46, 0.46, and

0.61 hertz at 6, 10, and 14 knots.

The qualitative conclusions are that the small-perturbation frequency-domain analysis is valid for short cables undergoing low-frequency oscilations and that

increas-ing towspeed can have a stabilizincreas-ing effect. Apparently as the forcing frequency

ap-proached 0.80 hertz, a resonant region was being apap-proached; frequencies high enough to be past such a region, because the inertia of the cable and body prevent excitation of their motion, might be beyond the range of practical interest.

ACKNOWLEDGMENTS

The assistance of Dr. Elizabeth Cuthill, Dr. J. Nicholas Newman, Mr. William E. Smith, Mr. George B. Springston, . Jr., and Dr. Henry T. Wang in conducting the research reported here and in writing the computer program s are gratefully acknowledged. Dr. Newman deserves special thanks for suggesting the small-perturbation

y

S=o

TOWED BODY

Figure 1 Quadrant I To ng Configuration CAB L E S= L

-TOWING SHIP TOWING SHIP

APPENDIX A

DERIVATION OF EQUATIONS OF MOTION.

Equations [1] [41, the equations of cable motion, can be derived by considering a

small element of the cable shown ¡n Fire or Figure .2. If x and y represent the coordinates of a point on the cable then x and y are functions of s and t which

satisfy ax = cos 0 as ay = sin Ø as 'òx ay smØ - cos = U

ay..

cosØ + -- sinØ = V

Take the derivative of Equation [A3] with respect to s, interchange the order of the s- and t-differentiations of x and y, and substitute from Equations [Al] and [A2] to

obtain

au

a a a

a.

ay.

a

cos0)smØ+

cosØ(sinØ)cosØ+ sinØ

which becomes Equation [1] after substitution of Equation [A41. In the saine way,

take the derivative of Equation [A4] with respect to s, ¡ntérchnge the order of differentiations, and substitute to obtain

av

a ax a a ay _ao

- = - (cos0)cosø

sinO - +

(sin0) sinØ+ cosØ

-This, with. the substitution of Equation [A3], is Equation [2].,

The apparent-momentum vector' of a segment of length s has components

,.6s(l+X)U, ¡isV, O in a right-handed orthogonal coordinate system With axes in the

directions of U, V, and upwards from the paper of Figure 1 or Figure 2. Here

u is the mass of the cable per unit length and ¡iX is the added mass in water for [Al]

motion in the direction of U. The angular velocity has components 0,0, aØ/at in the fixed coordinate system defined by the x and y axes and the direction up out of the paper in either figure. Therefore the rate of change of apparent momentum observed from the fixed coordinate system in either figure has components along the directions of LI and V equal

to and

au

p.s(l+X)

p&sV-at at E au

(1 +X) - - V

at

atJ

and

Therefore the dynamic equations of motion of the cable segment are

= resultant force on & in the direction of U [AS]

¡is

+ (I +X) U = resultant force on 6s in the direction of V [A6]

The right-hand sides of Equations [AS] and [A6] are easily computed and the results are, after division by 6s, Equations [31 and [4].

15

APPENDIX B

LISTING AND SAMPLE OUTPUT OF OMWAV

PROGRAM ONWAY ( INPUT,OUTPUT,TAPEb=INPUTiTAPE6OUTPUT)

OIMENSION PI1O(4OO) TO(4OO) U0(400). VO(400), PHRA(400).

PHRB(400), PHIA(400). PHIB(4001, TRA(400), TRB(400). TIA(400), TIB(400). URA(400). URB(400). UIA(400). UIBI400).VRA(400). VRB(400). VIA(400). VIB(400),PHR(400), P111(400). TR(400),

TI(400), UR(400)9 UI(400). VN(400), VI(400), FL(400), FG(400),

+ DFL(400), DFG(400), T(8)

COMMON DS,WUL,O,ALANO,ALAM1,ALAM2,BLAM1,

1BLAM2,AGAMO,AGAM1 !AGAM2,BGAN1 ,B(IAM2.

10 20M,ULM,APP,C,PHO,DFL,DFG,J9TO,PIIR9PHI,FL,FG

EXTERNAL DAUXDARÑ READ (5,11) TITLES

BLcBDBW'

4 ASH, DASS ANN, DAN, 0M,

BM, BMS, BMH,

+ WUL.. ULM, DS, N, C, ON, CH CD, AMP, ..

+ ALAMO, ALAM1. ALAM2. BLAMI. BLAM2

AGAMO, AGAM1, AGAM2, RGAM1, BGAM2

20 C WRITE TITLE AND INPUT DATA

WRITE (6,21) TITLE

WRITE (6,22) BL , BD , BW

WRITE (6,23) ASH, DAS, ANN, DAN, QN WRITE (6,24) BM BMS, 8H11

25 . WRITE (6,2b) WUL, ULM. .05, N WRITE (6,26) C, DN, CH,. CD, AMP

WRITE (6,27) ALAMO, ALAM1, ALAM, BLAM1,BLAM2 WRItE (6,28) AGAMO, AGAM1, AGAM,BGAM1, BGAM2

C INITIALIZE 30 D= Q.5*DN*C*C*CH*CD ANC 3.14Ï593*CHCt4DNAMP/4.O APP = 1.0 AMC/ULM CXO = 0.0 CYO = 0.0 35 CXR 0.0

CI

0.0

CYR0.0

Cyl = 0.0 JA = N-1

40 C SOLVE STEADYSTATE PROBLEM

ERR=.00001 ERA.001 CALL ATA(8o,Y,c) P1108 = ATAN2(BWBL,YIÇ) : 45 TOB = SQRT(8DBD.(BWBL)(BW-8L)) WRITE (6,31) WRITE (6,3e) P110(1) = PHOB

T0(i)T08

50 C MULTIPLY Tu FIND VELOCITY COMPONENTS

DO 101 J 1,N

U0(J) = CSIN(Pl1O(J))

VO(J) = C*CUS(Pi1O(J))

WRITE (6,33) P110(J), TO(J), UO(J)c VO(J), CXO,CYO. J

c INtEGRATE TO FIND TENSIONI ANGLE, ANO CABLE CONFIGURATION T(1)=TOLJ) T()PHO(J) 1(3) CX0 60 T(4) = CYO FIRST.0 ALS=FLOAT (J) *s CALL KUTMER(4IALS,T,ERR,DS,FIRSI,HCX,ERA9DAUX) TO(J+1)T(1) 65 PMO(J.1)=T(2) CXÔ = T(3) CYO 1(4) 101 CONTINUF.

C SOLVE DYNAMIC PROBLEM

70 WRITE (6,34)

WRITE (6,3e)

C COMPUTE VISCOUS-FORCE TERMS ALONG CABLE

DO 201 J = 1, JA FL(J) = D*(+ALAM0.ALAM1*COS(PH0(1!)ALAM2*C0S(2.O*0(J 75 BLÁM1*SIN(PHO(J))+BLAM2*SIN(2.U*PH0(J))) FG(J) = D*( AGAMO+AGAM1*COS(PHO(J)).AGAM2*COS(2.0*PH0(J)) BGAM1*SIN(PHO(J))+I3GAM2*SIN(2.0*PH0(J))) DFL(J) = O*(_ALAM1*SIN(PH0(J))_.0*ALAM2*SIN(2.0*PHO(J)) BLAM1*COS(PHO(J)),2.0*BLA142*CUS(2.0*PHO(J))) 80 DFG(J) = O*(_AGAM1*SIN(PH0(J)).0*AGAM2*SIN(2.0*PH0(J)) BGAM1*COS(PHO(J) )+2.0*BGAM2*CUS(2.0*PHO(J))) 201 CONTINUE

C COMPUTE MODAL SOLUTIONS

C COMPUTE A -MODE SOLUTION

85 C ASSIGN VALUES TO VELOCITY COMPONENTS AT BOTTOM

ORA(i) 0.000001

tuA(i) 0.0

VRAU) 0.0

VIA(i) = 0.0

90 C COMPUTE CORRESPONDING VALUES OF ANGLE AND TENSION AT BOTTOM

PHRA(1) OM*È_UIA(1)*(BM,BMS*SIN(PHO(1))*SIN(pH0(1))

BMH*COS(PH0(1))*COS(PN0(1))),V1A(1)(BMHM*SjN(O*CO

PHO(1)))/10(1)

PHIA(1) = _OM*(.URA(1)*(BM,BMS*SIN(P110(1) )5IN(PHO(1) )4

95 BMH*COS(PH0(1) )COS(PHO(1)))+VRA(1)*(BMH_BMS)*SIÑ(PH0(1) )COS(

+ PI-I0(Ï)))/T0(1)

TRA(i) = OM*(UIA()*(BMH_MS)*SIN(Ph0(1))*C0S(PHO(i

VIA(I)*(BM+BMS!COS(PH0t1))*COS(PH0(1)

BMH*SIN(PHO(i))*SIÑ(PKO(1))))

100 hA(i) = _OM*(_URA(1)*(BMHBMS)*SIN(PHO(1) )COS(PH0(1))

VRA(i)*(B14+BMS*CQS(PHO(1))*C0S(PH0(1)

)

+ BMH*SIN(PHO (1) )*SIN(PHO(i))))

C FIND VELOCITY COMPONENTS' ANGLE, AND TENSION BY INTEGRATING UP

C THE CABLE 105 DO 202 J = 1 JA T(1)PHRA(J) -T(?)PHIA(J) H t(3FTRA(J) . T(4)TIA.(J) 110 T(5)URA(J)

1.15 120 125 130 135 140 145 150 155 160 165 T(6)UIA(J) T(7)VRA(J) T(B)VIA(J) . . . FIRST.0 . ALS=FLOAT(J)DS . . CALL KUTMER(8,ALS,T,ERR,DScFIRSI,HCXcERA,DARN) PHRA(J.1)T(1) . PI-IIA(J+1)1(2) TRA(J+1)T(3) TIA(J1)=T(4) . URA(Ji)=T(5) . . UIAiJ+1)=T(6) VRA(J+1)=T(7) . . . . VIA iJ+1) (8) 202 CONTINUE

C COMPUTE. B MODE SOLUTION

C ASSIGN VALUESTO VELOCITY COMPONENT'S AT BOTTOM

uRB(i) = 0.0 UIB(l) = 0.0

VRB,(1) = 0.000001

VIB(i) = 0.0

C COMPUTE CORRESPONDING VALUES OF ANGLE AND TENSION At BOTTOÑ

PHRB1) = OM*(_UIB(1)*(BM.bMSDSIN(PHO(1I)*SIN(PHO(1)i+

BMH'COS(PHO(l) )COS(PH0(1) ) )sVIB(1)*(BMHeBMS)*SIN(PHO(1))*COS( PHO(1)))/10(l)

PiIB(1) = _OM*(URB(1)*(BM,8MS*SIN(PHO(1) )SIN(PHO(1)).

BMH*COS(PHO(1))*COS(PHO(1)))+VNB(1)'(BMH.BMS)*SIN(PHO(l))*COS( PHO(1)))/10(1) . TRB(1) = OM*(_UIB(.1)*(BMH_BMS)*SIN(PHO(1) )COS(PH0(1)) + VIB(1)*(HMBMS*COS(P,10(.1))*COS(PH0(1))I BMH*SIN(PHO(1))*SIN(PHO(1)))) TIB(1) '_OM*(URB(1)*(BMHBMS)*SIN(PH0(i))*COS(PHO(1)) VRB(1)*(BM.HMS*COS(PHO(1))*COS(PHO(1))+ BMH*SIN(PHO(1) )SIN(PH0(1)')'))

C FIND VELOCITY COMPONENTS, ANGLE, AND TENSION BY INTEGRATING UP

C THE. CABLE , . . . DO 203 JA' J(1)PHRB(J) T(2)PHIB(J) , T(3).TRB,(J) . T(4)t18(J) T(5)URB(J) . . T(6)U18(J) . T(7)VRH(J) . , 1(8) VI3 (j) FIRST.0 . . ALSFLOAT (J)*D5 . CALL KUTMER(8,ALS,T,ERR,DS,FINST,IICX,ERA,DARN) P'IRB(J.l)T'(l) . .. PHIB(J1)T(2) TRB(J+1)T(3) TIB(J1)T(4) . . ..' . .. URB(J+1)T(5) . UIB(J+1)t(6) : VRB(J+1)T(7) . .

VIB(J+1)=1 (6)

203 CONTINUE

C SUPERIMPOSE MOOAL SOLUTIONS

C COMPUTE VELOCITY COMPONENTS FORCEL) AT TOP

CSR=OM*ASM*SIN (DAS)

CSIOM*ASM*COS (DAS)

CHRQM*AHM*SIN (DAH)

CHI=OM*AI4M*COS (DAH)

UR(N) = CSR*SIN(PHO(N)) - Ci1RCOS(PH0(N))

UI(N) = CSI*SIN(PIiO(Ñ)) - CHIÓCOS(PHO(Ñ))

VR(N) = CSR*COS(PHO(N)) CHRSIN(PH0(N))

VI(N) CSI*COS(PHO(N)) CMISIN(PHo(N))

C COMPUTE CONS1ANTS NEEDED TO SUPERIMPOSE MODAL SOLUTIONS

DENOMP = URA(N)*VRB(N) - UIÀ(N)*VIB(N) - URB(N)*VRA(N)

+ UIB(N)*VIA(N)

DENOMI = UIA(N)*VWB(N) URA(N)VIB(N) UIB(N)*VRA(N)

-URB(N)VIA(N)

A1UMR = UR (N)°VRB(N) UI (Ñ)VIB(N) - URB(N)°VR (N)

UIB(N)*VI (N)

AJUMI UI (N)VRB(N) UR (N)VIB(N) UIB(N)VR (N)

-+ URB(N)*VI (N)

BNUMR = URA(N)VR (N) - UIA(N)*VI (N) - UR (N)aVRAjN)

+ UI (N)*VIA(N)

BNUMI = UIA(N)*VR (N) URA(N)*VI (N) - UI (N)*VRA(N)

+ UP (N)*VIA(N)

DENOM2 = _{DENOMR*DENOMR} OENOMIOENOMI

CAR = (ANUMR*DENOMR ANUMI*DLNOMI)/DENOM2:

CAl = (_ANUM*DENOMI+ ANYMI*DENOMR)/DENOM2

CBR (BNUMR*OENOMR + BNUMI*DÈNOP4I)/DENOM2

CBI = (_BNUMR*ONOMI BNUMIDENOMR)/DENOM2

C COMPUTE VELOCITY COMPONÈNTS ANGLE, AND TENSION ON CABLE

DO 204 J 1,N

UR(J) = CARURA(J) - CAIUIA(J) + CBR*URB(J). CBIUIB(J)

UI(J) CAI*URA(J) + CARUIA(J) CBI*URB(J) + CBR*UIB(J)

VR(J) = _{CAR*VRA(J) - CAI*VIA(J)} CBR*VRB(J) - CBI*VIB(J)

VI(J) = CAI*VHACJ) CARVIA(J) CBIVRB(J) CBI*VIB(J)

PHR(J) = CAR'PHRA(J) - CAI*PHLA(J) CBRPHR$(.J) CBIPHIB(J)

PHI(J) = CAI*PHRA(J) + CAR*PHIA(J) CBIPHRB(J) CBR*PHIB(J)

ÏR(J) CARÒTRA(J) - CAI*TIA(J) CBRTRB(J) - CBI*TIB(J)

TI(J) = CAITRA(J) + CARTIA(J) CBITRB(J) CBRTIB(J)

C COMPUTE PHYSICAL VELOCITY COMPONENISON CABLE

URP = _{C*PHR(J)*COS(PHO(J))} UR(J)

UIP = _{C*PHI(J)*COS(PHO(J))} UI(J)

VRP = C*PHR(J)*SIN(PÑO(J)) VR(J)

VIP = - C*PHI(J)*SIN(PHO(J)) +. VI(J)

C COMPUTE MAGNITUDE AN PHASE OF VELOCITY COMPONENTS ANGLE, AND

C TENSION ON CABLE

UM SORT (URP*URP UIP*UIP)

CALL ATA(URP YK)

DU = AtAN2(UIP, YK)

VM = SQRT(VRP*VRP + VIP*VIP)

CALL ATA(VNP, YK)

DV ATAN2(VIP, YIc)

PHIl SQRT(PHR(J)PHR(J) PHI(J)PHI(J))

OPtI ATAN2(PIII(J), .,i

TM = SORT(TR(J)TR(J) TI1J)TI(J)) CALL ATA(TR(J). Yk)

DT = ATAN2(TI(J) YI()

225 WRITE (6,36) PHM. DPH. TM. DT. UM, DU. VM. 0V, J

204 CONTINUE

WRITE, (6,37)

C COMPUTE MAGNITUDE AND PHASE OF DYNAMIC CABLE CONFIGURATION

00205 J=l!N

230 C*M = SQRT(CX.R*CxR,CxI*CxI)

CALL ATA(CXR. YK)

OCX AIAN2(CXI. YK)

CYM = SQRT(CYRCYR.CYICYI)

CALL ATA(CYR. (k)

235 DCY ATAN2(CYI. Yk)

WRITE(638) CXM. DCX. CYM. DCV. J

IF (J.EQ. N) GO TO 205

CXR - PHR(J)SIN(PHO(J))DS CXI - PHI (J)°SIN(PHO(J))°DS 240 CYR = PHR(J)COS(PHO(J))0S CVI PHÏ(J)COS(PMO(J)'L)S 205 CONTINUE

11 FORMAT t 10*. A6 / lOX, 3Flu.5/ lox, 5F10.5/ lOX, 3F10.5/

lox. 3Fb S. Ii0/1OX.bFbO.5 / lOX. 5Fl0 5/ lOX. 5Fl0.5)

245 21 FORMAT (1111. 10*' 57HSURFACE SHIP IN A SEAWAY TOWING A DEEPLY

SUBMENGED WEIGHT/ 1PI-. ¡Ox. A6///)

-22 FORMAT t ix. 8HBL = , FlO.b. 4*. 81180 F10.5. 4X.

BHBW FlO.5 II)

23 FORMAT iX. 8HASM = F1O.5, 4*, 8HDAS = . F1O.5. 4*, 250 BIIAHM FlO.5. 4*. 8HDAH = , Fl0.5. 6X 8110M =

FlO.5//)

24 FORMAT t ix. 8118M = , FbO.5. 4x. 8HBMS. =' FlO.5, 4x. BIIBI4H ,. FlOaS /1)

25 FORMAT t 1*. 8HWUL = , F10 5. 4*, 8HULM = . Fb0.5. 4*,

255 BHDS .= FlO.5. 4*, SPIN = , ¡10 1/)

26 FORMAT t IX. 8HÇ = . F1O.5, 4X. 8MDN = , FlOeS, 4*,

+811CM .F1O.5, 4x, 8HCD = ,FlO.5. 4X, 8HAMP = ,F1O.5//) 27 FORMAT t iX, BIIALAMO = .FiO.5. 4*. 8IIALAMI . , FiO.5,4*.

BHALAM2 = . FbO.5. 4x, 8HBLAM1 = . FlO.5 4x. 8HBLAM2 =

260 FiO.5 //) . .

V 28 FORMAT (IX. 8HAGANO = , FiO.5, ax. 8HAGAM1 . , F1Os59.4X,

BHAGAM2= Fl0.5, 4*. 8HBGAM1 = ' flO.5. -4X. 8HBGAM2 =

F10.5 //) V

31 FORMAT. t 22111STEADY-STATE SOLUTiON) .

265 32. FORMAT t 1MO, 1A.6HPHO(J), lox, 6HTOIJ) . lox. 61IUO(J)- , lox. 611V0(J) lOX. 6MCAOtJ), iox. bIICYO(J). lOX, 1HJ

//)

33 FORMAT t 2*. F9.6.. 7x. F13.6 3*. 2( F1l.6. 5x). 2tFl3.6. 3X). 14) 34 FORMAT t 17H1DYNAM1C SOLUTION)

35 FORMAT t 1110. ix. 6HPPIM(J). 9*. 6HDPH(J). 9x, 6HTM(J) . 9x.

270 + 6HDT(J) 9*, 6PIUM(J) 9*, 6HUU(J) 9*. 6PIVM(J) 9x. 6HDV(J) ,

9, 1H..) //) .

36 FORMAT t 2X. 2( F9.6, 6*). F12 6. 3*. F9 6. 6x. 2(FlO.6. 5*.

F9.6, 6*). 14)

37 FORMAT t lHl 1*. 6l1cXM(J).10x 6HDCX(J). lox. 6HCYM(J). lox.

275 6HDCY(J). lox. DU 1/) -, 38 FORMAT t x, 2( F12.6. 4*. F9.6. 7x). 14) STOP . . V END V

SUBROUTINI ATA(X,YIc)

THIS SUBROUIINE PREVENTS ATAN2 FROM DIVIDING BY ZERO WHEN

COÑPUTING TH ARCTANGENTS OF PI/E ANO 3*PI/2

ER = 0.0000.1 IF (ABSIX).GT.Ek) GO TO lo YP( = ER RETURN 10

YKX

RETURN END SUBROUTINE

C THIS SUBROUTINE PERFORMS FOURTHORDER KUTTAMERSON INTEGRATIONS

DIMENSION '(0(10) ,Y1(l0),Y2(lO)00(10) ,Fl (10) ,F2(10) IF (FIRST$20, 10,20 5 10 HÇH .KUTMOO5O IPLOC1 KUTMOO6O FIRST1. KUTNOO7O 20 LOCO . KUTMO080 HCXHC - KUTÑOO9O 10 30 CALL DAUX(1,YO,FO) KUTMO100 39 03 40 I1N KUTMO11O 40 Y1(I)=YO(I)(HC/3 )*FO(I) KUTMO12O

CALL DAUX(T.HC/3.,Y1,FI) KUTMO13O

DP 50 I=l,N KUTMO14O 15 50 Vi (1)O(l)+(HC/6.)FP(I),6*FlU) KUTNO15O CALL DAUX(1.HC/3..Y1,Fl) . KUTMO16O 00 60 11,N KUTMO17O 60 Yl(I)YO(I),HC/8.*FP(1),.375**F( . KUTMO18O

CALL DAUX(T+HC/2.,Y1,F2) KUTMO19O

20 DO 70 I=1.N

KUTMO200

70 Y1(I)Y0(I).HC/2 *F0(I)-1 5*HC*P1(I)+2 *HC*F2(I) KUTMO21O

CALL DAUX(T.HC,Y1,Fl) KUTMO22O

DO BO I1,Ñ KUT140230 80 y2(I)YO(I)+MC/6..FO(I)..666666b7* F2(1 ThC/6.)°F1 (I) KUTMO24O 25 INC=O KÚm0250 DO 110 I1'N KUTMO26O ZZZ=ABS(Y1U))A KUTMO27O IF(ZZZ) 85,87,87 . KUTMOZ8O 85 ERRORABS(.2*1(PY2(IW KUTMO29O 30 F(ERROR_Á)100,100,90 KUTMO300

87 ERRORABS( 2 2*Y2(I)/Y1(I)) KUTMO31O

IF(ERROREPS) 100,100,90 KUTMO32O 9 X128.*ABS(HC)ABSIH) KUTMO33O IFX) 91,95,95 KÜTHO340 35 91 WRIT(6'92) T,ERROR. KUTMO3SO

92 FORMAT(2111 RELATIVE ERROR AT 7= 1P1E12.3i3HIS F10.6. ) KUTMO36O

FIRST = 2. KUTMO37O RETURN KUTMO38O 95. HCHC/2. KUTMO39O 40 IPLOC=2 °IPLOC KUTMO400 L0C2 *LOC KUTMO41O HCXHC KUTMO42O GO TO 30 . KUTMO43O 100 IF(ERROR*64.EPS)11011O'l°' KIJTMO44O 45 101 INC1 KUTMO4SO 110 CONTINUE KUTMO46O 111 T=T.HC KUTMO47O DO 112 I1,N KUTMO48O 112 Y0(I)Y2(I) KUTMO49O 50 LOCLOC+1 KUTMOSOO IF(LOCIPLOC) 120,210,210 KUTMOS1O 120 IF(INC)210,130'210 KUTMO52O 13OE IF(LOC_(LOC/2)*2)210,1402lO KUTMO53O 140 IF(IPLOC1)'210.21°'200 KUTMOS4O KUTMO55O 5 lo

5 10 15 20 25 5 10 15 20 25 30 SUBROUTINE DAUX(Z,TX,F)

THIS SUBROUTINE PROVIDES THE INTEGRANDS FOR KUTHER WHEN IT

COMPUTES THE STEADY-STATE TENSION. ANGLE, AND CABLE CONFIGURATION

DIMENSION TX(S), F(8),PHO(400),DFL(400),DFG(400),T0(400),

PHR(200). PHI(200)

COMMON DS,WUL,D,ALAMO,ALANI.ALAM2,BLAM1,

1BLAM2,AGAMO ,AGAM1 ,AGAM2,BGAM1 .BbAM2,

20M,ULM,APP,C,PHO,DFL,DFG,J,T0,PHR,PHI TX1TX(1) 1X2TX(2) C1COS(TX2) S1SIN(TX2) C2C0S (2.*TX2) S2SIN(2.1X2).

FL=D* (ALAMO +ALÄM1C1 .ALAM2*C2.BLAM1*S1.BLAÑ2*S2)

FG D*(AGAMO.AGAM1*C1,AGAP42*C2.LSGAM1*S1,BGAM2*S2) JF(TXI.EQ.O.) GO TO 50 (WUL*C1FL)/TX1 = Cl = Si 7 F(1)= (WUL*S1,FG) RETURN 50 F(2)0.O GO TO 7 END SUBROUTINL DARN(Z,T,F)

C THIS SUBROUTINE PROVIDES THE INTEGRANDS FOR KUTMER WHEN _IT

C COMPUTES TH DYNAMIC ANGLE, TENSION. AND VELOCITY COMPONENTS

DIMENSION T(B), F(8),PHO(400),UFL(400),DFG(400).T0(400)c + FL(400),FG(400),PHR(400),PHI(400) COMMON DS,WUL,D,ALAMO,ALAM1,ALAM2,BLAM1, 1BLAM2,AGAMO,AGAM1,AGAM2,BGAM1,BUAM2, 20M,ULM.APP,C,PHO,DFL,DFG,J.TO,PHR,PHI ,FL,FG CI=COS(PHO(J)) S1=SIN(PHOIJ)) IF(TO(J).EQ.0.0) GO TO 50 F()(_OM*ULM*(APP*T(6)+(APP1.U)*C*T(2)*C1)T(3)*(PHO(J+1)PH0(J) +)/DS_(T(1)+iT(5)*C1_T(7)*S1)/C)*DFL(J) + _{_2.0*FL(J)*(T(5)*Sl+T(7)*C1)/CWUL*T(i)*S1)/T0(J)} F(2)=( OM*IJLM*(APP*T(5) +(App-1.0)*C*T(1)*C1)T(4)*(pHO(J+1)-pHO(J) ),DS_(T(2)+(T(6)*C1_i(8)*S1)/C)*OFL(J). 2.OáFL (J)(T (b)*S1.T (8)*C1)/CWUL*T(2)eS1)/TO (J) 7 F (3)=( OH*ULM*(T () +iAPP-l0)0C0T(2)S1) (1(1) .(T (5)C1-T(7)*S1)/C)*DFG(J)+2.0*FG(J)* (T(5)*S1+T(7)ØC1)/C+WUL.*T(1)*C1) F (4)=(-OMULM(T (7) (APP-l.0)C'T (1 )S1) (T(2) .(T(6)'Cl-T(8)*Si)/C)*DFG(J)+2.0*FG(J)* (T(6)*S1.T(8)*C1)/C+WUL*T(2)*CI) F(5)_OM*T(2),T(7)*(PHO(J+1)PHú(J) )/DS F(6)= OM*T(1)+T(8)*(PP$O(J+i)PHÚ(J) )/DS F(7)-T(5)*(PI1O(J+i)-PHO(J) )/DS F(8)-T(6)(PHO(J+1)-PH0(J) )/OS RETURN 50 Fil)0.O F(2)0.0 GO TO 7 -END LOCLOC/2 KUTMO56O IPLOCIPLOC/2 ICUTNO57O 210 IFIIPL.0C-L0C)30.22030 KUTNO58O 220 RETURN KUTMOS9O 60 END KUTMO600

SURFACE SHIP IN A SEAWAY TOWING A DEEPLY SUBMENGED WEIGHT 0M WAY = 0.00000 BD s 859.00000 BW 6.040.00000 = 0.00000 DAS = 0.00000 AHM; = 1.00000 DAM = 0.00000 0M 3.00:000 = 59:.0000:0: BMS = 1.55200 BÑH = 14.58000 = .6800.0 ULM = .02130 DS = 10.0000,0 N 20 = 23.6000.0 ON = 1.99040 CM = .06250 CD = 1.90.000 AMP 1.00000 .4750.0 ALAMa 0.00000 ALAM2 = .41500 BLAM1' .05000. BLAM2 = 0.00000 o.00000 AGAMi .05000 AGAM2 = 0.00000 BGAM1 0.0000.0 BGAM2 s 0.00000

STEADY-STATE SOLUTION OYNAHIC SOLUTION PHM(J) DPK(J) TM(J) DTIJ) UM'(J) DU(J) V14(J) DV(J) J .008842 .831901 301.489251 -2.814459 .331768 -.704693 1.489530 -1.362528 1 .005794 1.576827 302.732956 -2.807222 015479 1.880557 1.519121 1o264096 2 .005234 2.339834 303.882394 -2.800017 .237110 2.793122 1.494022 -1.209241 3 .005682 2.839033 305.049831 -2.792912 .372358 3.009067 1.450318 _1.:189389 . 4 .006049 3.106924 306.306609 -2.785766 .448404 -3.078312 1.406860 e.1-.193621 5 .006048 -3.060084 301.687693 -2.778396 .494144 -2.899734 1.364827 -1.209794 6 .005691 -3.061654 309.209469 -2.770663 .531310 -2.750316 1.312969 -1.227136 7 .005088 3o107297 310.884818 -2.762504 .571834 2.646450 1.235025 -1.240898 8 .004406 2.865046 312.730100 -2.753937 .619215 -2.599972 1.116202 -1.265479 9 .003893 2.470387 314.764821 -2.745055 .672496 -2.614224-.947426. -1.237930 1:0 .003873 1.946461 317.006731 -2.736011. .729948 -2.686214 ..27655 _L.210339 11 .004559 1.416073 319o464968 -2.727006 .791120 -2.810368 .465438 -1il33550 12 .005893 .972827 322.133194 -2.718269 .857743 -2.980820 .187316 -.792759 13 .007739 .611171 324.983988 -2.710038 .934251 3.091259 .172282 1.203138 14 .010034 .300525 327.965323 -2.702531 1.028245 2.845614 432074 1.589243 15 .012783 .019702 330.999525 -2.695925 1.150712 2.573098 .651397 1.683628 16 .016029 -.243020 333.984878 -2.690330 1.315609 2.283401 .782201 1.735628 17 .019836 -.493996 336.799172 -2.685773 1.538719 1.9.87242 .788414 1.783612 18 .024286 -.736794 339.309145 -2.68217.8 1.836279 1.694194 .639274 1.858182 19 .029477 -.973510 341.372911 -2.679355 2.224183 1.410826 .31980,7 2.116448 20 P110(J) T0(J) OU (J) VO(J) CXO(J) CYO(JI J 1.429525 6100.777082 23.364892 3.3292l 0.000000 0.000000 1 1.325929 6113.771030 22.896003 5.72*281 1.921487 9.809099. 2 1.227636 6129.825151 22.224023 7.940580 4.823133 19.374659 3 1.135983 6148.608529 21.403985 9.941299 8.619797 28.622110 4 1.051700 6169.768338 2ú.491126 11.70(850 13.214786 37.500561 5 .974994 6192.959022 19.533680 13.243691 Í8.509051 45.981229 6 .905677 62L7.861785 18.569499 14.64811 24.407890 54.053629 7 .843314 6244.194678 17.625688 15.693793 30.824960 61.721014 8 .787333 6271.715587 16.719984 16.65b394 37.684075 68.996069 9 .737111 6300.220740 15.862777 17.473761 44.919480 75.897365 10 .692026 6329.540867 15.059120 18.170936 52.415242 82.446679 11 651492 6359.536466 14.310620 18.766263 60.304203 88.667102 12 .614975 6390.093060 13.615734 19.276198 68.366795 94.581187 13 .581994 6421.116843 12.972702 19.714690 76.629882 100.213157 14 552127 6452.530925 12.378189 20.093293 85.065711 105.582455 15 .525004 6484.272167 11.828711 20.421596 93.651000 110.709508 16 .500302 6516.288590 11.320703 20.701527 102.366170 115.612642 17 .477742 6548.537265 10.850686 20.951638 111.194699 120.308688 18 .457079 6580.982604 10.415361 21.177352 120.122592 124a813037. 19 .638103 6613.594998 10.011657 21o371166 129.137953 129.139729 -20

CXM(J) DCX(.J) CYM(J) OCY(J) J 0,000000 o. ô00000 0.000000 0.000000 1 .087537 -2.303692 .012449 .837901 .056209

-j

.564765 .014046 1.576827 3 049291 -.801759 .017611. 2.3398i4 4 .051531 ..3Q2560 .023034 2.839033 5. .052518 -.034669 oo0O-i 3.10694 6 .050056 .08 1509 .0339:37 -3.0600b4 7. .044777 '07993 .035120 -3.0616b4 8 037998 -.034296 .033833 3.107297 9 .031219 -.276547 .031098 2.865046 10 .026169 -.671205 .Ó28826 2.470387 11

047

13 -1.195132 .029820 1.946461 12 .027646 t .25520 .036254 1.416073 13 033999 ?. 168766 .098134 97287 14 042539 53042? 064647 .611171 15 052629 -2.841067 .085431 .30055 16 .064073 -.121851 .110618 .019702 17 .076891 2 898573 .140647 -.243020 18 .091202 2.647597 .176153 -.493996 19 .107 182 2.404798; .217931 -.736794 20

5

10

15

20

APPENDIX C

LISTING AND SAMPLE OUTPUT OF OMFLO

PROGRAM OMFLO (INPUT.OUTPUT,TAPE51NPUT,TAPE6=OUTPUT)

DIMENSION PHO(2OO) TO(200), UO(402), V0(200). PHRA(200),

PHRB(200), PHIA(200)' PHIB(200), TRA(200),TRB(200), TIA(200),

TIB(200), URA(200), URB(2OO),.UIA(2OO) UIB(200), VRA(2O0)

. .VRB(200), VIA(OO), VIB(200),PMR(200), PHI(200). TR(200),

TI(200), IJR(200), UI(200), VR(200), VI(200), FL(200), FG(200),

DFL(2OO) DFG(2OO) 1(8), QPHO(200), QTO(200), QPMRA(200),

+ QPHRB(200)v QPHIA(200), QPHIB(OO), QTRA(200), QTRB(200),

QTIA(200), QTIB(2OO) QURA(200), QuRB200, OUIA(200),

QUIB(2O0) QVRA(2OO) QVRB(200), QV1A(200), QVIB(200)

COMMON DS,WUL,D,ALAMO,ALAM1.ALAM2,BLAM1, 1BLAI42,AGAMO,AGAM1,AGAM2,BGAM1 BGAM2t 2OMiULM,APP,C,PHO,DFL,DFG,K,.I ,TOPHR,PIf I FL,FG EXTERNAL DAUXDARN READ (5,11) TITLE, BLO, BDO, BW, BLM DBL, 6DM, DBD, 0M, . BM, BMS, BMH,

+ WUL ULM OS, N,

C, DN CH, CD, AMP,

ALAMO, ALAM1, ALAM2, 8LAMI., BLAM2,

AGAMO, AGAM1 AGAM2, BGAM1, BGAM2

C WRITE TITLE ANÒ INPUT DATA

WRITE (6,21) TItLE

25 WRITE (6,22) _BLD, BOO, 8W

WRITE (6,23) BLM, DBL, BON, DBD, 0M WRITE (6,24) BM. BNS, 8MM

WRITE (6,25) WUL, ULM, DS, N

WRITE (6,26) C ON, CH, CD, AMP

30 WRITE (6,27) ALAMO, ALAM1, ALAM2, BLAM1, BLAMe

WRITE (6!28) AGAMO, AGAM1, AGAM2, BGAM1, BGAMC

C INITIALIZE D O.5DNCCCH*CD -AMC = 3.1.4159*CHÖCH*DN*AMP/4.O 35 APP = 1.0 + AHÇ/ULM CXO 0.0 ÇYO = 0.0 .-CXR. = 000 -CXI-O.O -40 CYR = 0.0 CYI 0.0 KA = N-1 IA = Nel

C SOLVE STEADY-STATE PROBLEM

45 ERR.0000I ERA=.001 CALL ATA(BDO,YK) P1108 = ATAN2(BLO-BW, -YK) TOB = SQRT(8008P0 (BLO-BW)(BLO-BW)) 50 WRITE (6,31) WRITE (6,32) P110(1) = P1108 . TO(i) = TOB

C MULTIPLY TO FINO VELOCI1Y COMPONENTS

60 65 70 75 80 85 90 95 100 los 1 10 U0(K) = CSIN(PH0(K)) V0(K) C*COS(PHO(K))

J =K

WRITE (6,33) PHO(K), 10(K), U0(K) V0(K), CXO, CY0 J

1F (K.EQ.N) GO TO 101

Ç INTEGRATE TO FIND EÑSiÔN. ANGLE, AND CABLE CONFIGURATION

T(1)=T0(K) -T(2)PHOtK) 1(3) = CÁO 1(4) = CYO FIRST.O ALS=FLOAT(K)!ÒS CALL KUTMEIU4,ALS,T,ERR,DS,FIRSI9HCX,ERA,DAUX) T0(Kl)=T(1) P140 (K+1)=T (2) CXO T(3) CYO = 1(4) 101 CONTINUE

C SOLVE DYNAMIC PROBLEM

RtTE (6,34)

WRITE (6,35)

C REVERSE LISTS F PHO AND TO

bo 301 I = ioN QPHO(I) = PHO(I) 010(I) = TO(I) 301 CONTINUE DO 302 I 1,N -L = N+1I P110(I) QPHO(L) TO(I) = Q 0(L) 302 CONTINUE

C COMPUTE VISCOUSFORCE TERMS ALONG CABLE

DO 201 I = 1 IA -FL(I) = D*(,ALAMO_ALAM1*COS(PHO(I))ALAM2*C0S(2.0*PHO(I) + BLAM1*SIN(PH0(I))_BLAM2*SIN(2.U*PHO(I))) FG(I) = D*(_AGAMO+AGAM1*COS(P1f0(I))_AGAM2*COS(2.0*PHO(I)) _BGAM1*SIN(PHO(I))+BGAM2*SIN(2.O*PH0(I))) DFL(I) = O*(+ALAM1*SIN(PHO(I))_2s0*ALAM2*SIN(2.0*PHO(I)) BLAM1*COS(PHO(I))_2.O*BLAM2*ÇOS(2.0*PH0(I))) DFG(I) = D(_AGAM1*SIN(PH0(I))+s0*AGAM2*SIN(2s0*PH0(I)) + BGAM1*cos(PH0u) ) ,2.O*BGAM2*CUS(2.0*PHO(fl)) 201 CONTIÑUE

C COMPUTE MODAL SOLUTIONS

C COMPUTE A MODE SOLUTION

C SET VELOCITY COMPONENTS EQUAL TO ZERO AT BOTTOM

URA(1) 0.0

UIA(1) = 0.0 VRA(1) = 0.0 VIA(I) = 0.0

C ASSIGN VALUES TO ANGLE AND TENSION AT BOTTOM

PHRA(1) = 0.001

PHIA(1) = 0.0

TRAU) 0.0

TIA(i) 0.0

115 120 125 130 135 140 145 150 155 160 165 C THE CABLE DO 202 1= 1' IA T (1) =PHRA (I) (2) PHIA (I) T (3)=TRA C I) T(4)TIA(1) T(5)=URA(I) T(6)=UIA(I) T(7)=VRA(I) T(8)VIA(I) F I RST O ALS=FLOAT(I)DS CALL KUTMER(8,ALS,T,ERR,DS,FIRST,HCX,ERA,DARN) PHRA(I+1)T (1.) PHIA(:I.1) =1 (2) TRA(I1)=T(3) TIA(I.1)T(4) URA(i+1)T(5) UIAU+i)=T (6) VRA(i+1)T(7) VIA(.I+1)=T(8) 202 CONTINUE

C COMPUTE B -MODE SOLUTION

C SET VELOCITY COMPONENTS EQUAL TO LERO AT BOTTOM

URB(1) = 0.0 UIB(1) = 0.0

VRB(1) = 0.0

VIB(1) = 0.0

C ASSIGN VALUES TO ANGLE AND TENSION AT BOTTOM

PHRBU) = 0.0

PIIIB(1) 0.0

TRB(1) = 0.1

TIB(1) = 0.0

C FIND VELOCITY COMPONENTS, ANGLE, AND TENSION BY INTEGRATING UP

C THE CABLE DO 203 1=1, IA T(1)PHRB(I) T(2)PHIB(I) T(3)TRBU) T(4)=TIBU) T(5)URB(1) T(6)UIB(I) T(7)=VRB(Í) T(8)=VIB(I) FIRST.0 ALSFLOAT(I)'DS CALL. KUTMER(8,ALS,T,ERR,DS,FIRST,HCX,ERA,DARN) PHRB(I.1)ZT(1) PHIB(I.1)ZT(2). TRB(I+1)T(3) TIB(I1)T(4) URBU.Ì)f(5) uIBiI.Ï=T6) VRB(I.1)=T(7) VIB(I+1)zT(8)

i 70 175 180 185 190 195 200 205 210 215 220 203 COÑTINUE

C REVERSE LISTS 0F PHO , TO , AND MODAL SOLUTIONS

DO 303 K 1,N QPHO(K) = PH0K ÓTO(K) = 10(K) QURA(K) = URA(K) QUIA(K) UIA(K) QVRA(K) = VRA(K) QVIA(K) VIA(K) QPHRA(K) PHRA(K) QPHIA(K) PHIA(K) ÖTRÀ(K) tRA(K) QTIAIic TIA(K) QURB(K) = URB(K) QUIB(K) UIB(K) QVRB(K) = VRB(K) QVIB(K) = VIB(k) QPHRB(K) PHRB(K) OPHIB(K) = PHIB(K) QTRB(K) TRB(K) QTIB(K.) TIB(K) 303 COÑTIÑUE DO 304 K = 1, N

= N1K

PHO(K) QPHO(M) 10(K) = QTO(M) URA(K) = QURA(M) UIA(K) = QUIA(M) vRAUc = QVRA(M) VIA(K) = QVIA(M) PHRA(K) QPHIA(M) PHIA(K) QPHIA(M) TRA(K) = QTRA(M) TIA(K) = OTIA(Ñ) URB(K) = QURB(N) UIB(K) = QU1B(M) VRB(K) = QVRB(M) VIB(K) = QVIB(M) PHRB(K) = QPHRB(M) PHIB(K) = QPHIB(M) TRB(K) QTRB(M) TIB(K) = QTIB(M) 304 CONTINUE

C SUPERIMPOSE MODAL SOLUTIÒÑS

c COMPUTE FORCE COMPONENTS AT TOP

BDR =. _BDM*COS(BDB) BD.! BDMSIN(BDB) BLR BLM.COS (DeL) BL! = BLMSIN(bBL) PR = - BDRCOS(PH0(1)) + BLR*SIN(PHO(1)) PI = - BDICOS(PH0U)) BLI*SIN(PHO(1)) OR = BDRSIN(PH0(1)) BLR'COS(PHO(l)) Q! = BDI'SIN(PHO(l)) BLI*COS(PHÖ(i))

PAR TRA(i) - OM*(.(BMHBMS)*UIA(1)COS(PHO(1))SIN(PHO(1))

25 230 235 240 245 250 255 260 265 270 275 . VI*(1))

PAl n TIM1) - OM*(-(BMH-BMS)URA(1)COS(PHOt1))'SIN(PM0(1))

.. +(8M.BMSCOS(PH0(1) )COS(PH0(1) ).BMH'SINIPHO(l) )eSIN(PHO(1)))*

. VRA(1)) .

PBR z TRB(1) -. OM*(.ÇBHHBHS)*UIB(1)*ÇOPH0(1))'SIN(PHO(1))

. -(sM.BMs.coscpHo(1))'cos(pHo(1)).BMHesIN(pHo(1))'sIÑ(pHo(1))'

. VIB(1)) . .

PBI TIT1) - oM*((BMH_BPss)*uRs(1)*cos(pI.so(1.))eSIN(pHo(1))

.BNBMsecosoc1.cosPH01,BMH*S1N(PH0(:1J)*SIN(PH0(1)))* . VRB(1)) OAR = T0(1)*PilIA(i)OÑe(+(BM.W4S*SIN(PH0(1) )S1N(PH0(1)) e BMH*COS(PM0(1))ftCOS(PH0(1)))*UIA(1(BMH.BNS)*VIA(1)COS(PH0(1)) SIN(PHO(1))) QAI = TO(1)'*PHIA(1)-OM*(-(BM.MS'SIN(PH0(1))SIN(PH0(1)). BMH.COS(PH0(1))'COS(PHo(1)))'I»A(1r.(BMH-8MS)VRA(.i)C0S(PH0(1)) SINCPHO(1))) QBR = _{1o(1).PHRB(l)-oMec.(BM.BMs!sIN(pHo(l))sxN(pHo(1)).} BMH*COS(PHO(1))*COS(PHO(l)))*UIB(1)_(BMH_BMS)*VIB(1)COS(PHO(l)) fSIN(PHO(1))) OBI = TO(.l)*PHIB(1)0M*((BM+tSMSSIN(PHO(1))SIN(PHO(1))' BMH*COS(PHO(1))*COS(PHO(l)))*UMB(l)(BMHBMS)*VRB(1)'COS(PHO(l:)) SIN(PHO(l))) DENOMR = PAR*QBR-PAI'QBI-P8R*UAR.PBIOAI

DENOMI = PA! QBR.PARQBI-PBIQAR-PBRQAI

AiiIJMR = P RQBR-P I'QBI-PBRU R.PBI*Q I

AP4UMI = P 1.QBRpaQBI-PBIu R_PBRÓQ i

BPIUMR = PARQ RPAI*Q I-P R*UAR+P I*QAI

BÑIJMI = PAIQ R.PAR*Q I-P I*UAR-P RÒQAI

DENOM2 = DEN0MRDENOMR DENOMIDENOMI

CAR = (ANUMR'DENOMR. ANUMIDN0MI ) /DENOM2

CAl = (-ANUP4RDENOMI + ANUMIDNOMR)/ÔENOM2

-CBR (BNUMRDENOMR BNUMIÔÈNOMI)/DENOM2

cBI = (-BNUMRDENOÑI + BNUI4I*DÈNOMR)/DENOM2

C COMPUTE VELOCITY COMPONENTS ANGLES AND TENSION ON CABLE

DO 204 K 1 N

UR(K) = CAR*URA(K) - CAIUIA(K) CBR'URB(K) - CBIUIB(K)

UI(K) = CAI.'IJRA(K) CAR'UIA(Ic) + CBIURB(K) CBRUIB(K)

VR(K) = CARVRA(K) - CAIVIAÙc) CBRVRB(K) - CBI*VIB(K)

VI(K CAIVRA(K) CARVIA(K) +CBI*VRB(K) CBI*VIB(K)

PHRK CAR*PHIAtIO - CAI*PHIA(K) CBR*PHRB(K) - CBI*PHIB(K)

PHI(K) .CAI*PHRA(K) CAR*PHIA(K) CBI*PHRB(K) CBRPHIB(K)

TR(K) = CAR*TRA(K) - CAITIACIc) CBRTRB(K) -CBI'TIB(K)

TI(K) =. CAI*TRA(K) CARTIA(K) CBI*TRBK CBR*TIB(K)

C COMPUTE PHYSICAL VELOCITY COMPONENTS ÒN CABLE

URP = _{C*PHR(K).ØCOS(PHO(K))} + UR(K) «

tJIP = CPHI(K)ÇOS(PHO(K)) + UI(K)

VRP = - CPHR(K)SIN(PHÖ(K)) + VR(K)

VIP = - C*PHI(K)*SIN(PHO(K)) VI(K)

COMPUTE MAGNITUDE AND PHASE 0F VELOCITY COMPONENTS ANGLE, ANO

TENSION ON CABLE

UM SQÑT(URP*URP UIP*UIP)

CALL ATA(URP, YK)

DU = ATAÑ2(UIP YK)

VM = _SQRT(VRP*VRP VIP'VIP)

DV = ATAN2(VIP. Vic)

PHM SURT(PHR(tc)*PHR(K)+ PHI(K)'PHI(K))

CALL ATA(Phk.(PÇ), Yic)

DPH A1AN2(PHI(K), YK)

280 TM = SQR1(TR(Ic)TR(K) + TI(iTI(K))

CALL ATA(1R(). YK)

DT = ATAN2(TI(K). Y(()

J = -K

WRITE (6,36) PHM, DPH, TM. DT. UM, DU VM. DV J

285 204. CONTINUE

WRITE (6,37)

C COMPUTE MAGÑI LUDE AND PHASE OF DYNAMIC CA8LE CONFIGURATION

UO 205 Pc = l.N

CXM SQRT(CXR*CXH+CXI*CXI)

290 CALL ATA(CXR VI<)

DCX = A1AN2(CXI VI<)

CYM SQRT(CYNCYR+CYICYI)

CALL ATA(CYI, YK)

DCV ATAN2(CYI. VI<) 295 J = - K WRITE (6,38) CXM, DCX. CYM. DCV, J IF (K.EQ. N) GO TO 205 CXR. = PHR(K)*SIN.(PH0(K))*DS CXI PHI(K)*SIN(PHO(K))*DS 300 CYR = _PHRiPc)*CÔS(PHO(K))*DS CYI _PHI(K)*COS(PHO(K))*DS 205 CONTINUE 11 FORMAT ( 10*, A6 I lOX, 3F10.5/ 10*, 5F10.5/ 10*, 3F10.5/

lOX. 3F10.5. 110/ 10*, SF10.5 / lOXó SF10.51 ioX. SF10.5)

305 21 FORMAT (1HL. lox. 53HDEEPLY SUBMERGED SUBMARINE TOWING A FLOAT IN

GA SEAWAY .1 lH- 13X. A6///)

22 FORMAT lÀ, 8HBLO , FiO.5, 4X 8HBDO = , F1O.5. 4*, G 8HBW , FlOaS II)

23 FORMAF ( 1*, 8HBLM = F10,5, 4*. 8HDBL = , F10,5, 4*.

310 HBDM , F10,5, 4*. RHDBD , F10a5, 4*, 8140M =

F1O.5 /1)

24 FORMAT ( lX, SHBM , FlOoS. 4x, 8HBMS = , F1O.5. 4X,

G BHBMH FlO.5 II)

25 FORMAT ( 1*. RHWUL = . F10,5, 4*. 8HULM = , F10.5, 4*,

315 BHDS = . F10,5.. 4X, 814N = , 110 II)

26 FORMAT ( ix, 8HC = , Fi0.5. 4*. 8HDN = . F10,5, 4*,

+ BHCH = , F10.5ö 4*. 8HCD = . F10,5, 4*. 8HAMP = ,F1O.5//) 27 FORMAT ( lx. 8HALAMO = ,F10.5. 4*, 8HALAM1 F10.5.4X.

BHALAM2 = ,F10.5, 4*. 8HBLAM1 = . Fi0.5, 4*. 8HBLAM2. =

320 FiO.5 II)

28 FORMAT ( lX RHAGAMO = . F10,5, 4*, 8HAGAM1 = , F10,5. 4*. BHAGAM2 .= , F10,5, 4*, 8HBGAM1 = , F10.5, 4*, 8HBGAM2 =

FlOeS /1)

31 FORMAT C 22H1STEADY-STATE SOLUTION)

325 32 FORMAT. C 1HO. iX. 6HPHO(J). 10*. 6HTO(J) lOX, 6HUO(J) lOX,, + 6HVOCJ) lox! 6HCXO(J). i0X 6HCYOCJ). 10*. 1HJ

II)

33 FORMAT C 2*. F9,6. 7*, Fl3.6. 3)( 21 F11.6. 5*). 2(F13.6. 3*). 14)

34 FORMAT ( 17H1DYNAMIC SOLUTION)

35 FORMAT C 1KO, ix, 6HPHM(J). 9*. 6HbPH(J). 9*, 6HTM(J) 9*.

9X. 1HJ I/i

36 FORMAT ( 2A. 2( F9.6, 6*). F12 6, 3*, F9.6. 6*, 2(F10 6. 5*,

F9.6. 6*). 14)

37 FORMAT ( 1Ml. LX, 6HCXM(J). lOX. 6HDCX(J), lox, 6HCYM(J). 10*.

335 + 6HDCY(J). 10*. 1HJ /1)

38 FORMAI ( ¿X, 2( F12.6, 4*. F9.6, 7X). 14)

STOP

ED

-SUBROUTINE ATA(X.YK)

THIS SUBRÔUTINE PREVENTS ATAN2 IROM DÏVIDING BY ZERO WHEN

COMPUTING THE ARCTANGENTS OF P112 AND -3PI/2

= 0.00001 IF (ABS(X).GT.ER) -GO TO 10

Y=ER

RETURN lO YP(=X RETURN 10 END SUBROUTINE KUTMER(N,T,Y0.EPS,H,FIRST,NCX,A.DAUX)

C THIS SUBROUTINE PERFORMS FOURTH-ORDER KUTTA-MERSON INTEGRATIONS

DIMENSION YO (10) .Y1(10) .Y2(l0) ,F0(I0) F1 (10) .F2(1O)

IF(FIRST)20,10,20 10 HCH KUTMOO5O IFLoc1 KUTNOO6O FIRST1. UTM007O 20 L000 KUTNOO8O HCXHC KUTMO 090

30 CALL OAUX(T,YO.F0) _KUTMO100

39 DO 40 Ii.N _KUTMO11O

40 Yl (I)Y0(I).(HC/3.)FO(I) KUTMO12O

CALL DAU*(I.HC/3..Y1,Fl) KUTMOL3O

DO 50 11.N KUTMO14O

50 Yl (I)=Y0(I),(HC/6.)*F0(I.(HC/6.)eF1 (I) KUTHO15O

CALL DAUX(T.HC/3..Yl.Fl) KUTMO16O

DO 60 11,N KUTMO17O

60 Yl(I)=YO(I).HC/8.'F0(I)..375HCFj(I) KUTMO18O

CALL DAUX(T.HC/2..Yl.F2) _UTM0190

DO 70 11,N _KUTMO200

70 Yl(I)YO(I).HC/2.F0(I)-L.S.HC*F1(I).2..HC.F2(I) KUTMO2IO

CALL DAUX(T.HC.Yl,Fl) _KUTMO22O

DO 80 Il,N KUTMO23O 80 Y2(I)Y0(I).HC/6.!F0(I)..66666667.HC.F2(I).(Hc/6.)1F1(I) KUTMO24O INC=0 KÙTMO25O 00 110 11.N KUTh0260 ZZZ=ABS(Yl(I))-A _KUTÑO27O IF(ZZZ) 85,87.87 _KUTMO28O 85 ERROR=A8S(.2(Y1(I)-Y2(I))) KUTMO29O IF(ERROH-A) 100,100,90 KUTMO300 87 ERfiOR=ABS(.2-.2*y2(I)/y1(I)) _KUTMO31O

IF(ERROR-EPS) ioo,ioo.00 _KUTÑ032O

90 X=128.'ABS(HC)-ÀBS(H) _IWTMO33O

1rx 9is.os

KUTMO34O

91 WRITE(6.92 T,ERROR KUTMO3SO

92 FORMAT (Z1H RELATIVE ERROR AT T= 1P1E12.3,3HIS Fl0.6 ) KUTMO36O

FIRST 2a _KUTMO37O

RETURN _KUTMO38O

95 HCHC/2. _KUTMO39O

IPLOC=2 *IPLOC _KUTMO400

LOC2 LOC KUTMO41O

HÇXHC - KUTMO42O GO TO 30 _KUTM0430 100 IF(ERR0R64.-EPS)110.110i1Ol _KUTMO44O 101 INC1 _KUTMO4SO 110 CONTINUE _KUTMO46O 111 T=T.HC _KUTHO47O DO 112 11.N KUT140480 112 Y0(I)Y2(I) _KUTMO49O LOC=LOC+1 KUTHO500 ¡F(LOC-IPLOC) 120,210.210 _KUIM0510 120 IFUÑC)210.l-30,210. _KUTMOS2O 130 IF(L0C-(LOC/2)2)2l0.14O,2l0 _KUTMO53O 140 IF(IPLOC-1)210.210,200 _KUTMOS4O 200 HC2.HC KUTMOSSO 5 10 15 20 25 30 35 40 45 50 55

5 10 15 20 25 5 10 15 20 25 30 LOCLOC/2 I PLOC: I PLOC/2 210 IF(IPLOC-LOC)30,220,30 220 RETURN 60 END KUTMOS6O KÚTMO5TO KUTMO58O KUTMO59O KUTMO600 SUBROUTINE DAUX(Z,TX,F)

THIS SUBROUTINE PHOVIDES THE INTE(iRANDS FOR KUTMER WHEN IT

COMPUTES THE S1EAOY-STATE TENSION, ANGLE, AND CABLE CONFIGURATION

IMENSION _{tX(8), F(8)PH0(20O)'DFL(2O0),DFG(2O0)T0(200),} PHR(200), PHI (200) COMMON DSWÚL,D,ALAM0,ALAM1,ALAM2,BLAM1, 1BLAM2,AGAMU,AGAMÌ,AGAM2,BGAM1,BUAM2, 2OMULM.APP.C,PH0,DFL,DFG,K, I, T0,PHR,PHI TX 1TX (1) 1X2=TX(2) C1=COS(TX2) S1=SIN(TX?) C2COS (2.*1X2) S2=SIN(2.*1X2) FL=D* (ALAMO_ALAM1*C1+ALAM2*C2+RLAM1*S1_BLAM2*S2) FG=D* (_AGAMO+AGAM1*C1áAGAM2*C2_bGAM1*S1 .BGAM2S2) IF(TA1.EQ.Ú.) GO TO SO F (2) (WUL*C1FL)/TX1 = -Cl - Si 7 f(1)(WUL*S1,FG) RETURN 50 F(2)0.0 GO TÒ 7 END SUBROUTINE DARN(Z,T,F) V

THIS SUBROUTINE.PROVIDES THE INTEGRANDS FOR KUTMER WHEN IT

COMPUTES THE DYNAMIC ANGLE, TENSIOÑ ANO VELOCITY COMPONENTS

DIMENSION T(8), F(8),PHOi200).,OFL(200),DFG(2Ö0),TO(200), FL(200),FG(20O),PHR(200),PHI(200) COMMON DS,WULcD,ALAMOcALAM1,ALAM2BLAM1 1BLAM2.,AGAMOçAGAM1,AGAM2,BGAM1 ,BGAM2, 20M,ULM,APP,C,PHO,DFL,DFG,K,ITO,PHR,PI$I,FL.,FG C1COS(PHO(I)) s1=sIN(Po1:w IF(TO(I).EQ.0.0) GO TO 50

r(1)=(-OM.ULM*(APPeT(6) .(APF-i.0)*C'T(2)*C1)-T(3)*(PHO (1.1)-PMO (I) )1DS-(T(1)(T(5)C1-T(7)S1)/C)DFL(I)

_2.0FL(I)e(T(5)*S1+T(7)*C1)/CsWUL*TU)*S1i/T0 (I) F(2)=( OM*ULM*(APP*T(5),(APP_1.U)QCT(1)*Ci)uT(4)*(PH0(I,i)PH0(I)i

)/DS_ÇT(2).(T(6)*Ç1_T(8)*Sl)/Ç)*DFI..(I)

+ V 2.O*FL(I)ó(T(b)*S1+T(8)*C1)/CWUL*T(2)*S1)/T0(I),

7 F(3)=C OM*ULM*(T(8) (APP1.0)*C*T()*S1) +(T(i)

.(t(5)C1-T(7)*S1)ÌC)*OFG(I)+.O*FG(I)* (T(5)S1T(7)C1)/CWUL'T(1)C1) F(4)=(-OHÚLM(T (7)(APP1.0)*C*T(1)*S1)+ (1(2) +(T(6)'Cl-T(8)*Si)/C)*DFG(I).2.O*FG(Ì)*. + _{(T(6).S1+t(8)*C1)/C+WUL*T(2)*CI)} F15)_QM*T(2)+T(7)*(PHO(I+1)_PHO(I))/DS V F(6)= OM*T(Ï)+T(8)(PHO(I+1)PH0(I) )/DS V F(1)-T(5)(PHO(I+1)-PH0(I) )/DS F(8)=_T(6)(PH0(I.1)PH0(I))/DS RETURN V V 50

F(i)0.0

V F(2)=0.0 GO TO 7