Extreme ocean waves. Part I. The practical application of fully nonlinear wave modelling

Applied Ocean Research 34 (2012) 209–224

Contents lists available at ScienceDirect

Applied Ocean Research journal homepage: www.elsevier.com/locate/apor

Extreme ocean waves. Part I. The practical application of fully nonlinear wave modelling W.J.D. Bateman, V. Katsardi, C. Swan ∗ Department of Civil Engineering, Imperial College, London, SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 18 September 2010 Received in revised form 6 May 2011 Accepted 6 May 2011 Available online 2 October 2011 Keywords: Wave modelling Extreme ocean waves Nonlinear waves Wave spectra Directional waves Wave kinematics Wave–wave interactions

a b s t r a c t This paper concerns the description of large surface water waves in realistic ocean spectra and explains how fully nonlinear wave modelling can be applied in practical design calculations. In particular, the proposed methodology incorporates: (i) the transitory nature of individual waves or wave groups arising in a random or irregular sea involving a significant spread of wave energy across the frequency domain; (ii) the short-crestedness of waves arising due to the directionality of the underlying wave components; and (iii) the nonlinearity associated with the evolution of the largest waves. The procedures outlined build directly upon two advances in wave modelling, reported by Bateman et al. [1,2]. By combining these solutions with initial or starting conditions originating from statistical theory describing the most probable shape of a large linear wave, fully nonlinear descriptions of extreme waves in realistic sea states can be achieved. Furthermore, such calculations can be undertaken with relatively modest computational effort and are directly relevant to a wide range of design wave calculations. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction In both offshore and coastal engineering the accurate description of extreme surface water waves is critical to the design of both fixed and floating structures. For example, the setting of deck elevations and the assessment of over-topping rates are largely dependent upon estimates of the maximum crest elevation; the associated water particle kinematics (both velocities and accelerations) are required to determine the applied fluid loads; while the shape or steepness of a wave profile determines the possible occurrence of wave slamming. If these predictions are to be realistic, it is essential that the applied wave model fully incorporates the physical characteristics of the sea state. Field observations confirm that there are three principal features. First, real seas involve a significant spread of wave energy across the frequency domain. This is characterised by the underlying wave spectrum defining a random or irregular sea in which, due to the effects of frequency dispersion, individual waves evolve in both space and time. Second, the wave energy is distributed in terms of its direction of propagation. This again contributes to the evolution of individual waves and ensures that, depending on their directional spread, the waves are short-crested. Third, large individual waves may also be highly nonlinear, their properties differing significantly from the linear sum of the underlying wave components.

∗ Corresponding author. Tel.: +44 2075945999; fax: +44 2075945991. E-mail address: [email protected] (C. Swan). 0141-1187/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.apor.2011.05.002

The starting point for the present paper arises from advances in the numerical modelling of extreme surface waves. Indeed, there presently exists several fully nonlinear wave models appropriate to the description of extreme waves arising in directionally spread seas (see below). For example, taken together [1,2] provide computational procedures appropriate to the description of highly nonlinear surface water waves, involving a significant spread of wave energy in both frequency and direction. Although this combined wave model has been rigorously explained and validated against existing laboratory data, its application to the description of extreme waves in an ocean environment is yet to be considered. The present paper addresses this task and, in so doing, outlines a methodology which is appropriate to the application of a number of related wave models. In particular, it provides a detailed description of the required input conditions, explains how these relate to the evolution of extreme waves in a realistic sea state and outlines the modelling parameters and domain discretisation necessary to optimise the accuracy achieved. With the application of the model explained, a number of sample calculations are provided to demonstrate the practical importance of incorporating both the nonlinearity and the directionality of a specified sea state. In effect, the paper seeks to provide a practical guide concerning the application of a state-of-the-art wave model. The motivation behind this work lies in the increasing awareness of the potential importance of nonlinear effects, both from a regulatory perspective [3] and from recent scientific investigations [4–8]. In considering the practical application of a fully nonlinear, three-dimensional wave model, it will become clear that the

210

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

approach builds upon developments in the specification of extreme waves. Both the A.P.I. and I.S.O. design guidance notes readily acknowledge the benefits of applying deterministic solutions based upon the most probable shape of a large linear wave arising in a specified sea state, thereby avoiding the need for long time-domain simulations. Such solutions were first proposed in [9,10], and popularised in an industrial context in [11]; the latter solution commonly referred to as NewWave, but essentially being identical to that given in [9]. Using this solution to define the initial or starting conditions, well in advance of the focal event, for the wave models outlined in [1,2], the present paper will demonstrate how design waves, with a given return-period, may be modelled in a fully nonlinear sense, taking due account of all the underlying physics. Section 2 commences by highlighting the inadequacies of the existing design solutions and provides a brief review of the present state-of-the-art in terms of nonlinear wave modelling. Section 3 concerns the representation of ocean spectra in the spatial domain. This is an essential step in the successful application of the available nonlinear wave models and builds upon the most probable shape of a large linear wave [9]. Indeed, it is the success of this latter approach that largely facilitates the use of exact wave evolution models based upon the process of time-marching the nonlinear free-surface boundary conditions. Section 4 concerns the practical application of the wave model provided in [1] and demonstrates why computational efficiency becomes fundamentally important when dealing with sea states that are broad-banded in both frequency and direction. Sample calculations involving JONSWAP spectra, with varying directional spreads, are provided in Section 5. The practical implications of these results are discussed in Section 6 together with some concluding remarks concerning the application of the proposed method.

2. Background 2.1. Existing design solutions The commonly applied design wave models may be sub-divided into two categories. The first attempts to simplify the problem by assuming that an individual wave crest may be represented by a wave of permanent form. This approach neglects entirely the underlying wave spectrum and seeks to model a large wave either by calculating an ‘equivalent’ wave height, H, and wave period, T, or by fitting a velocity potential or stream function solution [12] to a measured time-history of the water surface elevation, (t). Irrespective of which method is adopted some nonlinearity can be included but the essential unsteadiness, caused by frequency dispersion, is neglected. Alternatively, the second category neglects the nonlinearity, makes some assumptions about which spectral components are freely propagating, and applies a linear random wave theory. Although this solution includes some unsteadiness the applied linearisation causes additional difficulties, not least those associated with the prediction of the water particle kinematics immediately beneath a wave crest. In this case, the kinematics can only be predicted by applying some form of empirical or stretched wave model, the most widely adopted being that outlined in [13]. Although such models are widely used, they do not satisfy the governing field equation representing mass continuity and irrotationality. Further discussion of this point is given in [14]. In describing these two categories of wave solution, no mention has been made of directionality. Unfortunately, this poses further difficulties, particularly in the case of the nonlinear steady wave models (category 1, above). As a result, directionality is either neglected or introduced in the form of a crude velocity reduction factor, applied uniformly to all frequency components.

2.2. Wave evolution models Given the obvious limitations of the commonly applied wave models, it is perhaps surprising that early wave evolution models have not been more widely adopted. In general, these solutions can again be divided into two categories. The first concerns solutions to the nonlinear Schrödinger equation (see, for example, [15]) or its extension proposed in [16]. Although such solutions can be rapidly calculated, they are applicable to weakly nonlinear waves arising in narrow-banded frequency spectra. As a result, they may be employed to describe the modulations of a regular wave train over a large number of wave cycles, but are not appropriate to the description of extreme waves arising in realistically broad-banded sea states. They will not, therefore, be considered further. The second category of solutions was envisaged in the work of Zakharov [17]; the first detailed calculations being provided in [18]. For uni-directional waves propagating in a two-dimensional fluid assumed to be incompressible, inviscid and irrotational, the velocity potential, , satisfies Laplace’s equation, ∇ 2  = 0, and the kinematic and dynamic free-surface boundary condition applied on z =  can be re-arranged to give t = z − x x , t = −g −

1 2

(1) 2

(x ) + (z )

 2

.

(2)

Within this description, g is the acceleration due to gravity, the subscript denotes partial differentiation with respect to the variable involved, t indicates time and (x, z) define the usual Cartesian coordinates, in which x defines the direction of wave propagation and z is measured vertically upwards from the mean water level. In formulating Eqs. (1) and (2) all time dependence arises on the left hand side. Accordingly, [18] showed that given a spatial representation of  and  at some initial time t = to , (x, to ) and (x, z, to ), Eqs. (1) and (2) can be time-marched allowing the evolution of the wavefield to be calculated for all times. The elegance of this approach lies in its simplicity. From a practical perspective it is important to note that there are no fundamental limitations on the bandwidth of the underlying spectrum. To date numerous examples of this approach have been applied to uni-directional waves, the most important examples of which include: [18–23]. Although such solutions can be highly accurate, they have seldom been applied in engineering design. The reasons for this are three-fold: (i) Even the most efficient formulations require significant computing resources and, until recently, were not easily undertaken on a standard PC. (ii) With the exception of [20], the earlier models noted above were limited to uni-directional wavefields and could not easily be expanded to include the effects of directional spreading. (iii) All such models require initial conditions specified in the spatial domain at some initial time, t = t0 . At first sight this information is seldom available, particularly in the case of field data which typically provides temporal traces at one or more fixed locations ((t) and (t)). Historically, these issues were sufficient to limit the practical application of these wave models. However, recent developments are such that this now needs to be re-considered. First, there has been a relentless improvement in computing speed and resources. This has largely overcome the restrictions noted in (i) above, at least in respect of the most efficient models. Second, a number of fully nonlinear, three-dimensional wave models have been proposed, examples of which include [1,20,24–30]. In considering the practical application of these models, their computational efficiency is of fundamental importance. This is not sought for its own sake but

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

is, in fact, absolutely essential to achieve the resolution necessary to accurately model a wavefield that is broad-banded in both frequency and direction. Given this latter requirement, an accurate wave model must necessarily incorporate a large number of wave components; hence the need for computational efficiency. In the uni-directional wave models, discussed above, the most efficient (and accurate) schemes involve a dimensional reduction so that the entire computation can be undertaken on the water surface, z = . However, with  only defined on z = , the vertical derivative of  present in both Eqs. (1) and (2) cannot be evaluated directly. In the wave models noted above a variety of methods have been applied to overcome this difficulty. Notably, the application of Green’s theory in [18], perturbation methods in [19,21], Cauchy integral theory in [22] and a Dirichlet–Neumann operator in [23]. Unfortunately, several of these approaches are not readily expandable to directional wave models. However, notable exceptions to this include the perturbation methods adopted in the Higher Order Spectral model of [19,21] and the Dirichlet–Neumann operator employed in [23]. Adopting these approaches, fully nonlinear directional wave models have been achieved; examples of the former being provided in the original work of [20] and, more recently, in [24]; while the latter approach has been adopted in [1,25]. Comparisons between these alternative procedures are provided in [31]. In the context of the present study [31] highlights three important points. First, having expressed these solution procedures in a comparable form, it is shown that the differences are minimal. Second, from the perspective of practical applications, it is argued that the recursive formulation of [20] is preferable since it helps to avoid multiples of identical sub-calculations; the same advantage being obtained a posteriori in [1] by their optimised or “quick form”. Third, comments are made concerning the accuracy of the solution procedures, noting that whilst all methods experience difficulties with the convergence of the perturbation and Taylor series expansions, [1] successfully models episodic waves (herein referred to as focused waves) with very high nonlinearity, close to their breaking limit. Schäffer [31] argues that the success of this model in this respect is due to the rapid and local evolution of these wave forms, and hence the lack of time and space for any inaccuracies to contaminate the overall solution. The present authors agree with this view, adding that the recent calculations undertaken in [5], using the model outlined in [1,24], using a solution based upon [20], provide further confirmation of the success of these models when describing this practically important category of extreme wave events. Although there exists a number of fully nonlinear threedimensional wave models, the calculations presented in the sections that follow are based upon the application of [1,2]. The first of these contributions describes the wave evolution model in which an approximation of the Dirchlet-Neumann operator is used to time-march the water surface elevation, (x, y), and the velocity potential on that surface, (x, y, ), using the nonlinear free surface boundary conditions. The second contribution adopts the values of  and  given by the first and uses a related operator to calculate the water particle kinematics arising throughout the flow field. With the present paper principally concerned with a methodology allowing the practical application of a fully nonlinear three-dimensional wave model, a full theoretical description of the model is unnecessary. However, some aspects of the modelling procedure are fundamental to its practical application. Accordingly, a brief summary is provided in Appendix A; Appendix A.1 concerning the wave evolution and Appendix A.2 the kinematics predictions. At this stage it is also relevant to note that the methodology that follows is equally appropriate to the application of other comparable wave models, [20,24,25]. Indeed, given the findings outlined in

211

[31] it is clear that the application of these models should lead to near identical results. 2.3. A methodology for fully nonlinear design calculations Eqs. (1) and (2), together with their three-dimensional equivalents given in Appendix A, confirm that the initial conditions appropriate to the commencement of the nonlinear calculations must involve a spatial description of both  and  at some initial time t = t0 , irrespective of the details of the model employed. This directly affects the form of the proposed solution and the way in which it is applied in practice. It is now widely accepted that linear dispersion and wave component superposition accounts for the evolution of the largest waves in deep water involving realistic, broad-banded, seas. This is not to say that slow wave modulation effects, noted by Hassemann [32,33] and Phillips [34], are not important, rather that the evolution of the most extreme waves occurs very rapidly, involving time-scales of a few tens of wave periods. In such cases the constructive interference of freely propagating wave components, involving the summation of wave crests at one point in space and time, provides the obvious underlying mechanism. Indeed, this assertion forms the basis of the statistical models [9–11] mentioned in Section 1. Although these models are linear, extensive comparisons with field data [35,36], suggest that the underlying methodology is correct. This has immediate implications for the problem in hand. If it is possible to identify the wave frequency components and their relative phasing that contribute to the formation of an extreme event, it is possible to define these components at some initial time, t = t0 , well in advance of the extreme (or focused) event. If this initial time is judicially chosen so that the wave energy is widely dispersed (no large waves are present), a linear wave model will provide an accurate representation of both (x, y, t0 ) and (x, y, , t0 ). In this way it is possible to overcome the difficulties noted in (iii) above. The solution procedure outlined above is central to the alternative design solution we are seeking to provide. In practice, [9–11] describe the most probable shape of a large linear wave. This can be matched to a given design condition by incorporating the correct underlying wave spectrum and scaling the solution to achieve a specified wave height (or crest elevation) corresponding to a given probability of exceedence. By adopting the underlying linear wave components, detailed at an earlier (t = t0 ) dispersed condition, as the input to the nonlinear wave model, a fully nonlinear description of a large design wave can be achieved. If the solution is again matched to a target wave height, the nonlinear change in the associated crest elevation can be predicted. Furthermore, the nonlinear water particle kinematics can be predicted in a fully nonlinear sense, without the ambiguity in the existing solutions outlined in Section 2.1. In adopting this approach it is important to stress that we are not seeking to match the nonlinear calculations to some measured wave profile. In practice, such a profile is seldom available and, if it is, there are other solutions specifically developed for this task [12,37,30]. In contrast, the present procedure seeks to provide a fully nonlinear description of a large design wave arising in a specified sea state and having a given probability of exceedence. In the calculations that follow (Section 5) the input conditions simply define the underlying frequency spectrum (including the spectral peak period), the directional spread and a target or design wave height; the latter based upon a 10−4 annual probability of exceedence. 3. Spatial representation of ocean spectra To apply the wave models outlined above both (x, y, t) and (x, y, , t) must be defined at some initial time, t = t0 , and at discrete spatial locations evenly distributed over the entire computational

212

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

domain (0 ≤ x ≤ (x) , 0 ≤ y ≤ (y) ); (x) and (y) defining some large fundamental wavelengths in the x and y directions over which the solution is assumed periodic. In normal circumstances, particularly those concerning field data, an evolving wavefield is recorded as a time-history, (t), at a single point fixed in space. As a result, it is characterised by a frequency spectrum, S (ω), where ω = 2/T. Converting this to a spatial representation, S (k), where k is the wave number, is relatively straightforward provided it is assumed that all the wave components are freely propagating and hence satisfy ωi2 = gki tanh(ki h),

(4)

where (t) is the autocovariance function 1

(t) =



∞

S (ω) cos(ωt) dω,

2

(5)

0

and the variance is defined by



∞

2 =

S (ω) dω.

(x, t) =

(3)

where h is the constant water depth and the index i denotes a particular wave component. Unfortunately, this process is not in itself sufficient since the adopted Fourier series for both  and  also require the surface values to be periodic in both x and y (see Appendix A.1). Achieving spatial representations of the surface parameters that are both periodic and representative of the underlying frequency spectrum is altogether more difficult. We have already noted that the model application will build upon the statistical predictions of [9], referred to in [11] as the NewWave model. In its original form this model describes the expected or most probable shape of the water surface elevation, (t), conditional on the presence of a large crest elevation,  = A, at some arbitrary time t = . Accordingly, if ()t= = A, (t )t= = 0 and (tt )t= < 0, where the subscript t denotes a partial derivative with respect to time, the most probable shape for the water surface elevation is given by (t) = A(t),

about max = A and which is approximately periodic. To enhance this periodicity, the ends of the trace may be smoothed to ensure they converge to the same value, typically zero. Having identified an appropriate segment of the trace, Fourier analysis is used to compute the amplitudes of the underlying wavenumber components, ân . Re-applying Eq. (3) the frequencies, ωn , corresponding to the uniformly spaced wavenumber components, kn , can be determined and the resulting water surface elevation given by

(6)

0

In effect, Eq. (4) provides a linear realisation of the extreme wave we are seeking to describe nonlinearly using the wave model outlined in Appendix A.1.

ˆ N 

aˆn cos(kn x − ωn t).

(9)

n=0

Fig. 1 outlines the various stages in this analysis and contrasts its success with other methods outlined below. These results are based on a JONSWAP spectrum [38]. This is widely applied in design calculations and forms the basis of the sample calculations given in Section 5. Within this description ˛g 2 S (ω) = 5 exp ω



−ˇ

ωp4 ω4







exp

−(ω−ωp )2 /2ω2 2 p



,

(10)

with = 0.07, = 0.09,

ω ≤ ωp ω > ωp ,

where ˛ = 0.0081, ˇ = 1.25, ωp is the frequency of the spectral peak or 2/Tp where Tp is the spectral peak period and  is the peak enhancement factor, which is taken to be 2.0. Fig. 1(a) depicts the JONSWAP spectrum S (ω), while Fig. 1(b) presents the most probable shape of a large linear wave ((t) from Eq. (4)) non-dimensionalised with respect to A. Fig. 1(c) concerns the corresponding spatial profile, (x), and identifies the segment of the trace, of length (x) , to which the Fourier analysis is applied. Fig. 1(d) depicts the outcome of the analysis in terms of an amplitude spectrum of the discrete wavenumber components, ân (kn ). In Fig. 1(b)–(d) the results of the present approach are labelled as being based on Fourier analysis. An alternative but closely related procedure adopted in [26] employs a least squares fit (LSF) to determine the amplitudes of the wavenumber components ân , consistent with a good fit to the timehistory, (t), defined in Eq. (4). Setting x = 0 in Eq. (9) and matching it to Eq. (4) gives ˆ N 

3.1. Uni-directional simulations

aˆn cos(ωn t) = A(t).

(11)

n=0

Adopting the linearly predicted wave profile given in Eq. (4), the amplitudes, ai , of the frequency components, ωi , may be identified by Fourier analysis such that ai =

A 2

S (ωi ) ω,

(7)

where ωi = i ω. Assuming that the linear focal event corresponding to a crest elevation of  = A occurs at t = 0 (i.e. setting  = 0), the corresponding spatial representation of the water surface elevation, (x) at t = 0 is given by (x) =



ai cos(ki x),

(8)

i

where the wavenumber components, ki , are based on the evenly distributed frequencies, ωi , and are derived from the dispersion equation given in Eq. (3). Unfortunately, the spatial distribution given in Eq. (8) is not periodic. To overcome this difficulty a number of different approaches have been applied. The first approach selects from Eq. (8) a section of the trace with the desired fundamental wavelength, (x) = 2/k0 , which is centred

This approach has the benefit of being able to limit the creation of high wavenumber components, while retaining the best fit to the desired (t) (Eq. (4)) using available terms. This is particularly relevant to the reproduction of experimental data where the underlying spectrum is inevitably truncated due to the limitations of the laboratory wave generation system [39,26]. Given the earlier discussion concerning computational efficiency, a method that limits the total number of wave components has obvious advantages. Furthermore, the unnecessary inclusion of high wavenumber components, having little or no energy, can impact adversely on the stability of even the most robust numerical schemes. The results produced by this method and denoted by LSF are also indicated in Fig. 1(b)–(d). In both this and the previous method (Fourier analysis), the description of the desired wave profile is good. However, away from the central maximum the spatial records, (x) in Fig. 1(c), show some significant and unrealistic high-frequency fluctuations; the small insert providing a close up of the region of interest. Further evidence of this is seen in the amplitude spectra, ân (kn ) in Fig. 1(d), where the tail of the spectrum does not vary as smoothly as one might expect; the magnitude of the effect again being highlighted

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

a

b

10

213

1.2 NewWave

8

Fourier

0.8

η(t) / A

Sηη (ω)

LSF 6 4 2

Continuous

0.4 0.0

-0.4

0 0.0

0.5

1.0

ω (rad/s)

1.5

2.0

2.5

-0.8 -100

-80

-60

-40

-20

0

20

40

60

80

100

t (s)

c

d 1.2 1.0

0.05

Fourier

Fourier 1.0

LSF Continuous

0.6

an(kn) / a max

η(x) / A

0.00

/ NewWave -0.05 -1000

-800

-600

0.2

Continuous

0.8

/ NewWave

0.6 0.4 0.0 0.05

-0.2 0.2 -0.6 -1000

-800

-600

-400

-200

0

200

400

600

800

1000

LSF

0.2

0.0 0.00

0.05

0.10

0.10

0.15

0.15

0.20

0.25

kn (rad/m)

x (m)

Fig. 1. Representing an ocean spectrum in the spatial domain (a) JONSWAP spectrum, S (ω), with Tp = 11.4 s and  = 2.0, (b) (t), (c) (x), (d) amplitude spectrum in space, â(k).

in the insert. Given the inevitable concerns regarding numerical stability, this effect is clearly undesirable. To overcome this difficulty a third approach is proposed. Instead of seeking to model (t) based on Eq. (4), an alternative approach is sought that directly transforms the power spectrum in the frequency domain into a power spectrum in the wavenumber domain. The essence of this transformation is given in Fig. 2. In adopting such an approach (ω(k + (1/2)ık) − ω(k − (1/2)ık)) , ık

wave crest at x = , i.e. ()x= = A, (t )x= = 0 and (tt )x= < 0. The solution is similar to that presented in Eqs. (4)–(6), except that x and k replace t and ω respectively. The resulting surface elevation, corresponding to the most probable shape of a large wave in space, is now periodic. As a result, discrete amplitudes can be calculated by taking a Fourier transform of (x). However, given that the autocovariance function included in Eq. (4) essentially contains an inverse Fourier transform, the amplitude spectrum is simply given by

(12)

ˆ = a(k)

where ık is the discretisation of the new spectrum in the wavenumber domain, such that k = nık. Taking the limit as ık → 0, Eq. (12) yields

where

S (k) = S (ω(k))

∂ω , S (k) = S (ω(k)) (13) ∂k which defines a continuous spectrum. A modified version of the theory outlined in [9,10,11] can now be applied to this solution, based on the conditional presence of a large

A ˆ n2



S (k),

(14)

S (k)dk.

(15)

∞

ˆ n2 = 0

The required discrete wavenumber spectrum is now formed by integrating Eq. (14) over the discretisation width, k0 , centred about each discrete wavenumber component



(n+(1/2))k0

aˆn = 80

Sηη (k)

60

10

Sηη (ω)

40 20

8

0

6

0

0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2

k (rad/m)

4 2 0 0

0.2

0.4

0.6

0.8

ω (rad/s)

ˆ a(k)dk.

(16)

(n−(1/2))k0

1

1.2

1.4

Fig. 2. Schematic diagram showing transformation from S (ω) to S (k).

With ân defined, Eq. (9) yields both the temporal and the spatial profiles of the water surface elevation. The results achieved using this approach are referred to as the continuous method and are again presented in Fig. 1(b)–(d). Contrasting the different approaches presented in Fig. 1 suggests that all three methods (Fourier, LSF and Continuous) provide a good description of (t) predicted by Eq. (4) (Fig. 1(b)). However, in the case of the first two methods this is to be expected since they are based upon an initial fit to this record. In contrast, the continuous method provides a good description of (t) and, more significantly, a substantial improvement in the description of (x) (Fig. 1(c)). In particular, there are no high-frequency oscillations present at the ends of the record (arising far from the maximum crest elevation) and the corresponding amplitude spectrum, ân (kn )

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

in Fig. 1(d), is smooth in the tail of the spectrum. Indeed, in Fig. 1(c) and (d) the continuous method is identical to the spatial prediction of the modified linear prediction on which it is based. The continuous method also has the advantage that it is only dependent upon the initial choice of k0 (see below). In contrast, the LSF method is dependent upon the length of the temporal profile selected and the choice of the components included. Similarly, the Fourier transform method is also subject to uncertainties depending on the extent of the smoothing applied to the ends of the (x) record. In considering these approaches it is important to note that if a very large number of wave components were included, the difficulties associated with the Fourier and LSF methods would reduce, and in the limit they would disappear. However, given the practicality of applying a fully nonlinear wave model, a substantial increase in the number of wave components included within the solution is not feasible, particularly since the directional spread is yet to be considered (see below). Furthermore, whilst the Fourier method provides a logical first approach, and the LSF method has been previously applied [26], the continuous method gives the best results. Given that the success of any subsequent wave calculation is dependent upon the accuracy of the initial or starting conditions specified in the spatial domain, (x, y, t0 ) and (x, y, , t0 ), this is an important aspect of the model application.

a 5

4

m=l / l 0

214

3

2

1

0 0

1

2

3

4

5

6

7

8

9

10

n=k / k 0

b

3.2. Introducing directionality It is common practice in analytical and statistical descriptions of directionally spread wavefields to assume that the sea state may be represented by a frequency spectrum coupled with a frequency independent directional weighting function. Two widely applied examples of this weighting include the Mitsuyasu [40] spreading parameter, s, whereby

   2

(17)

and the application of a normal distribution having a standard deviation of  so that



E  =

ˇ exp 



− 2 2 2



,

(18)

where E() defines the proportion of the total energy propagating at an angle  to the mean wave direction and ˇ is a normalising coefficient dependant upon the model applied. If ˛() defines the proportion

of the total amplitude propagating in a given direction, ˛()˛ E() and the constant of proportionality simply ensures that the sum over all directions gives the total amplitude associated with any individual frequency (or wave number) component. To put these distributions into perspective, field data recorded at the Tern Platform in the northern North Sea [36] suggest that the directional distribution associated with severe storms is typically represented by a normal distribution with a standard deviation of  = 30◦ . This corresponds closely to a Mitsuyasu distribution with a spreading parameter of s = 7. Irrespective of which weighting function is applied, the first difficulty to overcome involves the representation of a uni-directional wave spectrum propagating at any angle  to the x-axis, where the latter is taken to define the mean wave direction. If, for arguments sake, we set  = 20◦ , it is clear from Fig. 3(a) that within a discrete wavenumber spectrum (k = nk0 , l = ml0 , where k and l are the wave numbers in the x and y directions respectively such that k0 = 2/(x) and l0 = 2/(y) ) there is no sequential line of components that travel at exactly  = 20◦ . Indeed, this condition only exists for  = 0◦ , 45◦ and 90◦ . Under such circumstances an obvious approach is to apportion an appropriate amount of energy to the components that lie closest to where the energy should ideally reside. This process is

c 50

Angular Resolution ∆θ (deg)



E  = ˇ cos2s

40 30 20 10 0 0

5

10

15

20

25

30

35

40

45

50

K / K0 Fig. 3. Incorporating a directional spread. (a) Diagram of a discrete 3-D spectral grid. The thick solid line indicates the ideal location for all energy consistent with a wave travelling at 20◦ to the x-axis. The solid grid lines represent the edges of discrete cells, the centres of which are defined by the intersection of the dashed lines. These centres represent the discrete positions, in wavenumber space, to which energy may be ascribed. (b) 3-D view of the amplitude distributions in wavenumber space corresponding to a uni-directional JONSWAP spectrum propagating at 20◦ to the x-axis. (c) Variation in angular resolution with radial distance.

highlighted by the small arrows in Fig. 3(a) that indicate the locations on the underlying wavenumber grid where the wave energy is to be located. The outcome of this approach is indicated in Fig. 3(b). This provides a spectral view of the placement of wave amplitudes propagating at 20◦ to the x-axis. This process naturally results in a collection of waves travelling on slightly different directions. However, preliminary calculations presented in [1] confirm that this

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

approach allows uni-directional waves propagating at any angle across the computational domain to be successfully modelled. In many respects the process noted above is akin to the discretisation of a continuous spectrum â(k) into ân . The principle difference, however, is that the angular resolution (based on the angle subtended by adjacent centroids) is non-uniform across the wavenumber domain. Indeed, it increases rapidly with distance from the origin n = m = 0. As a result, it becomes imperative when modelling a directionally spread sea that the dominant wave energy is positioned away from the origin. Fortunately, only a small increase in the radial distance is needed to achieve a dramatic improvement in resolution as indicated in Fig. 3(c). The mechanism for distributing the energy as described above is comparatively simple. In Fig. 3(a) there are a number of open dots along the 20◦ line. These represent the intersection between the line and the underlying grid. The radial distance between successive dots is taken to represent the section of wave amplitude from a unidirectional spectrum to be placed at the centre of a

specific grid.

aˆ (r) dr,

(19)

Ka

where the uni-directional spectrum aˆ (r) is defined in Eq. (14). To create a spectrum representing a directionally spread sea the process described above is repeated several thousand times, typically 100 iterations per degree, through an appropriate range of angles. At this stage the amplitude factor, ˛() from either Eq. (17) or Eq. (18), is applied to weight the individual amplitudes as they are added to the overall spectrum. An example of this process is given in Fig. 4(a). This describes a three-dimensional spectrum corresponding to the most probable shape of a large linear wave ((t) in Eq. (4)) arising in a JONSWAP spectrum with  = 2.0, Tp = 11.4 s and a normal spread of  = 30◦ applied over the range −90◦ ≤  ≤ + 90◦ . Fig. 4(b) and (c) also concern this case and demonstrate respectively that: (i) At any arbitrary angle , in this case  = 20◦ , the required form of the normalised spectrum is recovered, and (ii) the overall directional spread corresponds closely to the chosen distribution.

1.0 Uni-Directional

+ l2 ,

Kb

aˆnm =

b

Directional Slice

0.8

a(k) / a max



a

0.6

0.4

0.2

0.0 0.00

0.02

0.04

0.06

0.08

0.10

K (rad/m)

c

1.0 Target Actual

0.8

a(θ) / a max

the With the radial distance to any point defined by K = wave amplitude that should be placed in the grid square spanned by two successive points Ka and Kb , is given by k2

215

0.6

0.4

0.2

4. Model application 0.0

4.1. Initial conditions Analytical models that capture the essential unsteadiness of a wavefield are limited to (i) a linear random wave theory and (ii) a weakly nonlinear second-order model; the latter based upon the sum of the two wave interactions first identified in [41]. Given this constraint, the initial conditions representing the required wavefield must be specified well in advance of the focal event so that the wave energy is fully dispersed and, as a result, there are no large waves present within the computational domain. The methodology underlying this process is highlighted in Fig. 5. These calculations concern a JONSWAP spectrum with Tp = 11.4 s and  = 2.0, focused to produce a large linear wave, (t) in Eq. (4), with max = 7.8 m at the focal location. Fig. 5(a) provides a linear description of the maximum crest elevation occurring anywhere across the computational domain for times prior to the focal event, at t = 0. At this point we will only consider those results relating to a uni-directional wavefield, the data for which is indicated by the solid line. Using figures of this form, it could be argued that the initial time should be chosen so as to minimise the maximum crest elevation, or at least a time at which the minimum value is asymptotically approached (t ≤ − 2000 s in Fig. 5(a)). However, this approach leads to unnecessarily long computations, with the

-80

-60

-40

-20

0

20

40

60

80

θ (deg) Fig. 4. A directionally spread JONSWAP spectrum. (a) 3-D contour plot of amplitude distribution in wavenumber space, â(k, l), (b) Wavenumber spectrum, â(k), for  = 20◦ , (c) directional spectrum, â(), for k = 0.03 (rad/m).

possibility of reduced accuracy due to the cumulative numerical errors arising from a very long time-stepping procedure. A preferred approach, pioneered in [26], is to consider the magnitude of the second-order correction when describing the largest wave crest anywhere within the computational domain. After repeated studies of several laboratory wave cases [26] suggests that provided this second-order correction is less than 2%, a linear wave model is sufficient to describe the initial conditions, (x, t = to ) and (x, , t = to ). Calculations of this type, relating to the previous JONSWAP spectrum, are outlined in Fig. 5(b). In applying this approach, Fig. 5(b) suggests that to = − 2000 s would still be appropriate. Following the developments outlined in [1], and the potential to incorporate many more interacting wavenumber components, additional sensitivity studies suggest an improved criteria in which provided the second-order correction is not larger that

216

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

a 12 Uni-directional Directional

η max (m)

9

6

3

0 -2000

-1000

-800

-600

-400

-200

0

t (s) Most Linear

Nonlinear

b 2nd Order Correction (%)

16 Uni-directional 12

Directional

8

4 Second Order Linear 0 -2000

-1000

-800

-600

-400

-200

0

t (s)

c

Fig. 6. Spatial profile of the maximum water surface elevation, (x, y). (a) Unidirectional (  = 0), (b) long-crested (  = 10◦ ), and (c) short-crested (  = 30◦ ).

10

Start times -800s -600s -400s

η(t) (m)

6

2

-2

-6 -90

-80

-70

-60

-50

-40

-30

t (s) Fig. 5. Defining the initial conditions (a) max vs. t. (b) 2nd-order correction to max vs. t. (c) (t) for different start times.

required to time-march the entire surface over several peak periods using the fully nonlinear wave model outlined in [1]. Fortunately, in directional wavefields the additional dispersion due to directionality is significant. As a result, the water surface and the velocity potential on the surface is well represented by a linear wave model within only a few peak periods of the focal event. Evidence of this is seen in Fig. 5(a), where the maximum crest elevation associated with the directional wave case reduces very rapidly. In undertaking directional wave calculations the start-time appropriate to the equivalent uni-directional wave case clearly provides an upper limit for to . However, the preferred approach is to undertake a number of simulations to confirm that the final result, usually (t) at the focal location, is independent of the chosen initial conditions. 4.2. Model parameters

3%, a second-order solution could be used to generate the initial conditions. This condition corresponds to to = − 800 s in Fig. 5(b). To illustrate the effect of applying these criteria and changing the initial start time, to , Fig. 5(c) describes the focal event generated with initial conditions specified at to = − 800 s, −600 s and −400 s. These cases are all uni-directional and correspond to the JONSWAP spectrum described earlier. From a practical perspective there is virtually no difference in the results from the two earliest times. Indeed, even the solution based upon the latest start time, to = − 400 s, is fairly consistent, despite the fact that in this case the second-order correction is in excess of 4%. The application of this procedure to uni-directional wave trains is, as indicated above, comparatively straightforward. Unfortunately, its application to directional wavefields is more difficult. In these cases the available second-order analytical models, [42] in infinite depth and [43] in finite depth, require a surprisingly large computational effort. For example, in a typical multi-directional problem the number of surface points might be Nx × Ny = 256 × 256 (see below). In this case the analytical models require computations of order (Nx × Ny )2 to compute the water surface elevation, (x, y). This level of computation would typically exceed the effort

Having repeatedly stressed the need for good resolution to successfully model realistic broad-banded spectra, guidance is yet to be given on the appropriate number of surface points (Nx × Ny ) or their physical spacing necessary to optimise the corresponding fundamental wavenumbers (k0 , l0 ). In several respects the available options are limited by the need to maximise computational efficiency and maintain realistic calculation times. For the wave model outlined in [1] and briefly summarised in Appendix A.1, the majority of the computational effort is associated with the large number of Fourier transforms that need to be applied. For these to be undertaken efficiently, both Nx and Ny should be an integer power of two. Assuming that fast Fourier transforms are first applied in the x-direction, followed by the y-direction, the resulting ˆ x and N ˆ y wavenumber components in the x and spectrum contains N ˆ x = Nx /2 + 1 and N ˆ y = Ny + 1. In y directions respectively, where N ˆ y includes the latter case it should be noted that the dimension N components travelling in both the positive and the negative y directions. In the case of a uni-directional wavefield, sample calculations suggest that for a realistic broad-banded frequency spectrum

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

217

Fig. 7. The evolution of a large short-crested wave event in a JONSWAP spectrum, Hmax = 30.2 m, Tp = 16.6 s,  = 2.5 and  = 30◦ .

256 ≤ Nx ≤ 2048 and Ny = 1, where the lower limit for Nx may be reduced further depending on the steepness of the wave profile. As directionality is increased, Ny must obviously increase and Nx may be reduced. Considering calculations involving a given computational effort, possible sets of values for (Nx , Ny ) are (32768, 2), (16384, 4), (8192, 8), (4096, 16), (2048, 32), (1024, 64), (512, 128) and (256, 256). Of these, only the last two could sensibly be applied in a directional wave simulation. In general terms the resolution in the physical domain must be sufficient to capture the detail of an evolving wavefield, particularly the steepness of any part of the wave profile. Conversely, too many wavenumber components, each associated with very small contributions to the overall solution greatly increases the risk of numerical errors and may lead to an unstable solution. In the present paper all the uni-directional results are based upon Nx = 1024, while the multidirectional results employ Nx = 512 and Ny = 256. ˆ x, N ˆ y ) defined, the fundamental wavenumbers (k0 , l0 ) With (N are crucial to the description of the underlying spectrum. In effect, these parameters should be sufficiently small to ensure that there are enough wavenumber components to describe the most energetic parts of the spectrum, but at the same time they must be sufficiently large to ensure that the high wavenumber components can effectively represent the energy transfers associated with the nonlinear wave–wave interactions. This latter point is particularly important for the description of the largest waves. If the high wavenumber components are not sufficiently high, the effective

bandwidth will be restricted and wave energy may be lost due to the truncation of the wavenumber spectrum. Such errors generally lead to an under-estimation of the properties of the largest waves. To overcome this difficulty the weighted-mean position of the spectrum in either coordinate direction should be approximately 1/10 th the size of the domain so that



2

aˆn nk0

n



2 aˆn

≈

1 ˆx k0 N 10

(20)

≈

1 ˆ y. l0 N 10

(21)

n

˙m aˆ 2m ml0 ˙m aˆ 2m

Having adopted these definitions, the calculation procedure requires the selection of an appropriate time-step, t. Again, there is a balance to be achieved in choosing this value. If t is too large, the complexity of the wavefield will not be modelled and cumulative numerical errors will result. However, if t is too small the calculation becomes excessively long. Given the nature of the governing equations (Appendix A.1), it is clear that the optimum value of t is dependent upon the nonlinearity of the wavefield, particularly the steepness of individual waves. Within the present calculations the majority of the time-marching was undertaken using a fourth-order Adams–Bashford predictor coupled with a fifth-order Adams–Moulton corrector. Using this scheme repeated

218

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

20

a Linear Nonlinear

η [m]

10 0 −10 −20

−10

0 time [s]

10

20

b 2

du/dt [m/s ]

u [m/s]

12 8 4 0 −4 −20

−10

0 time [s]

10

20

6 4 2 0 −2 −4 −6

2

du/dt [m/s ]

u [m/s]

6

−20

−10

0 time [s]

10

20

d 2

du/dt [m/s ]

u [m/s]

6 4 2 0 −2 −4 −6

6

4 2 0 −2 −4 −6

−20

−10

0 time [s]

10

20

−20

−10

0 time [s]

10

20

−10

0 time [s]

10

20

f

c 6 4 2 0 −2 −4 −6

e

g

4 2 0 −2 −4

−20

−10

0 time [s]

10

20

−6

−20

Fig. 8. Properties of a large wave arising in a uni-directional JONSWAP spectrum (Hmax = 30.2 m, Tp = 16.6 s,  = 2.5 and  = 0). (a) (t), (b) u(t) on z = (t), (c) u(t) on z = 0, (d) u(t) on z = − 12 m, (e) du/dt(t) on z = (t), (f) du/dt(t) on z = 0, (g) du/dt(t) on z = − 12 m.

calculations involving a large number of realistic ocean spectra and different time-steps suggest that highly accurate and stable results may be achieved by setting

t ≤

1 Tp , 200

(22)

where Tp is the spectral peak period. Finally, having defined all the model parameters, it is necessary to note that with calculations undertaken on a single 2.4 GHz processor, Nx = 512, Ny = 256 and t = 0.05 s, the simulations of a near breaking wave, involving 17000 time steps, takes approximately 20 h of CPU. 5. Numerical calculations This section provides calculations relating to three large design wave events, the purpose of the results being to demonstrate the successful application of the wave model and to highlight both the importance of directionality and nonlinearity. All three cases are propagating in a constant water depth of h = 150 m and arise in a JONSWAP spectrum having a peak period of Tp = 16.6 s and a

peak enhancement factor of  = 2.5. The first example concerns a uni-directional wave train (s =∞ or  = 0), while the latter two correspond to directionally spread wavefields in which the standard deviation of the normal distribution, Eq. (18), is given by  = 10◦ and 30◦ respectively. In considering these cases, the second corresponds to a relatively long-crested sea state, while the third is short-crested. In each of these three wave cases the maximum wave height is defined by Hmax = 30.2 m, giving a nominal wave steepness of (1/2)Hmax kp = 0.22. The surface profiles, (x, y), corresponding to the maximum crest elevation associated with each of the three wave events are presented in Fig. 6; the difference between the records highlighting the influence of directionality and the development of progressively shorter crest lengths. The evolution of the short-crested wave event (  = 30◦ ) is considered in Fig. 7. This provides a series of snap-shots describing (x, y) taken at specific instances in time, up to and including the occurrence of the largest wave crest. Comparisons between these profiles show the progressive convergence of wave energy leading to the occurrence of a large, near-focused, wave event at t = 0.2 s in Fig. 7(h). Although Figs. 6 and 7 provide

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

219

a Linear

20

Nonlinear

η [m]

10 0 −10 −20

12

−10

0 time [s]

10

20

b du/dt [m/s2]

u [m/s]

8 4 0 −4 −20

−10

0 time [s]

10

20

6 4 2 0 −2 −4 −6

2

du/dt [m/s ]

−10

0 time [s]

10

20

d 2

du/dt [m/s ]

u [m/s] u [m/s]

6 4 2 0 −2 −4 −6

−20

−20

−20

−10

0 time [s]

10

20

−20

−10

0 time [s]

10

20

−10

0 time [s]

10

20

f

c 6 4 2 0 −2 −4 −6

e

−10

0 time [s]

10

20

6 4 2 0 −2 −4 −6

6 4 2 0 −2 −4 −6

g

−20

Fig. 9. Properties of a large wave arising in a directional JONSWAP spectrum (Hmax = 30.2 m, Tp = 16.6 s,  = 2.5 and  = 10◦ ). (a) (t), (b) u(t) on z = (t), (c) u(t) on z = 0, (d) u(t) on z = − 12 m, (e) du/dt(t) on z = (t), (f) du/dt(t) on z = 0, (g) du/dt(t) on z = − 12 m.

an overview of the evolving wavefields, they do not allow rigorous comparisons. These are provided in Figs. 8–10 and involve time-histories of the water surface elevation, (t), the horizontal component of the water particle velocity, u(t), and the horizontal component of the unsteady water particle acceleration, ∂u/∂t(t). In each case the time-histories relate to conditions at the focal location, or the location of the largest wave crest, and kinematics data is provided at three vertical elevations. By comparing the results of the fully nonlinear computations with equivalent linear calculations based upon the maximum wave height (Hmax ), the importance of the wave nonlinearity is clearly identified. Furthermore, comparisons between the three wave cases allow the importance of directionality to be assessed. Fig. 8 concerns the steep uni-directional wave case. Fig. 8(a) describes (t) and contrasts the results of the fully nonlinear model with linear predictions (Eq. (4)) corresponding to a wave of equal height; the latter based upon a zero up-crossing analysis. To facilitate comparisons between these data sets, the time-base associated with the nonlinear calculations has been shifted so that the largest wave crest arises at the linearly predicted focus time, t = 0; a similar approach having been adopted on all subsequent plots. In this case the importance of nonlinearity is clearly demonstrated in

respect of the increased crest-trough asymmetry and hence the larger maximum crest elevation. Indeed, the nonlinearity accounts for an additional 2.0 m of crest elevation, representing an increase of approximately 11.5% over the linearly predicted value. Fig. 8(b)–(d) concern u(t) at z = (t), z = 0 (or the mean water level, MWL) and z = − 12 m respectively. At the second elevation (z = 0 in Fig. 8(c)), the velocity trace is intermittent reflecting the passage of a wave trough where (t) < 0. At each elevation the model predictions are compared to the linearly predicted kinematics arising beneath the linear wave profile (Eq. (4)) of equal height (Fig. 8(a)). In calculating the linearly predicted kinematics the well-known problems associated with high-frequency contamination are overcome by adopting a stretched or empirical wave model based on [13]. Although this procedure is commonly adopted in design practice, comparisons with the fully nonlinear model clearly demonstrate the inadequacy of this approach. Indeed, if one considers the maximum velocity arising at the water surface the nonlinear effects lead to an increase of 51% (or 4.15 m/s) above the linearly predicted value. Although the difference between the linear and nonlinear results reduces with increasing distance beneath the water surface, the nature of the nonlinear interactions are such that they will not necessarily converge. For example, at z = 0 the nonlinear increase

220

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

20

a Linear Nonlinear

η [m]

10 0 −10

6

0 time [s]

10

20

b

6 2

8 6 4 2 0 −2 −4

−10

du/dt [m/s ]

u [m/s]

−20

−20

−10

0 time [s]

10

c du/dt [m/s2]

u [m/s]

0 −2

6

−20

−10

0 time [s]

10

d 2

du/dt [m/s ]

u [m/s]

0 −2 −10

0 time [s]

10

20

−20

−10

0 time [s]

10

20

−20

−10

0 time [s]

10

20

−10

0 time [s]

10

20

f

4 2 0 −2

6

2

−20

0 −2

−4

20

4

−4

2

6

2

−4

4

−4

20

4

e

g

4 2 0 −2 −4

−20

Fig. 10. Properties of a large wave arising in a directional JONSWAP spectrum (Hmax = 30.2 m, Tp = 16.6 s,  = 2.5 and  = 30◦ ). (a) (t), (b) u(t) on z = (t), (c) u(t) on z = 0, (d) u(t) on z = − 12 m, (e) du/dt(t) on z = (t), (f) du/dt(t) on z = 0, (g) du/dt(t) on z = − 12 m.

in the maximum horizontal velocity is 13% (or 0.73 m/s), while at z = − 12 m this increase reduces to 3% (or 0.12 m/s). Nevertheless, significant differences in the velocity profile, u(t), remain. A similar sequence of results concerning ∂u/∂t(t) at z = (t), z = 0 and z = − 12 m are presented in Fig. 8(e)–(g). The data presented on these figures concern time-histories of the unsteady component of the horizontal acceleration beneath the large wave event. Given that the water particle acceleration is strongly dependent upon the steepness of the wave profile, the (t) traces presented in Fig. 8(a) suggest that the nonlinearity will significantly increase the maximum acceleration and that it will occur at a higher elevation above the MWL. Fig. 8(e) confirms both of these trends. First, the nonlinearity of the wavefield increases the linearly predicted maximum acceleration by 111.5%. Second, the acceleration maximum (both positive and negative) move closer to the focal position. As a result, there is a very rapid reversal from the maximum positive to the maximum negative acceleration. A similar trend, albeit of reduced magnitude, is also observed at lower elevations, z = 0 in Fig. 8(f) and z = − 12 m in Fig. 8(g).

Figs. 9 and 10 provides a similar sequence of results relating to the long-crested and the short-crested directionally spread wavefields respectively. In these cases the underlying spectrum, S (ω), is as defined above and the directional spread (  = 10◦ and 30◦ respectively) applied uniformly to all frequency components. It has already been noted (Section 3.2) that the directional spread adopted in the short-crested case,  = 30◦ , is believed to be representative of severe storm conditions arising in the northern North Sea [36]. In both directional cases, the maximum wave height was held constant at Hmax = 30.2 m. In considering this value, it is interesting to note that in the uni-directional wave case (Fig. 8) Hmax = 30.2 m represented a near-limiting wave case. However, the two directionally spread wave cases were not at their breaking limit; an increase in the input amplitude sum, ân , leading to substantially larger waves before the onset of wave breaking. In part, this provides the first evidence of the importance of directionality. If the wavefield has even a small directional spread, the nonlinearity rapidly reduces, allowing larger limiting wave heights. This observation was first reported by Johannessen and Swan [39] in

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

their analysis of experimental data based upon idealised laboratory spectra. However, their explanation – that irrespective of wave height an increase in directionality leads to reduced in-line wavefront steepness (even in a linear sense) – is valid for all underlying spectra. Indeed, the present calculations are consistent with recent contributions included in [5,6,29]; all emphasising the importance of directionality when seeking to describe large waves arising in realistic ocean spectra. Figs. 9(a) and 10(a) concern (t) and again contrast the nonlinear predictions with linear calculations (Eq. (4)) based on a wave of equivalent height. The difference between these predictions are similar to those observed in the uni-directional case (Fig. 8(a)); the nonlinearity leading to a significant enhancement of the cresttrough asymmetry and hence larger maximum crest elevations. However, with the introduction of directionality, and a corresponding reduction in the local wave steepness (∂/∂x), the magnitude of the nonlinear change is reduced. For example, in the longer crested case (  = 10◦ in Fig. 9) the maximum crest elevation is increased by 1.78 m or 10% of the linearly predicted value. In contrast, in the short-crested case (  = 30◦ in Fig. 10), the maximum crest elevation increases by 1 m or 5.5%. Figs. 9(b)–(d) and 10(b)–(d) concern the corresponding timehistories of the inline component of the horizontal velocity at z = (t), z = 0 and z = − 12 m. Once again, comparisons between the linear and nonlinear predictions confirm that there is a large nonlinear increase in the maximum horizontal velocity occurring immediately beneath the largest wave crest. For example, in the long-crested case (Fig. 9) the nonlinear solution gives umax = 10.5 m/s on z = (t) which corresponds to a 31.3% increase over the linearly predicted value. Similar comparisons in the shortcrested wave case (Fig. 10) indicate that umax = 7.75 m/s on z = (t) corresponding to a 9% increase relative to linear theory. The predicted horizontal velocities are significant in two respects. First, they indicate that the nonlinear increase in the near-surface waveinduced velocities can be very large, despite a relatively modest increase in crest elevation. Second, they confirm that the extent of this nonlinear increase is critically dependent upon the directionality of the sea state. These findings are consistent with earlier laboratory observations; uni-directional wave data being presented in [44–47] and multi-directional wave data in [39]. Figs. 9(e)–(g) and 10(e)–(g) provide a similar set of comparisons relating to the in line component of the unsteady water particle acceleration. Having already noted that the introduction of directionality reduces the in line wave-front steepness, the difference between the linear and nonlinear water particle accelerations might be expected to be somewhat smaller. Although this is indeed the case, significant differences in the magnitude of the unsteady acceleration remain. Evidence of this is provided in Figs. 9(e) and 10(e) in which the nonlinear increase in the maximum acceleration varies from 68% to 19%; the former relating to  = 10◦ and the latter to  = 30◦ .

6. Concluding remarks The present paper has considered the description of extreme ocean waves and has demonstrated how advances in wave modelling can be applied to provide a fully nonlinear representation of large transient waves arising in realistic ocean spectra. The physical characteristics that govern the evolution of the largest waves in an ocean environment are three-fold: (i) unsteadiness, due to the underlying frequency spectrum; (ii) nonlinearity, due to the steepness of the wave profile; and (iii) short-crestedness, due to the directional spread. Typical design procedures incorporate (i) or (ii), but never both, and seldom address (iii), except through the possible application of a crude velocity reduction factor applied

221

uniformly to all frequency components. In contrast, the present procedure allows (i), (ii) and (iii) to be fully included and provides full-field (both spatial and temporal) descriptions of the evolving water surface, (x, y, t), together with any aspect of the associated water particle kinematics, u(x, y, z, t). The proposed method builds upon developments in design wave philosophy [48] and could be used in conjunction with a number of existing fully nonlinear three-dimensional wave models, notably [1,20,24,25]. Given a valid set of initial conditions, based on the most probable shape of a large linear wave [9–11], the calculation procedures employed in the model are essentially exact and have been used to provide sample calculations for three cases involving JONSWAP spectra in relatively deep water. The first example concerns a uni-directional wave, whilst the other two include realistic directional spreads. In each case the difference between the linear and nonlinear predictions of the wave profile and the water particle kinematics are significant from a design perspective; the former being relevant to air-gap predictions and the latter load calculations. Two key elements of the proposed model concern its versatility and its computational efficiency, the latter being essential to maintain the former. Although the sample calculations relate to relatively deep water, the model is in no sense restricted to this flow regime. It can be successfully applied in a very wide range of water depths and is capable of incorporating any form of frequency and directional spectra. For example, in oceanographic circles there is an increasing trend towards the separation of wind sea and swell components and the application of frequency dependent directional spreading [49]. There is no reason why such descriptions should not easily be incorporated within the proposed model. Indeed, the only restrictions placed on the model arise from the assumed Fourier series representation. This requires both  and ()z= to be single-valued functions and periodic over some large fundamental length-scales, (x) and (y) . The former condition prohibits the modelling of over-turning waves, while the latter precludes the description of waves propagating over a variable bottom topography which is non-periodic. Although these restrictions are less important in deep water, they may become significant in shallow water coastal engineering applications. To overcome this difficulty, research is presently considering alternative calculation procedures in which  and ()z= can be both multi-valued and non-periodic. Finally, it is important to note that the comparisons contained within this paper contrast linear and nonlinear predictions for waves having an identical wave height. In many respects this is consistent with traditional design procedures in which the identification and description of an appropriate design wave is undertaken in two distinct stages. The first involves the statistical analysis of field data to determine the properties (usually Hmax and Tz ) of a design wave having a given return-period, perhaps 1 in 100 or 1 in 10,000 years. The second involves the application of a wave model to determine the corresponding crest elevation and the associated water particle kinematics. Clearly, the calculations presented in Section 5 address this second stage. Although this two-stage procedure is well established, the availability of a fully nonlinear wave model capable of predicting the evolution of extreme waves in realistic wavefields allows an alternative interpretation. In the first stage the statistical extrapolation involves field measurements of individual wave heights, the vast majority of which will be linear or near-linear. This is particularly true where realistic directional spreads are concerned. The extrapolation of this data to predict very rare events, up to and including the 1 in 10,000 year event, may not therefore adequately reflect the impact of nonlinearity on both Hmax and Tz . Indeed, it could be argued that the Hmax and Tz values so derived are only representative of the underlying linear sea state. If this is the case they

222

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

remain appropriate to the calculation of the initial conditions, necessary to start the numerical model, but no restrictions would be imposed on the final wave height resulting from the evolution of the wave field or the focusing of the wave components. If calculations of this form are undertaken, previous experimental observation (admittedly involving simplified laboratory spectra, [39]) suggest that Hmax would increase and Tz would reduce. These nonlinear effects would act in addition to those noted in Section 5 and might be expected to produce further large differences in the predictions of the maximum crest elevation and the underlying water particle kinematics. The statistical aspects of this problem have, in part, been addressed in [7]. Further consideration of the deterministic aspects, particularly relating to kinematics predictions, will be addressed in a follow-up part II paper.

difficulty in evaluating z . However, it has the disadvantage that whilst  is expressed in terms of a Fourier series,  is not. Accordingly, whilst derivatives of  can be rapidly calculated using fast Fourier transforms, a similar approach cannot be applied in the case of . As a result, the evaluation of the unknown coefficients requires large matrix inversion. This process is very costly in terms of computational effort and therefore restricts the number of wave components and hence the resolution that can be achieved in both the wavenumber and directional domains. To overcome this difficulty, the velocity potential given in Eq. (A.5) is evaluated at the water surface, z = , to yield

Appendix A.

Similarly, z , again evaluated on the water surface is given by

A.1. An efficient wave evolution model

(z )z= =

∞ ∞  

()z= =

akl cosh(K( + h)) ei(kx+ly) .

(A.6)

k=−∞l=−∞

∞ ∞  

akl K sinh (K( + h)) ei(kx+ly) .

(A.7)

k=−∞l=−∞

Within a three-dimensional fluid domain the wave motion is assumed to be irrotational and the velocity potential (x, y, z, t) is defined so that the velocity vector u = (u, v, w) = ∇ , where ∇ = (∂x , ∂y , ∂z ) and (x, y) defines a horizontal plane at the mean water level or z = 0. With the fluid assumed to be incompressible, mass continuity is expressed by Laplace’s equation, ∇ 2  = 0. This applies throughout the fluid domain, which is assumed to be bounded by a horizontal bed at z = − h and the free surface at z = (x, y, t). If the horizontal bed is impermeable, it follows that z = 0 on z = −h.

(A.1)

At the free surface (z = ) the kinematic and dynamic boundary conditions given in Eqs. (1) and (2) must be expanded to include the additional horizontal dimension, y. After some rearrangement these yield t = z − x x − y y , t = −g −

(A.2)

1 2 [] . 2

(A.3)

For reasons of computational efficiency (see below), both  and z= should be represented by Fourier series in the horizontal plane. Assuming that the wave motion is periodic over some large fundamental wavelengths, (x) and (y) in the x and y directions respectively, it follows that (x + (x) , y + (y) ) = (x, y) and (x + (x) , y + (y) ) = (x, y). Within this rectangular domain one class of solution is given by (x, y, t) =

+∞ +∞  

Akl ei(kx+ly) ,

(A.4)

with the corresponding velocity potential, satisfying both Laplace’s equation and the bottom boundary condition Eq. (A.1), given by ∞ ∞  

akl cosh(K(z + h)) ei(kx+ly) ,

k=−∞l=−∞

(z )z= = G()()z= .

(A.8)

Comparing Eqs. (A.6) – (A.8) the transformation implied by the Goperator simply involves the “multiplication” by K tanh (K( + h)) for all values of k, l, x and y. Unfortunately, the evaluation of tanh (K( + h)) is difficult because it contains information in both the wavenumber domain and the physical domain. To avoid this the cosh and sinh terms in Eqs. (A.6) and (A.7) are expanded using a Taylor series expansion about  = 0. Applying these expansions, “multiplications” appropriate to the two domains can be resolved, with Fourier transforms used to convert between the wavenumber and physical domains. Furthermore, substituting these expansions into Eq. (A.8) allows the G-operator to be evaluated at various orders. When attempting to express the solution in terms of surface parameter alone, it is important to note that the time derivative of the velocity potential on the surface is given by (z= )t = (t + z t )z= , where the second term on the right reflects the vertical motion of the water surface profile. Adopting this result and applying the transformation given in Eq. (A.8), allows the partial derivatives of  evaluated at z =  to be defined as (x )z= = (z= )x − x G()()z= ,

(A.9)

(y )z= = (z= )y − y G()()z= ,

(A.10)

(t )z= = (z= )t − t G()()z= .

(A.11)

Using these results the nonlinear free-surface boundary conditions given in Eqs. (A.2) and (A.3) can be re-written as

k=−∞l=−∞

(x, y, z, t) =

If the three-dimensional G-operator transforms ()z= into (z )z= , it follows that

(A.5)



where Akl and akl are functions of time only and K = Within a Fourier system k and l are integer multiples of the fundamental wavenumber components in the x and y directions respectively, so that k = nk0 and l = ml0 where k0 = 2/(x) , l0 = 2/(y) and n and m are integers. At this stage one could progress the solution by substituting Eqs. (A.4) and (A.5) directly into Eqs. (A.2) and (A.3) and calculate the time evolution of the unknown coefficients (Akl and akl ). This approach is akin to the method adopted in [26] and has the advantage that with  explicitly defined in terms of z there is no (k2

+ l2 ).

t = G()()z= − (x x − y y )z= , (z= )t = t G()()z= − −

(A.12)

1 2 (G()()z= ) 2

1 2 2 (x ) + (y ) − g. z= 2

(A.13)

In adopting this approach the required dimensional reduction is achieved and the solution formulated so that the evaluation of the surface derivatives given in Eqs. (A.12) and (A.13), a task typically undertaken twice per time-step, can be rapidly achieved using fast Fourier transforms. This latter point is fundamentally important since it ensures that the computational effort increases as O(N ln N), where N is the number of surface points or twice the number of wave components. This is in stark contrast to other directional

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

wave models [26], where the necessary inclusion of large matrix inversion means that the computational effort increases as O(N3 ), with consequent restrictions on the maximum number of wave components that can realistically be included. A.2. Calculation of the water particle kinematics From an engineering perspective methods to predict the evolution of the water surface elevation represent only part of the required solution: a practical wave model must also provide a corresponding description of the water particle kinematics. Regarding the model outlined in [1], and briefly reviewed above, the calculation of the underlying kinematics is complicated by the necessary dimensional reduction. With the velocity potential only defined on the water surface, (x, y, , t), it is not immediately clear how this should be used to calculate the kinematics within the interior of the flow −h ≤ z ≤ . Although there are a number of possible procedures, those that involve large matrix inversion must be avoided (for the reasons noted above), as must any solution whose accuracy or reliability becomes uncertain close to the water surface. This latter point is particularly important when seeking to define the fluid motion associated with large highly nonlinear waves. For example, if one adopts a Fenton and Rienecker [21] type solution for the velocity potential the amplification of small numerical errors, due to rounding and/or alaising, can lead to very misleading kinematic predictions close to the water surface. To overcome this difficulty, [2] proposed a two-stage transformation. In the first stage any arbitrary kinematic quantity calculated on the water surface, (x, y, ), is transformed onto progressively flatter surfaces by the repeated application of the so-called H2 -operator H2 (, 2 )(x, y, ) = (x, y, 2 ),

(A.14)

where 2 defines the flatter surface onto which the kinematic quantity is transformed. Having repeated the process until  is transformed onto a sufficiently flat surface (full details of which are given in [2]), the so-called H-operator is applied to define the global parameters, ˛kl , that describe the kinematic quantity everywhere within the fluid domain. The second stage is defined by H() s =

∞ ∞  

˛kl sinh(Kh)ei(kx+ly) ,

(A.15)

k=−∞l=−∞

where the superscript s denotes the value of the quantity evaluated on the surface n resulting from (n − 1) applications of Eq. (A.14) such that  s (x, y, t) = (x, y, n (x, y, t), ) and all other symbols are as defined previously. Contrasting Eqs. (A.14), (A.15) and (A.8) it is clear that both the H- and H2 -operators can be written in a form similar to the G-operator. Indeed, although the H- and H2 -operators undertake different tasks, their structure and formulation are closely related to the three-dimensional G-operator proposed in [1]. It is also clear that both the H- and H2 -operators could be used individually to calculate the value of the kinematic quantity at any or all elevations. However, if the H-operator is applied directly to the kinematic quantity on the water surface (i.e. there is no prior application of the H2 -operator such that n =  in Eq. (A.15)), then significant errors may arise at the water surface. This occurs because of the amplification of small numerical errors that become progressively larger as the nonlinearity of the wave form increases. In many respects this effect is akin to the process of high-frequency contamination discussed in [14]. Alternatively, if the H2 -operator is applied alone, an excessive number of repeated calculations would be required to determine the value of the kinematic quantity at appreciable depths beneath the water surface. Such calculations would again

223

involve the accumulation of small numerical errors. In contrast, if the H- and H2 - operators are used in tandem (as is proposed in [2]), the method provides a stable, accurate and highly efficient procedure to determine the value of any kinematic quantity based on knowledge of (x, y, t) and (x, y, , t). When this approach is taken together with the wave evolution model outlined in Appendix A.1, the procedures represent an extremely efficient, fully nonlinear, directional wave evolution and wave kinematics model. References [1] Bateman WJD, Swan C, Taylor PH. On the efficient numerical simulation of directionally-spread surface water waves. J Comput Phys 2001;174:277–305. [2] Bateman WJD, Swan C, Taylor PH. On the calculation of the water prticle kinematics arising in a directionally spread wavefield. J Comput Phys 2003;186:70–92. [3] ISSC. Report of the Environment Committee. In: 14th International Ship and Offshore Structures Congress. 2000. [4] Walker DAG, Taylor PH, Eatock-Taylor R. The shape of large surface waves on the open sea and the draupner new year wave. Appl Ocean Res 2004;26: 73–83. [5] Gibbs RH, Taylor PH. Formation of walls of water in ‘fully’ nonlinear simulations. Appl Ocean Res 2005;27:142–57. [6] Gibson RS, Swan C. The evolution of large ocean waves: The role of local and rapid spectral changes. Proc R Soc Lond A 2007;463:21–48. [7] Gibson RS, Swan C, Tromans PS. Fully nonlinear statistics of wave crest elevation calculated using a spectral response surface method: applications to unidirectional sea states. J Phys Oceanogr 2007;37(1):3–15. [8] Katsardi V, Swan C. The evolution of large non-breaking waves in intermediate and shallow water. Part I. Numerical calculations of uni-directional seas. Proc Roy Soc Lond A 2011;467:778–805. [9] Lindgren G. Some properties of a normal process near a local maximum. Ann Math Stat 1970;41:1833–70. [10] Bocotti P. Some new results on statistical properties of wind waves. Appl Ocean Res 1983;5(3):134–40. [11] Tromans PS, Anaturk AR, Hagemeijer P. A new model for the kinematics of large ocean waves. In: Proc. First Intern. Confer. on Offshore and Polar Engineering, vol. 3. 1991. p. 1–64. [12] Dean RG. Stream function representation of non-linear ocean waves. J Geophys Res 1965;70:4561–72. [13] Wheeler J. Method for calculating forces produced by irregular waves. In: Proc Offshore Tech Conf, vol. 1. 1970. p. 71–82. [14] Sobey RJ. Wave theory predictions of crest kinematics. In: Water wave kinematics. Kluwer; 1990. [15] Lo E, Mei C. A numerical study of water-wave modulation based on a higher order Schrödinger equation. J Fluid Mech 1985;150:395–416. [16] Dysthe KB. Note on the modification of the nonlinear Schrödinger equation for application to deep water waves. Proc R Soc Lond A 1979;369:105–14. [17] Zakharov VE. Stability of periodic waves of finite amplitude on the surface of a deep fluid. J Appl Mech Tech Phys 1968;9:190–4 [English translation]. [18] Longuet-Higgins MS, Cokelet ED. The deformation of steep surface water waves on water. A numerical method of computation. Proc Roy Soc Lond A 1976;350:1–26. [19] Dommermuth DG, Yue DKP. A high-order spectral method for the study of nonlinear gravity waves. J Fluid Mech 1987;184:267–88. [20] West BJ, Brueckner KA, Janda RS, Milder DM, Milton RL. A new numerical method for surface hydrodynamics. J Geophys Res 1987;92(C11): 11803–24. [21] Fenton JD, Rienecker MM. A Fourier method for solving nonlinear waterwave problems: application to solitary-wave interactions. J Fluid Mech 1982;118:411–43. [22] Dold JW, Peregrine DH. Steep unsteady water-waves: An efficient computational scheme. In: Proc. 19th Conf. Coastal Engng., vol. 1. Houston: ASCE; 1984. p. 955–67. [23] Craig W, Sulem C. Numerical simulation of gravity waves. J Comput Phys 1993;108:73–83. [24] Ducrozet G, Bonnefoy F, Le Touzé D, Ferrant P. 3-D HOS simulations of extreme waves in open seas. Nat Hazards Earth Syst Sci 2007;7:109–22. [25] Xu L, Guyenne P. Numerical simulation of three-dimensional nonlinear water waves. J Comput Phys 2009;228:8446–66. [26] Johannessen TB, Swan C. On the nonlinear dynamics of focused wave groups in two and three dimensions. Proc Roy Soc Lond, Ser A 2003;459:1021–52. [27] Fructus D, Clamond D, Grue J, Kristiansen Ø. An efficient model for threedimensional surface wave simulations. Part I. Free space problems. J Comput Phys 2005;205:665–85. [28] Clamond D, Fructus D, Grue J, Kristiansen Ø. An efficient model for threedimensional surface wave simulations. Part II. Generation and absorption. J Comput Phys 2005;205:686–705. [29] Toffoli A, Benoit M, Onorato M, Bitner-Gregersen E. The effect of third-order nonlinearity on statistical properties of random directional waves in finite depth. Nonlinear Process Geophys 2009;16:131–9. [30] Grue J. Computation formulas by FFT of the nonlinear orbital velocity in threedimensional surface wave fields. J Eng Math 2010;67:55–69.

224

W.J.D. Bateman et al. / Applied Ocean Research 34 (2012) 209–224

[31] Schäffer HA. Comparison of Dirichlet–Neumann operator expansions for nonlinear surface gravity waves. Coast Eng 2008;55:288–94. [32] Hasselmann K. On the nonlinear energy transfer in a gravity wave system. Part 1. J Fluid Mech 1962;12:481–500. [33] Hasselmann K. On the nonlinear energy transfer in a gravity wave system. Part 2. J Fluid Mech 1963;15:273–81. [34] Phillips OM. The dynamics of the upper ocean. 2nd ed Cambridge University Press; 1977. [35] Rozario JB, Tromans PS, Taylor PH, Efthymiou M. Comparisons of loads predicted using NewWave and other models with measurements on the Tern structure. Adv Underwater Tech Ocean Sci Offshore Eng 1993;29: 143–58. [36] Jonathan P, Taylor PH. Irregular nonlinear waves in a spread sea. ASME Trans J Offshore Mech Arctic Eng 1996;119:37–41. [37] Baldock TE, Swan C. Numerical calculations of large transient water waves. Appl Ocean Res 1994;16:101–12. [38] Hasselmann DE, Dunckel M, Ewing J. Directional wave spectra observed during JONSWAP 1973. J Phys Oceanogr 1980;10:1264–80. [39] Johannessen TB, Swan C. A laboratory study of the focusing of transient and directionally spread surface water waves. Proc Roy Soc Lond, Ser A 2001;457:1–36.

[40] Mitsuyasu H. Observations of directional spectrum of ocean waves using a cloverleaf buoy. J Phys Oceanogr 1975;16:750–60. [41] Longuet-Higgins MS, Stewart RW. Changes in the form of short gravity waves on long waves and tidal currents. J Fluid Mech 1960;8:565–83. [42] Longuet-Higgins MS. The effect of non-linearities on statistical distributions in the theory of sea waves. J Fluid Mech 1963;17:459–80. [43] Sharma JN, Dean RG. Second-order directional seas and associated wave forces. Soc Petrol Eng J 1981;4:129–40. [44] Kim CH, Randall E, Boo SY, Krafft MJ. Kinematics of 2-d transient water waves using laser doppler anemometry. J Waterway, Port, Coastal Ocean Eng 1992;118:147–65. [45] Baldock TE, Swan C. A laboratory study of non-linear surface waves on water. Philos Trans R Soc Lond, Ser A 1996;354:649–76. [46] Grue J, Clamond D, Huseby M, Jensen A. Kinematics of extreme waves in deep water. Appl Ocean Res 2003;25:355–66. [47] Grue J, Jensen A. Experimental velocities and accelerations in very steep wave events in deep water. Eur J Mech B/Fluids 2006;25:554–64. [48] Smith D, Birkenshaw M. Iso (draft) offshore design codes: Wave loading requirements. Trans Inst Marine Eng 1996;108:89–96. [49] Ewans K. Directional spreading in ocean swell. In: Proc. WAVES 2001, vol. 1. San Francisco, USA: ASCE; 2001. p. 517–29.