The application of padé approximants in the time domain simulation of the non-linear motion of a semi-submersible excited by waves

W.G. Price * E.H. Twizell + Yigong Wang *

Summary

This paper discusses one of a family of numerical methods suitable for the time domain simulation of systems described by non-linear second order

differential equations. The basic approach involves Pads approximants [1] which are capable of solving problems exhibiting oscillatory and growing solutions, that is,instabilities. The particular Pads approximant described is the (1,1) formulation which, combined with predictor-corrector procedures, produces a mathematically stable technique allowing accurate and convergent solutions to be determined. The technique permits the evaluation and description of the

behaviour of deterministic and random time varying non-linear systems and the method is applied to simulate the responses of a semi-submersible excited by waves.

1. Introduction

THE APPLICATION OF PAD g APPROXIMANTS IN THE TIME DOMAIN SIMULATION OF THE NON-LINEAR MOTION OF A

SEMI-SUBMERSIBLE EXCITED BY WAVES

* Department of Mechanical Engineering, + Department of Mathematics & Statistics,

Brunel University of West London Uxbridge, Middlesex, UB8 3PH.

-TECHNISCHE UNIVERSITET

Laboratorium voor

Scheepshydromechanica

Archief

Makelweg Z

2628 CD

Delft

VOL- 016

men

-Fax.015- 781838

U.K.

Although much physical insight into

fluid-structure

problems may be deduced from a linearised mathematical model involving small amplitude motions, when interest turns to large motions the mathematical model inevitably becomes non-linear in one form or another. The actual model chosen is always open to question but once adopted a suitable method of solution must be sought which reflects the correct physical characteristic of the problem.

Many different methods are available to solve non-linear equations [2-6]

(i.e. perturbation expansions, equivalent linearisation, Runge-Kutta numerical procedures, etc), each having advantages and disadvantages. For example, the equivalent linearisation approach is relatively easy to implement but it is surprising that after a detailed discussion on the development of the non-linear mathematical model describing the system, the resultant set of non-linear equations

kr

is reduced to a set of linearised equations for solution. At times this approach may be very satisfactory but the technique introduces uncertainties since some legitimate solutions may have been discarded in its implementation. Again, the

fourth order Runge-Kutta numerical method is widely used by engineers to create time domain solutions describing the behaviour of non-linear systems. If the time step increment size is incorrectly chosen [5], a mathematical instability may be introduced into the solution which is not necessarily a physical

characteristic of the system. However when analysing a non-linear system which displays instabilities through growing oscillatory solutions, doubts about the numerical procedures adopted can cause extra work in verifying that the actual unstable solution is real and not an idiosyncrasy of the mathematics. Further,

a restriction on the maximum incremental time interval produces an increase in the number of computations and also the cost of solution.

To overcome these restrictions and difficulties a numerical method was sought which satisfies the following conditions

computationally efficient, producing accurate and convergent solutions,

the method is mathematically stable and is able to generate growing oscillatory solutions exhibiting instability in the linear or non-linear system,

no severe restriction on the time step interval,

independent of the time description of the functions involved, i.e. sinusoidal, random, etc.

To this end a numerical method was devised involving the (1,1) Pad4 approximant combined with a predictor-corrector procedure. This approach

satisfies all the

above requirements and the technique is briefly illustrated

with reference to the evaluation of the non-linear motions of a semi-submersible excited by waves.

Numerical Technique

In ship and offshore dynamic problems, non-linear descriptions of the fluid actions may be introduced into the mathematical model in many different ways. For example, in the expressions describing damping forces, restoring forces, wave actions, etc. As an example to illustrate and demonstrate the (1,1) Pad4

approximant technique, developed to solve typical time dependent, non-linear equations encountered in marine dynamics problems, let us consider a second order non-linear initial value problem described by the generalised equation

A P(t) + B fp(t) + C p(t) + D IP(t) P(t) = c(t) (1) P(0) - Po ,

P(o) - Po

399 (i) =

where A, B, C, D are N _{x N matrices independent of the column}

N-vector respons ,

p = p(t) = _{[pit, p2(t),...,p (t)jT}

, with superscript T denoting

transpose. The column N-vector E(t) _{represents a generalised external}

applied force due to

-the seaway which may be sinusoidal or random in form. _{In more specific problem} equation (1) may be modified such that the N _{x N matrices are time dependent,} the function ip(t)1 _{L(t) is replaced by a term y(t)}

ly(t)1, where y(t) = p(t)

-.

.

-

-

_-

-C(t) is the relative velocity

between structure and wave (C(t)), the inclusion

-of higher order polynomials

(i.e. pn(t),p (t)) and combinations of such.n _terms, etc.

A family of numerical methods based on Pad-6 _{approximants [1] was developed}

to solve a system of

general nonlinear second order ordinary differential equations of order N.

This system was transformed into a system of first _order equations of order 2N and numerical methods derived _{by replacing the exponential} term in a recurrence relation

by a number of approximants. _{It was found that} the (1,1) Pade approximant

version satisfied the requirements stated previously, i.e. conditions (i-iv).

The substitutions u = u(t) = p(t), _{v = v(t) = i5(t) transforms equation (1)}

- _{_}

-into a system of 2N first order ordinary differential

equations represented by the matrix form

q(t) = a q(t)

+ a (q(t)}

_{-,- fB(t)} (2) T

Tr

c(o) =

[

00 , po . q =1.

-0

_ where q q(t)

[11T, v1T

[0T,

Ivl

(A-1 D

v)Tr

3(t) [0T, (A-1 , (t))

1T

are vectors of order 2N and

I a = 1 -A-1 C A -B is a matrix of order _{(2N x 2N).}

Furthermore, by defining the (2N x 2N) diagonal differential matrix operator L in the form

E I d/dt = _{d/dt 0} 0 0 _d/dt 0 . . d/dt 400 L

equation (2) may be written as

Lq=aq+afg) +

(t)

q(0) = q _o

If T denotes a small incremental time step, then the Taylor series _expansion of q(t + T) about q(t) is q(t + T) = q(t) + T -d- q(t) +

<12-dt2

q(t) + dt 2! 2 = + T L + 1 T L2 ...1 q(t) = exp(T L) q(t) (3) (4)

and a family of numerical methods can now be derived based on an approximation made to the exponential operator term exp(T L). In particular it can be shown that the exponential operator may be expressed as [7,8]

(I - E T L) exp(T L) = (I + (1-E) T L)

-1 ,

Or exp(T L) = (I - E

T L)

L I + (1-E)

T L ]

where E is some parameter lying in the range 0 5 5 1. The values E = 0,

1,-2 are equivalent, respectively, to the (0,1), (1,0), (1,1) PadS approximants to exp(T L) and it is the last of these which is now discussed.

2.1. The (1,1) Pads aPproximant

Using equation (5), with E = 0.5, in equation (4) gives

T L) Q 1 T L)

(I n+1 = (I +

2 2

where upper case Q has been introduced to distinguish the theoretical solution of a numerical method from the theoretical solution of equation (2) _{and Qn = Q(n T)} for n = 0,1,... _{Using equations (3) and (6) gives the implicit method}

1 , n 1 n 1 ,n

(7)

n+1 1 n+1 1 n+1

(I - 1 T a) Q ₂ T a_- T = (I + a) Q +

a

+ T

2 - 2 - 2

-which leads to the set of algebraic equations

n+1 _{n+1 n} 1

TV

= U + T V 2 -1 n+1 +

(A + T B) V

1 n+1 1 i n+1 n+1 ₁ _n+1

TCU

T IV I

DV

- T 2

1n

1 n 1 nl n = - C U

+(A -

T B) V - T IV i

DV+

2 401 {I +

-' - - -+ 1 (8) (9) (.5) (6)

Here the solution vector T2n+1 and the vector Vn+1 are both determined implicitly using a linear algebraic solver such as the Newton-Raphson method for a non-linear algebraic system [4].

2.2. Numerical accuracy

The local truncation error vector associated with equation (6) is defined to be = -1 T If) Sr(t+T) -

(I ±

1 T (10)

and using Taylor series, this may be expressed as

1 3 E = - T L3 q(t) + 0(T4) 12 so that 1 3 (3) 1E1 T 2m Now, if Q and Q -(1) -(2) + o(T')

where is the maximum value of the 3rd derivative of every component of q(t) at time t and the numerical method is of 'second order' with error constant

This local truncation error is connected with a single application of the numerical method between t = nT and t = (n+1)T. In integrating m times, say,

from t = 0 to a fixed time T = m T, all the local truncation errors build Up

into the global truncation error defined as [1,3]

(1) ₌ 1 _2 _ L e +

0(T3)

12 -2 where e 2 is independent of T.

-Suppose that the time interval is discretised into 2m subintervals each of width 1

2 T (j = 0,1,...,2m), then the global truncation error is

G(2) = - T2 e + 0(T3) .

48 -2

are the solutions obtained at time T for the two discretisations, then the globally extrapolated solution is (riven by

in = y Q2m + (1 -

y)

Qin -(E) -(2) -(1) 402 (12) (15) q(t) 1 12 i.e. points

m

(13)

and the globally extrapolated error function is

The initial value problem described in equation (2) is of the general non-linear form

Cr(t) = f (t, q(t)) ; q(0) = q -o

and the stability of the adopted numerical method to solve this equation is usually examined with reference to the single model equation [5]

q(t) = X q(t) , q(0) =

qo

In this expression X is taken to be the real part (assumed negative to ensure

a bounded solution) of an eigenvalue of the Jacobian 3f/aq.

For the scalar equation (18) equation (4) simplifies to

q(t + T) = exp (XT) q(t) (19)

and it follows from equation (6) that n+1

(1

+ 0.5 ),

T)

= S, Qn 1 - 0.5 A T

where the amplification factor or symbol

G(E)

= y G

(2) + (1 - y) G(1)

where y is some suitably chosen parameter. In fact, when y = 4/3, then

G(E) = 0(T3) providing a third order numerical method [1]. The introduction of

3m

a third, even finer, discretisation which yields a solution

Q(3) at time T,

2m 3m

and then taking a linear combination of

Q()'

(2) and Q(3), will give a fourth order approximation to q(T).

2.3. Numerical stability

(1/-T) + 0.5

(1/-T) - 0.5

is such that

IS, I 5_ 1 for all values - m < "1--

< 0

lim-1

S1/2

These results imply that there is no restriction on the time steps used in the marching sequence and the method is referred to as being A-stable. The

( = XT 21) 403 (16) (20) (22) =

method is suitable for tackling problems

having oscillating solutions _{[1,3] and} also problems with oscillating and growing solutions. _{In the latter case the} amplification symbol S½ must be allowed to exceed unity _{in modulus so that}

ie+1

exceeds

WI

for n = 0,1,2,...

To this end, let E = 0.5 - _E , where 0

_{< E}

< 0.5 _{is arbitrarily small.} Substituting this value into equation _{(5), the equivalent of equation (7) is} now of the form

E) T a ] Qn+1 2 ₂ 1 _{E) an+1} E) isn+1 [ I + (--2 + E)

T a] Q

+ ( + T CX + + E) T 1 - n 1 - n 1 (23)

which has local truncation _error

E 2 L2 1

= E T q(t) + (

--

)

2 12 T3 L3 q (t) + (24)

and an amplification factor

S 7 E

(1/T)

+ 0.5

+ E

(1/T)

-

0.5 +

The method is reduced to first order with error constant _{C2 = E and possesses a} globally extrapolated error function _{of order T2 if the parameter y = 2.} In the limit T

_-m,

_{the amplification factor}

reduces to (0.5

+ E

) --

_{= -1- 4E}

- 8 E2 - 16 E3

-1/2 - E (-0.5

+ E)

and

-1

-

4c

1-2E

so that the numerical solution _{oscillates and grows.}

When E 0, the previous analysis, i.e. E = 0.5 is _{reiterated and the growth of}

the oscillating solution depends on E _{(>0) which can be arbitrarily}

chosen.

2.4. _{Predictor - corrector combinations}

In order to reduce computation

time, but retaining numerical accuracy,

_a

predictor-corrector combination procedure may be developed. _{The numerical} method devised relies on the following algorithm [1].

404 (25) = (26) S 1 ½

-(a) Predict: The starting point of the sequence relies on the prediction n+1

from the (0,1) Pade approximant method [1]. That is, C = 0 in equation (5) and it may be shown that

n+1

: Q =

ri

+ T a) Qn + T

a

(Qn) + T $n

(b) Evaluate the derivative L Qn+1 from equation (3), i.e.

: L Qn+1 = a Qn+1

n+1 n+1

+

a(Q

) +

Correct Qn+1 from equation (6) by writing

= Qn n+1

C : Q + T (

2n +n+1)

where L Qn+1 is evaluated in (b).

Repeat steps (b) and (c) until each element of successive iteration vectors agree with each other to the appropriate chosen accuracy.

Evaluate the derivative L Qn+1, using the final value of Qn+1 obtained in step (d), a final time using

n+1 n+1

: L Q = a Q +

a(gn+1) +n+1

(f) n: = n+1 ; go to (a)

If step (d) is repeated r times then the algorithm is used in P(EC)rE mode, where r is the number of applications of the corrector required

to correct for convergence. This usually takes 2 or 3 iterations but as r co the stability properties of this predictor-corrector algorithm take on the properties of the corrector [5]. Namely, the (1.1) Pads approximant method used in step (c). Thus the algorithm is second order accurate and A-stable (i.e. unconditionally or absolutely stable), but may be further improved by adopting a globally extrapolated solution as given in equation (5).

2.5. Comment

A stable numerical method has been devised satisfying conditions (i-iv). To illustrate its ability to generate a growing oscillatory solution consider the coupled linear equations

[61.

pl + _{w12 p1} - d cos _{(W1

-W2)0

p2 P2 + w22 1072 = 0 405 Q (c1

e)

2

with initial conditions

101(0) = d/(8 (.02

(wi

- (.02)1, P1(0) = 0 , p2(0) = 1 , D2(0) = 0

where W > 0,

w2 > 0, w1 w2 and d 0 0. The theoretical solution is

1

p (t)

= dt

sin

w

t

+ d cos Ow1-2w2)1 t

1 4w 1

1 8

(A)2 (w1-w2)

p2(t) = cos w t

and pi(t) increases with increasing time. _{Figure 1 illustrates a comparison} between the theoretical and numerical (PECE) solutions _{for d =}

0.5,w1

=3.0rad/s,

U2 = 1.0 rad/s. As can be seen the results generated are accurate.

3. _{Semi-submersible excited by waves}

To predict the behaviour of a _{floating semi-submersible excited by large} amplitude sinusoidal or random waves, it is commonly _{argued that non-linear} effects must be included in the mathematical model if reality is to be represented. _{Therefore the developed equations of motions}

describing the six degrees of freedom of the rigid structure contain non-linear terms associated with inertias, viscous effects, fluid actions and wave forces due to the motions excited in the structure _{or by the waves passing along the}

semi-submersible. _{In this way, the forces and moments acting on the semi-submersible} in its

instantaneous

_{position relative to the wave are determined.}

Such models have been discussed extensively elsewhere and the procedures involved are not reiterated here but accepted without question [9-11].

In this paper, the semi-submersible model (see _{figure 2), the general} mathematical model and equations of motion developed by Naess and Hoff [10] are adopted, but minor changes of detail are introduced, for example, in the method used to evaluate the hydrodynamic coefficients, etc. _{Time domain simulations} of the responses are determined using the (1,1) Pad 4 approximant method _{in the} predictor-corrector combination described previously.

The set of equations describing the behaviour of the semi-submersible _are of the form [9-11]

(27)

(t) = v(t) at -M(t) = E M. (t) + E

-1

i=1 j=1 -d-t- ( t ) = flt) d r m v(t) 1 = F(t) , dt I 1P(t) = M(t) - (t) x I

where the column vectors fl(t), e(t) represent the translations (surge, sway, heave) and the. three Euler angles of rotation respectively. The (3 x 3) matrix

E describes a transformation between an inertial fixed axis system and a body -V

fixed axis system. The (3 x 3) matrices m and I contain the total (body plus -hydrodynamic) mass and total moments and products of inertia, respectively.

The remaining column vectors F(t) and M(t) represent the total external forces and moments acting on the semi-submersible.

The approach developed by Naess and Hoff assumes that the total force and moments are the sum of independent contributions acting on the M pontoons and N columns such that

M N

F(t) = E FY (t) + E F.c(t)

i=1 j=1

1.1,c (t)

--]

Each contribution

FY,

F.c etc may be suitably defined and through these

-i -3

terms non-linearities of the forms relative velocity !relative velocity',

time dependence in the coefficients, etc, are introduced into the analysis. The interested reader may consult Paulling [9], Naess and Hoff [10], etc for the details of the approach and it can be shown that the resultant equations may be cast into the typical forms given in equation (1).

4. Numerical simulations

To illustrate the numerical approach discussed previously a limited series of time domain simulations were performed for the unrestrained semi-submersible model floating in level trim at a draught of 0.25m. A time step interval

T = (0.1 x period of wave) was used in the simulation which was started from rest, allowing the exciting forces and responses to build up gradually and therefore reducing the effect of initial transients.

The semi-submersible was excited into heave, pitch and surge motions by

407

t

individual regular head waves of amplitude a(= 0.02, 0.04,

0.1,

0.15m) having frequencies La in the range 0.35 rad/s. Figure 3 illustrates a

selection of these motion time simulations for various wave amplitude - frequency combinations. For wave amplitude a = 0,04m, the motion amplitudes are nearly symMetric above and below the equilibrium position and

at

illustrated in figure A the predicted response amplitudes are in close agreement with results derived from a linear three, dimensional analysis [12], and the finds of Huang, Hoff and Naess Asitiewave amplitude increases, the non-symmetry characteristic, of the: amplitudes

of

motion becomes more evident and the envelope amolitude

curve displays a small variation with time.

As expected, no mathematical Or numerical difficulty was experienced in the generation of these, time domain responses.

'Conclusions

A numerical procedure, based

on

the (1,1): Pads approximant method, has been developed which is shown, to be mathematically stable and provides accurate

solutions to non-linear differential equations of the form arising in ship and offshore dynamic problems. This has been demonstrated by the generation of an oscillatory growing unstable solution as well as generating realistic simulations of the responses

of

a floating, unrestrained semi-submersible excited by large amplitude regular sinusoidal waves. The latter may be replaced by a randot wave

for without introducing any new concepts into the numerical method. Whether the excitation is deterministic or random, any instability observed in the

simulation

of

the responses, arises, from the developed mathematical model and not ,froth. the numerical approach..

References

(i1) TWIZELL, E.H., PRICE, W.G. and WANG, Y: "Application of Pads. approximants

to the simulation of fluid-structure interactions in irregular seaways".

Int. Conf. in Numerical Methods in Engineering: Theory and Applications. University College, Swansea, 1987.

PRICE, W.G. and BISHOP, R.E.D. Probabilistic theory of Ship dynamics. Chapman Hall, London, 1974.

TWIZELL, EH,.. Computational methods for partial differential equations. Ellis Horwood/Wiley, Chichester, 1984.

A)

BURDEN', R.L. and FAIRS, J..D. Numerical analysis (third edition), Prindle, Weber and Schmidt, Boston, 1985.

w 8.5

[13].

LAMBERT, J.D. Computational methods in ordinary differential equations. Wiley, Chichester, 1973.

CESARI, L. Asymptotic behaviour and stability problems in ordinary differential equations. Springer-Verlag, Berlin, 1959.

GOURLAY, A.R. and MORRIS, J. Ll. "The extrapolation of first order methods for parabolic partial differential equations, II". Siam Journal Numerical Methods, 19, 1980, 641-655.

GOURLAY, A.R. and MORRIS, J. Ll. "Linear combinations of generalized Crank Nicolson schemes". IMA Journal of Numerical Analysis, 1, 1981, 347-357.

PAULLING, J.R. "Time domain simulations of semi-submersible platform motion with application to the tension leg platform". Proc. STAR Symposium, 1977.

NAESS, A. and HOFF, J.R. "Time simulation of the dynamic response of heavily listed semi-submersible platforms in waves". Norwegian Maritime Research 1, 1984, 2-14.

PAULLING, J.R. and SHIN, Y.S. "On the simulation of large amplitude motions of floating offshore structures". Proc. Int. Symposium on Ocean Space Utilization 185, Nihon University, Tokyo, 1, 1985, 235-243. INGLIS, R.B. and PRICE, W.G. "A three dimensional ship motion theory-calculation of wave loading and responses with forward speed". Trans. Roy. Instn. Nay. Arch. 124, 1982, 183-192.

HUANG, X., HOFF, J.R. and NAESS, N. "On the behaviour of semi-submersible platforms at large angles. Offshore Technology Conference, 1982, paper OTC 4246, 193-197.

1-7

00

Fig.l. Comparison between analytical Fig.2. Sketch and particulars of semi-submersible model (---) and numerical (0) growing

[10, 13]. oscillatory solution pi(t) (see

equation (27)) for d=0.5, co1=3.0

rad/s' w2=1.0 rad/s. (6)

(1) 9.98 0.9 s, e SI

;

-2.15 1.1 V IS 000co .ptx C.11 0 0 .L

Fig. 3.

(i)

_(ii)

Fig.4.

Comparison between numerical and experimental _{data for}

(i) _{the heave per unit wave amplitude and (ii) the pitch} per unit wave slope. _{The semi-submersible is in level}

Heave and pitch motion responses excited in the semi-submersible model b: sinusoidal head waves of amplitude a(m) and freque w(rad/s). (i) a = 0.04, W = 4.5; (11) a = 0.10, w = 5.5; (iii) a = 0.15 w = 5.5. (ii .3 3 111.0