TIME-DOMAIN SOLUTION FOR SECOND-ORDER W A V E DIFFRACTION by Kwok Fai Cheung B . A . S c , University of Ottawa, 1985 M . A . S c , University of British Columbia, 1987 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH C O L U M B I A April 1991 ©KwokFa i Cheung, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of C i v i l Engineering The University of British Columbia Vancouver, Canada Date A p r i l 12. 1991 DE-6 (2/88) Abstract A numerical method based on potential flow theory is developed for simulating transient, second-order interactions of ocean waves with large fixed bodies of arbitrary shape in two and three dimensions. The physical problem is represented by a mathematical model composed of a fluid domain bounded by the body surface, the still water surface, the seabed, and a control surface truncating the infinite fluid region. The nonlinear free surface boundary conditions defined on the instantaneous free surface are expanded about the still water level by a Stokes expansion procedure. The flow potential to second order is thereby defined with respect to a time-independent boundary which includes the still water surface, and its solution involves a numerical discretization of an integral equation. With the potential separated into incident and scattered components, the Sommerfeld radiation condition applied to the scattered potential is modified to incorporate a time-dependent celerity to account for the transient and second-order effects. The free surface boundary conditions and the radiation condition are then satisfied to second order by a numerical integration in time. An alternative second-order solution is derived based on a different expansion procedure in which the nonlinear free surface boundary conditions and an integral equation defined on the instantaneous free surface are both expanded by a Taylor series, and terms up to second order are retained. The two approaches give rise to identical first-order problems, but give rise to second-order problems which are apparendy different. The discrepancy arises from the second-order integral equation in which additional second-order terms are retained. The physical interpretations and limitations of these terms are explored and their effects on the evaluations of wave forces are assessed. Applications of the present method are made to studies of regular wave diffraction around a fully submerged and a semi-submerged circular cylinder in two dimensions, and around a bottom-mounted surface-piercing circular cylinder in three dimensions. The stability and ii numerical accuracy of the proposed solution and the treatment of the radiation condition to second order are examined. Comparisons of computed wave forces and runup are made with previous theoretical and experimental results and these indicate favourable agreement. iii Table of Contents Abstract ii List of Tables vii List of Figures viii List of Symbols xi Acknowledgement xv 1. INTRODUCTION 1 1.1 General 1 1.2 Literature Review 2 1.2.1 Perturbation Methods 3 1.2.2 Time-Stepping Methods 5 1.3 Objectives of Present Investigation 7 2. THEORETICAL FORMULATION 8 2.1 Full Nonlinear Method 8 2.2 Stokes Second-Order Method 11 2.2.1 Perturbation Expansion 11 2.2.2 Three-Dimensional Problem 14 2.2.3 Two-Dimensional Problem 21 2.2.4 Time-Integration and Initial Conditions 22 2.3 Alternative Second-Order Method 24 2.3.1 Second-Order Expansion of Integral Equation 24 2.3.2 Alternative Second-Order Integral Equation 26 iv 2.4 Wave Runup and Forces 27 2.4.1 Free Surface Elevations and Runup 27 2.4.2 Wave Forces and Moments 29 3. N U M E R I C A L F O R M U L A T I O N 33 3.1 Field Solution 33 3.1.1 Stokes Second-Order Solution 33 3.1.2 Alternative Second-Order Solution 40 3.2 Temporal Solution 43 3.2.1 Radiation Condition 43 3.2.2 Free Surface Boundary Conditions • • • • 47 4 . R E S U L T S A N D D I S C U S S I O N 49 4.1 Introduction 49 4.2 Computational Considerations 50 4.2.1 Discretization Schemes 50 4.2.2 Computing Cost. 52 4.2.3 Accuracy and Stability of Matrix Solution 53 4.2.4 Initial Conditions 56 4.2.5 Stability of Time-Stepping Procedure 58 4.2.6 Numerical Accuracy of Discretization 60 4.2.7 Effectiveness of Radiation Condition 63 4.2.8 Irregularity of Second-Order Solution 67 4.3 Submerged Circular Cylinder in Two Dimensions 68 4.3.1 Free Surface Profiles 68 4.3.2 Wave Forces 69 4.3.3 Applicability of Two Methods 72 v 4.4 Semi-Submerged Circular Cylinder in Two Dimensions 73 4.4.1 Free Surface Profiles and Runup 73 4.4.2 Wave Forces 74 4.5 Surface-Piercing Circular Cylinder in Three Dimensions 76 4.5.1 Free Surface Profiles 76 4.5.2 Wave Runup 78 4.5.3 Wave Forces and Moments 80 5. CONCLUSIONS 82 5.1 Theoretical and Numerical Procedures 82 5.2 Numerical Studies 84 5.3 Recommendations for Further Study 86 References 88 vi List of Tables 1. Minimum T/At as a function of L/AS for stable solutions of regular wave diffraction at first and second order around a semi-submerged circular cylinder for ka = 1 and deep water 93 2. Amplitudes of wave force components at first and second order on a submerged circular cylinder as functions of Nb for ka = 0.4 and deep water 93 3. Amplitudes of free surface elevations at first and second order on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water 94 4. Amplitudes of wave force components at first and second order in the x direction on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water. 95 5. Amplitudes of wave force components at first and second order in the z direction on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water . 95 6. Amplitudes of free surface elevations at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1 96 7. Amplitudes of wave force components at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1 96 8. Amplitudes of wave moment components at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1 . 96 vii List of Figures 1. Definition sketch 97 2. Examples of surface discretization in two dimensions, (a) fully submerged circular cylinder with Nb = 30, N c = 30, N 0 = 76 and N = 136. (b) semi-submerged circular cylinder with N b = 15, N c = 30, N 0 = 64 and N = 109 97 3. Examples of surface discretization of a bottom-mounted surface-piercing circular cylinder in three dimensions, (a) Nb = 144, N c = 312, N Q = 1358 and N = 1824. (b) N b = 220, N c = 312, N 0 = 1402 and N = 1934 98 4. C P U time as a function of the rank of the matrix equation for the two-dimensional problem 99 5. C P U time as a function of the rank of the matrix equation for the three-dimensional problem 99 6. Condition number as a function of the rank of the matrix equation for the two-dimensional problem 100 7. Condition number as a function of the rank of the matrix equation for the three-dimensional problem . 100 8. Free surface profiles to first and second order around a semi-submerged circular cylinder obtained by different modulation times T m for ka = 0.6, H /L = 0.1, t/T = 5 and deep water 101 9. Celerity and flow potential of a bichromatic wave train at x = 0 as functions of time for ai/a2 = 3, C1/C2 = 0.5 and deep water, (a) celerity, (b) flow potential 102 10. Celerity and flow potential of a bichromatic wave train at x = 0 as functions of time for ai/a2 = 2, C1/C2 = 0.8 and deep water, (a) celerity, (b) flow potential 103 11. Development with time of oscillatory force components at first and second order on a semi-submerged circular cylinder subject to various treatments at radiation boundaries for ka = 0.6 and deep water, (a) components due to ty^. (b) components dueto<J>2 104 12. Development with time of free surface profiles to first and second order around a submerged circular cylinder for ka = 0.5, h/a = 2, H /L = 0.1 and deep water 105 13. Incident and developed wave profiles to first and second order around a submerged circular cylinder at t/T =8 for ka = 0.5, h/a = 2 and H/L = 0.1 106 14. Development with time of wave force components to first and second order on a submerged circular cylinder for ka = 0.5, h/a = 2, H /L = 0.1 and deep water . . . . 106 15. Amplitude of first-order oscillatory force and second-order steady drift force on a submerged circular cylinder as functions of ka in deep water, (a) first-order force amplitude, (b) second-order vertical drift force 107 viii 16. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of ka in deep water 108 17. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of ka in deep water 109 18. Amplitudes of second-order oscillatory force components on a submerged circular cylinder as functions of ka for kd = 2. (a) x component, (b) z component 110 19. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of h/a for ka = 0.35 and deep water I l l 20. Amplitudes of second-order oscillatory force components on a submerged circular cylinder as functions of h/a for ka = 0.35 and kd = 2 I l l 21. Development with time of free surface profiles to first and second order around a semi-submerged circular circular for ka = 0.6, H/L = 0.1 and deep water 112 22. Development with time of free surface elevations to first and second order adjacent to a semi-submerged circular circular for ka = 0.6, H / L = 0.1 and deep water. (a) upwave surface, (b) down wave surface 113 23. Runup to second order and amplitudes of free surface elevations at first and second order on a semi-submerged circular cylinder as functions of ka for H / L = 0.1 and deep water, (a) upwave surface, (b) down wave surface 114 24. Development with time of wave force components to first and second order on a semi-submerged circular cylinder for ka = 0.6, H /L = 0.1 and deep water 115 25. Amplitudes of first-order oscillatory force components on a semi-submerged circular cylinder as functions of ka in deep water, (a) x component, (b) z component 116 26. Amplitudes of second-order oscillatory force components on a semi-submerged circular cylinder as functions of ka in deep water, (a) x component, (b) z component 117 27. Second-order steady force components on a semi-submerged circular cylinder as functions of ka in deep water 118 28. Amplitudes of second-order oscillatory force components on a semi-submerged circular cylinder as functions of ka in deep water, (a) x component, (b) z component 119 29. Views of free surface to first order around a surface-piercing circular cylinder for ka = 1 and kd = 1. (a) 9/2TC = 0.85. (b) 0/27C = 0.35 120 30. Views of free surface to second order around a surface-piercing circular cylinder for ka = 1, kd = 1 and H/L = 0.08. (a) 8/271 = 0.85. (b) 9/2TC = 0.35 121 ix 31. Runup and rundown to second order and amplitudes of free surface elevations at first and second order on a surface-piercing circular cylinder as functions of azimuthal angle, (a) ka = 0.481, kd = 1.332 and H / L = 0.0697. (b) ka = 0.684, kd = 1.894 and H/L = 0.091. (c) ka = 0.917, kd = 2.536 and H/L = 0.1004 122 32. Runup to second order and amplitudes of free surface elevations at first and second order at 9 = % on a surface-piercing circular cylinder as functions of ka for d/a = 1 and H/L = 0.1 124 33. Maximum amplitude (over 9) of free surface elevation at second order and location of maximum on a surface-piercing circular cylinder as functions of ka for a/d = 1. (a) free surface elevation, (b) location of maximum 125 34. Development with time of wave force and moment to first and second order on a surface-piercing circular cylinder for ka = 1, kd = 1 and H/L = 0.1 126 35. Wave force and moment to second order and components at first and second order on a surface-piercing circular cylinder as functions of ka for d/a = 1 and H/L = 0.1. (a) wave force components, (b) moment components 127 36. Amplitudes of second-order oscillatory force and moment components on a surface-piercing circular cylinder as functions of ka for d/a = 1 . (a) wave force components, (b) moment components 128 x List of Symbols The following symbols are used in this thesis: a = radius of circular cylinder; A = incident wave amplitude; aj, bj = matrix coefficients to calculate potential at an interior point; Ay , Bjj = matrix coefficients for the simultaneous equations; c n = celerity in the normal direction; Cy, Dij , Ejj, Fy = additional matrix coefficients for the alternative second-order method; d = water depth (see figure 1); D = domain bounded by Sb, S c and SQ; 3 = domain bounded by S w , S c and Sf, F = force; Fx. F-x = force maxima to second order in the positive and negative x directions; F m = modulation function defined in (4.1); g = gravitational constant; G(x,t;) = Green's function; h = depth of submergence measured from the cylinder axis to the still water level in the two-dimensional problem; H = incident wave height; J = Jacobian; k = wave number, L = wavelength; LQ = length of domain in the two-dimensional problem; M = moment; and number of Gauss points in one direction; My, M . y = moment maxima to second order about the y axis in the clockwise and counter-clockwise directions; xi n = normal vector given by (n x , n 2) in the two-dimensional problem and by (nX ) n y , n z) in the three-dimensional problem; n = distance in the direction of normal; n' = defined in (2.60); n x , n y , n z = direction cosines of the normal vector with respect to the x, y and z directions; Nj = shape function defined in (3.8); N = total number of facets; and rank of the matrix equation; N ' = number of facets on Sb and S c; N b , N c , N Q = numbers of facets on S b , S c and S 0 ; p = hydrodynamic pressure; P a = atmospheric pressure; r = radial distance measured from the body; and distance between x and!;; R = runup to second order; R ' = rundown to second order, (s, t) = local coordinates at a facet; S = total surface composed of S w , S c and Sf, S' = total surface composed of Sb, S c and S 0 ; Sb = body surface below still water level; S c = control surface (see figure 1); Sf = instantaneous free surface (see figure 1); S 0 = still water surface (see figure 1); S w = instantaneous wetted body surface (see figure 1); t = time; T = wave period; T m = modulation time; u = fluid particle velocity; xii U = Ursell parameter, w - distance measured along the water-line contour, w 0 = intersections between the still water surface and body surface in the two-dimensional problem; and still water-line contour in the three-dimensional problem; w m , w n = weight functions for the Gaussian integration; x = position vector given by (x, z) in the two dimensional problem and by (x, y, z) in the three-dimensional problem; (x, y, z) = Cartesian coordinate system (see figure 1); (3 = Courant number, 8 = dirac-delta function; An = distance of the order of a facet diameter in the direction of the normal; AS = facet diameter in the two-dimensional problem; and facet area in the three-dimensional problem; At = time-step size; e = perturbation parameter, <(> = velocity potential; y = angle between the radial direction and the normal on the control surface; r = integral given by (3.11) in the two-dimensional problem and by (3.5) in the three-dimensional problem; T| = free surface elevation (see figure 1); 9 = azimuthal angle measured horizontally from the x axis in the three-dimensional problem; and phase angle of incident waves; p = fluid density; co = angular frequency; \ , \|/, £ = distances in the x, y and z directions; and xiii £ = position vector given by (£, £) in the two-dimensional problem and by Y> 0 in the three-dimensional problem. The following subscripts are used in this thesis: 0 = steady component at second order, 1 = component at first order, 2 = component at second order; i, j = indices for field and source points; k = index for image point; and order of the problem; m, n = indices for Gauss points; and x, y, z = components in the x, y and z directions. The following superscripts are used in this thesis: (1) = component due to potential at first order, (2) = component due to potential at second order, s = scattered component; and w = incident component xiv Acknowledgement The author would like to express his sincere appreciation and gratitude to his thesis adviser, Prof. M . Isaacson, for his advice and guidance throughout this research. Thanks are also extended to Mr. T. Nicol for his advice in developing and running computer codes on the I B M 3090/150S computer at the University of British Columbia. The assistance of Mr. S. Prasad was instrumental in animating numerical results on a Macintosh IIx computer. The author also wants to thank his wife Wendy for her support and encouragement. Financial support from the University of British Columbia in the form of a University Graduate Fellowship and from the Natural Sciences and Engineering Research Council of Canada in the form of a research assistantship is gratefully acknowledged. xv CHAPTER 1 INTRODUCTION 1.1 General Numerical modelling of nonlinear wave diffraction around large offshore structures has been the subject of investigation for a number of years. The original motivation of such studies arises primarily because of the need to obtain more accurate wave force and runup predictions than those of linear cuffraction theory which are based on the assumption of infinitesimal wave heights. With the advancement of theories on nonlinear wave diffraction, an increasing number of practical problems inherent in the design of ships and large offshore structures can now be investigated theoretically. Applications of practical significance which cannot be accounted for by linear diffraction theory include the low and high frequency excitations and responses of offshore structures in ocean waves. Although the magnitude of these nonlinear components is generally small, they act at frequencies which are significantly different from those of the incident wave spectrum, and may be critical when such excitations are near the natural frequencies of the body motions or when the damping forces are small. Typical examples are low frequency drift motions of moored offshore structures or tankers (e.g. Pinkster & Wichers, 1987; and Matsui, 1989) and high frequency excitation of tension-leg platforms (e.g. K i m & Yue, 1988). For some structural geometries such as a triangular ship section and a conical platform, the second-order component has been found to be dominant in the total wave force by theoretical and experimental studies (e.g. Kyozuka, 1981b; Kyozuka, 1982; and Jamieson et al., 1985). For surface-piercing bodies in waves, it is well known that nonlinear effects account for a significant portion of the total wave runup. In some regions, the design criteria for wave runup and overtopping can be significant in determining the design freeboard height of an offshore 1 structure (e.g. Wang et al. 1983). For a vertical circular cylinder in regular waves, Kriebel (1987) (also in Kriebel, 1990) has shown that in some cases both the measured runup and second-order solution exceed the linear prediction by 50 percent When the body sidewalls are outward sloping towards the bottom, such as in the case of a vertical cone, the runup component at second order was found to be greater than the first-order runup for moderately steep incident waves (e.g. K im & Yue, 1989). Despite their well-known practical importance, theoretical and numerical studies of nonlinear wave diffraction are mainly confined within the academic domain. Existing nonlinear models are either algebraically so complicated or computationally so intensive that general applications are somewhat restricted. In the present study, a relatively simple, efficient and versatile numerical method is developed to provide solutions for a variety of practical problems in the industry. 1.2 Literature Review Wave force and runup predictions for large offshore structures are generally made on the basis of linear diffraction theory, which is formally valid for small amplitude sinusoidal waves. (For a review of the methods used see Sarpkaya & Isaacson, 1981; and Chakrabarti, 1987.) Due to the demand of a more accurate design procedure, recent research has been directed towards the development of nonlinear models to account more realistically for the effects of large wave heights. During the last two decades, a considerable number of papers and research reports have been published on the subject. These methods involve the solution of a potential flow problem and can generally be classified into two categories. One is a second-order solution obtained by a perturbation procedure in the frequency domain, while the other is a full nonlinear solution obtained by a time-stepping procedure. 2 1.2.1 Perturbation Methods In the perturbation methods, linear diffraction theory for regular waves is extended to a second approximation on the basis of the Stokes expansion procedure in which the variables describing the flow are developed as power series in terms of a perturbation parameter, analogous to the development of Stokes second-order wave theory in the frequency domain. However, the method involves two inhomogeneous free surface boundary conditions and an inhomogeneous radiation condition for which an exact solution has been extremely difficult to obtain. Most of the pioneer work failed to satisfy the inhomogeneous free surface boundary conditions or to recognize the co-existence of second-order forced and free waves in the solution and in the radiation condition (e.g. Yamaguchi & Tsuchiya, 1974; and Raman & Venkatanarasaiah, 1976). Most subsequent studies seemingly attempted to correct these deficiencies. These methods generally involve the decomposition of the potential at second order into free and forced components satisfying respectively the homogeneous and inhomogeneous terms in the free surface boundary conditions, and the radiation condition at second order is also modified accordingly in order to account for the forced term. If interest is restricted to the wave loading on bodies, the wave forces at second order can be obtained without evaluating the potential at second order explicitly by making use of the Haskind reciprocal relationship. A n asymptotic form of the potential at second order which is sometimes known as the "weak" radiation condition is used to provide a far-field closure of the solution. This approach, applicable in principle to three-dimensional structures of arbitrary shape, was independently proposed by Lighthill (1979) for infinite water depth and for small piles; and by Molin (1979) for finite water depth. Refinements of the method were made to treat regular wave diffraction around a bottom-mounted surface-piercing circular cylinder in three dimensions (e.g. Rahman, 1984; Demirbilek & Gaston, 1985; and Eatock Taylor & Hung, 1987), and around a truncated circular cylinder in three dimensions (Abul-Azm & 3 Williams, 1988). Rahman (1987) and Ghalayini & Williams (1989) have applied the method to bodies of arbitrary shape in three dimensions, and Molin & Marion (1986) have extended the method to include the effects of body motions. The application of the method to an array of vertical cylinders in three dimensions has been made by Williams et al. (1990) and Masuda & Nagai (1991), while the application of the method to a submerged cylinder in two dimensions has been made by Wu & Eatock Taylor (1990). Although this approach has recently been extended to calculate local pressure at second order on a vertical surface-piercing circular cylinder by Eatock Taylor et al. (1989), and to account for the diffraction of subharmonic waves by Matsui (1989), the second-order wave kinematics, which may be just as important, can not be evaluated by these methods. On the other hand, the potential at second order may be obtained explicitly by satisfying the free surface boundary conditions and radiation condition at second order. The wave forces and kinematics at second order such as the free surface profile and runup can be obtained subsequently from the potential. The three-dimensional problem involving a vertical circular cylinder has been studied extensively by Hunt & Baddour (1981), Hunt & Williams (1982), Chen & Hudspeth (1982), Kriebel (1987) and Shimada (1987), while Garrison (1984) has developed a method to treat problems involving structures of arbitrary shape. Huang et al. (1989) have examined the use of an asymptotic form of the potential at second order in solving the second-order diffraction problems in two and three dimensions. With the use of the asymptotic solution as described by Molin (1979), Kim & Yue (1989) have solved the second-order diffraction solution for axisymmetric bodies and subsequently extended the method to include bichromatic wave diffraction and body motions (Kim & Yue, 1990). More recently, Newman (1990) has derived an approximation to the diffraction potential at second order for the case of large depth. The two-dimensional diffraction problem in infinite water depth has been treated by Kyozuka (1980), Vada (1987), Inoue & Kyozuka (1988), Miao & Liu (1989) and Mclver & Mclver (1990), while the two-dimensional radiation problem in infinite water depth has been treated by Yamashita (1977), Kyozuka (1981a) and Papanikolaou (1984). This 4 method has been extended to include the effects of body motions in waves by Kyozuka (1983), and to include bichromatic wave diffraction by Kyozuka (1988). Kriebel (1987) compared most of the second-order solutions prior to his work and observed that these solutions exhibit noticeable differences, and apparently depend on the way in which the free surface boundary conditions and the radiation condition are satisfied. Comparisons among most of these solutions are inadequate and are further obscured by the fact that most of these earlier results were presented and compared in terms of the total second-order forces rather than the force components at second order which themselves are one order smaller. On the basis of the asymptotic expansion of the potential at second order in three dimensions as described by Molin (1979), solutions at second order obtained by Molin & Marion (1986), Eatock Taylor & Hung (1987), Abul-Azm & Williams (1988) and K i m & Yue (1989) have been shown to agree with each other. In two dimensions, Wu & Eatock Taylor's (1990) analytic solution has also been shown to agree with Vada's (1987) numerical solution. 1.2.2 Time-Stepping Methods A second approach has been to solve the full nonlinear boundary value problem numerically by a time-stepping procedure, in which the free surface boundary conditions are applied on the instantaneous water surface and a new system of simultaneous equations are generated and solved at each time step as the free surface moves to a new position. Since the introduction of the method by Longuet-Higgins & Cokelet (1976), considerable progress has been made on the accurate calculation of two-dimensional flows, including breaking waves in the presence of a fixed obstacle (e.g. Vinje & Brevig, 1981; and Jagannathan, 1988), and breaking waves generated by a wave-maker (e.g. Lin et al., 1984; and Grosenbaugh & Yeung, 1989). Most formulations of this general method have been developed on the basis of a complex potential and are restricted to applications in two dimensions. Restricting the movement of the free surface in the vertical direction, the three-5 dimensional problem applicable to bodies of arbitrary shape in waves was treated by Isaacson (1982) on the basis of an integral equation method based on Green's theorem, and the corresponding two-dimensional problem with a fully submerged circular cylinder was treated by Stansby & Slaouti (1984) on the basis of a vortex sheet method. However, these methods require a substantial computational effort which limits their applications, and possible numerical instabilities or errors and the effects of waves reflected from the control surface provide further limitations to the length of simulation time. It is difficult to achieve a stable, steady state solution over many wave periods by these methods. The manner of treating the radiation condition with a domain of finite size requires particular attention. In the three-dimensional radiation problem, corresponding to waves generated by the motion of a body in otherwise still water, the amplitude of radiated waves as well as the effects of free surface nonlinearities decrease with the distance away from the body, so that the nonlinear solution in the neighbourhood of the body may be matched to a linear solution on the control surface (e.g. Lin et al., 1984; and Dommermuth & Yue, 1987). This procedure is untenable in two-dimensional flow since the nonlinearity of the radiated waves persists to the far-field and the matching boundary condition is problematic (e.g. Vinje et al., 1982). Furthermore, such an approach cannot be used to treat the diffraction problem, since the flow then remains markedly nonlinear throughout the computational domain. On the other hand, Orlanski (1976) has suggested the use of the Sommerfeld radiation condition to calculate the time-dependent celerity on the boundary of a hyperbolic flow which in turn may be used to estimate the boundary conditions on the control surface for the next time step. Applications of Orlanski's method to the radiation problem of full nonlinear waves were made by Brorsen & Larsen (1987) with a constant celerity calculated from the linear dispersion relation, and by Chan (1977), Yen & Hall (1981) and Jagannathan (1988) with a time-dependent celerity calculated near the control surface. Although this numerical radiation condition was originally developed for hyperbolic flows, its applications in parabolic problems have also been shown to be effective by Han et al. (1983). 6 1.3 Objectives of Present Investigation Most of the previous studies on second-order wave diffraction have been found to be algebraically complicated and are restricted for applications to regular or bichromatic wave diffraction in the frequency domain. Study such as diffraction of a deterministic irregular wave profile to second order which is of practical interest is somewhat complicated by the stochastic analysis in the bi-frequency domain. On the other hand, the full nonlinear diffraction in the time domain is computationally intensive in which a new system of simultaneous equations must be generated and solved at each time step. Numerical errors introduced into the procedure in solving the system of simultaneous equations and which accumulate through the time-stepping provide a further Umitation to the length of simulation time. In addition, the absence of a proper radiation condition for the full nonlinear diffraction problem makes it more difficult to achieve a stable, steady state solution over many wave periods. If the wave height is moderate and wave breaking is not considered, the second-order approach to the nonlinear diffraction problem, although less accurate in principle, appears to be more effective and reliable. Despite the recent advances in the theories of second-order wave diffraction, there is still a need to verify existing second-order solutions and to develop a more versatile approach for practical problems applicable to industry. The overall objective of the present investigation is to develop a second-order diffraction method in the time domain which avoids the complicated algebra required by most conventional second-order methods and allows studies on deterministic irregular wave diffraction without solving the problem at each time step. More important, the present method may be extended to include the effects of arbitrary body motions or manoeuvres without any serious modification to the approach. It is hoped that this thesis wil l help provide designers of offshore structures a simple, effective and comprehensive numerical method for nonlinear analyses of regular or irregular wave diffraction around large offshore structures. 7 CHAPTER 2 THEORETICAL FORMULATION 2.1 Full Nonlinear Method The full nonlinear boundary value problem defining the fluid motion is first considered. With reference to figure 1, the three-dimensional problem is defined with a right-handed Cartesian coordinate system (x,y,z), in which x and y are measured horizontally and z is measured vertically upwards from the still water level. The body located at the centre of the domain is rigid and fixed in space, and the body surface is impermeable. The incident waves are two-dimensional in the x-z plane and are progressive in the positive x direction. In the special case of a two-dimensional body in the x-z plane, the flow does not depend on the y axis and so the y coordinate can be omitted. Let t denote time and T| the free surface elevation above the still water level. The seabed is assumed impermeable and horizontal along the plane z = -d . With the fluid assumed incompressible and inviscid, and the flow irrotational, the fluid motion can be described by a velocity potential <) such that the fluid particle velocity u at any point in the fluid domain 3is given by • t o * <2i> and that the velocity potential <)> satisfies the Laplace equation within the fluid domain: The potential <|) is subject to boundary conditions on the seabed, the instantaneous wetted body surface S w and the free surface Sf at z = rj, which are given respectively as 8 at z = -d (2.3) 8j> 3n = 0 on S w (2.4) d§ drj 3<j> dr\ 8<j> dr) 3z - dt ~~ dx dx - dy 3y on Sf (2.5) f + g T l + i | V ( t , | 2 = P a on Sf (2.6) Here g is the acceleration due to gravity; P a is the atmospheric pressure; n denotes distance in the direction of the unit normal vector n = (nx,ny,nz) directed outward from the fluid region. Equations (2.3) and (2.4) correspond respectively to the kinematic boundary conditions on the impermeable seabed and body surface where the normal fluid velocity component vanishes. Equation (2.5) is the kinematic free surface boundary condition requiring that the fluid particle velocity normal to the free surface is equal to the velocity of the free surface itself in that direction, while (2.6) is the dynamic free surface boundary condition stating that the pressure at the free surface expressed in terms of the Bernoulli equation is equal to the atmospheric pressure P a thereby neglecting the effects of surface tension. In the present application, the atmospheric pressure can be arbitrarily set to zero without inducing any loss of generality. In addition, a suitable radiation condition is to be applied on the control surface S c some distance from the body as indicated in figure 1. The major difficulties associated with this problem arise primarily from the poorly defined radiation condition involving transient waves and from the two nonlinear free surface boundary conditions applied on the instantaneous free surface which itself is unknown a priori. A boundary integral method involving a Green's function may be used to solve the boundary value problem. Since the potential <|> is a harmonic function, the second form of Green's theorem may be applied over a closed boundary containing a fluid region, relating the 9 values of $ and its normal derivative 3<j>/9n along the boundary. The potential <}>(x) at a point x is expressed as <Kx) = ± J[G(x,^)^) - «KS)^k.W dS s (2.7) Here \ represents a point (£, \|/, Q on the surface S over which the integration is performed and n is measured from the point The constant a depends on the location of !; and on the dimensions of the problem. G(x,£) is a Green's function which can be interpreted as the potential at any field point x due to a source of unit strength at % and which satisfies the Laplace equation in the domain except at the point x = where 8 a dirac-delta function. In the present context the surface S would comprise of the wetted body surface S w , the instantaneous free surface Sf, the control surface S c surrounding the body and the seabed, as indicated in figure 1. If initial conditions are chosen to correspond to still water in the vicinity of the body and a known incident wave train approaching it, the development of the flow together with the associated fluid forces can be obtained by a time-stepping procedure in which the flow at each time step is calculated by solving the integral equation (2.7). The formulation so far provides the basis for an exact solution to the boundary value problem and a solution procedure has been described by Isaacson (1982). Since the initial waves in the incident wave train are unsteady, this approach requires a relatively large number of time steps for the fully developed flow to be established at the cylinder. A n alternative initial condition is to treat the incident wave train as initially regular everywhere in the computational domain so that the body imposes its presence on the fluid just after the initial time (Isaacson & Zuo, 1989). Normally, the steady state solution for the wave force calculation is obtained in less than a wave period. V2G(x,S) = 8(x-4)8(y-v)8(z-C) in 3 (2.8) 10 Although in principle, the full nonlinear method would give a more accurate solution, it requires a lengthy numerical calculation in which a new system of simultaneous equations must be generated and solved at each time step as the free surface moves to a new position. Numerical errors introduced into the procedure in solving the system of simultaneous equations and in predicting the new position of the free surface, accumulate through the time-stepping procedure and provide a further limitation to the length of simulation time. In addition, the absence of a proper radiation condition for the full nonlinear diffraction problem makes it more difficult to achieve a stable, steady state solution over many wave periods. 2.2 Stokes Second-Order Method 2.2.1 Perturbation Expansion When the water depth is not small compared with a typical wavelength and when the wave height is moderate, it is possible to reduce the two nonlinear free surface boundary conditions originally derived on the instantaneous free surface to conditions evaluated at the still water level by the Stokes second-order expansion procedure. The variables describing the flow are developed as power series in terms of a perturbation parameter. It is expected that under the required conditions the series will converge and provide a good approximation to the full nonlinear problem as the number of terms increases. As the location of the free surface is unknown a priori, the full nonlinear free surface boundary conditions (2.5) and (2.6) are first expanded about the still water level with a Taylor series expansion: + ... = 0 (2.9) z=0 + ... = 0 (2.10) 11 The problem may now be cast in a time-independent geometry which includes the still water surface S 0 in place of the instantaneous free surface Sf. To proceed with the Stokes second-order expansion, components at first and second order in the formulation are separated by introducing perturbation expansions for <() and TJ: $ = e2<|>2 + - (2.11) 1 1 = 6 ^ ! + e 2 r |2 + - (2.12) where e is a perturbation parameter related to the wave height which is small and the subscripts 1 and 2 indicate components at first and second order respectively. Each subscripted variable appearing in the series is taken as having the same order of magnitude, and thus each additional term in the series represents a quantity smaller than the preceding one by a factor of e. It should also be noted that the free surface boundary conditions at first and second order are now applied on the still water surface. The potentials and free surface elevations at first and second order are further separated into incident and scattered components: <|> = eCtf + tf) + e2(<^ + ^ ) +••• (2.13) •TI = e (nT" + n si) + e 2 (Tl2 +T12 ) + - (2-14) The superscript w indicates the incident components which satisfy the Laplace equation and all boundary conditions except on the body surface, while the superscript s indicates the scattered components which are to be evaluated. For the case of regular incident waves with height H , length L and period T, Stokes second-order wave theory can be applied to give the incident potentials at first and second order respectively as w 7iH cosh [k(z+d)] . .. . ,~ . C N ^ = k T sinh(kd) " ° 0 » - < t t ) <2-15> 12 (2.16) and the free surface elevations at first and second order respectively as -~- cos (kx-cot) (2.17) H ,7CrL cosh (kd) 8 L ^ sinh3 (kd) [2 + cosh (2kd)] cos [2(kx-cot)] (2.18) where k = 27C/L is the wave number and co = 27t/T is the angular frequency. The velocity potential at second order also contains an arbitrary function of time. Since this does not contribute to wave forces and kinematics at second order, it need not be considered here. The Ursell parameter U = HL 2 / d 3 is useful as an indication of the relative importance of second-order effects and the validity of the solution at second order. It is noted that the incident wave components at second order act at twice the first-order wave frequency and have a wavelength equal to one-half the first-order wavelength. These second-order wave components propagate at the same speed as the first-order waves and are phase-locked with the first-order wave system to produce an overall wave profile with steeper crests and flatter troughs. Substitution of the power series representation for <]> and rj into the Laplace equation and boundary conditions, including the expanded free surface boundary conditions, gives rise to separate boundary value problems for each of the e and e2 terms in the power series. Although solutions to any order can be formulated in principle, in the present study only solutions to second order are considered. With the incident potential and free surface elevation at k-th order given in terms of the perturbation parameter as £k<j>k and ekT|k respectively, the solutions for the scattered components should then correspond to &fa and £kT|k respectively. Since the perturbation parameter in the form of ek is attached to every term in the governing equations at k-th order, for simplicity the perturbation parameter is not written explicitly in the governing equations and from now on $ and fa are used to denote e 1 ^ and £*fa respectively. 13 2.2.2 Three-Dimensional Problem The formulation in three dimensions is first considered in this section, and simplifications are subsequently made for the special case of a two-dimensional problem in the next section. Since the incident potential already satisfies the Laplace equation and most of the boundary conditions, it is excluded from most of the governing equations. Its effects are accounted for by the body surface boundary conditions at first and second order and by the free surface boundary conditions at second order, which provide the basis for the development of the scattered potential. The solution to the flow problem is then obtained by a superposition of the incident and scattered potentials. The boundary value problems at first and second order are now applied to a time-independent domain D bounded by the seabed, the body surface below the still water level Sb, the control surface S c and the still water surface S Q . At first order, the Laplace equation, the seabed boundary condition, the body surface boundary condition, the kinematic and dynamic free surface boundary conditions are given respectively by V% = 0 in D (2.19) at z = -d (2.20) ,s w 9n 8n on Sb (2.21) (2.22) s (2.23) 14 The first-order problem is characterized by the two linear homogeneous free surface boundary conditions (2.22) and (2.23) applied on the still water surface, and the external excitation to the hydrodynamic system is applied through the body surface boundary condition (2.21). The solution to the scattered potential is homogeneous and is composed of free waves at the incident wave frequency. For the case of regular wave propagation, the angular frequency and wave number have to satisfy the linear dispersion relation co2 = gktanh(kd) (2.24) analogous to the free vibration of a single-degree-of-freedom system. At second order, the Laplace equation, the seabed boundary condition, the body surface boundary condition and the kinematic and dynamic free surface boundary conditions are given respectively as V2<)>! = 0 in D (2.25) ^ = 0 at z = -d (2.26) f = o „ S b (2.27) "3T-~3r = -3T + "3r + "3T"ST + WW ~ n'~s* ° ( ' ^ + CTt2 = - ^ - g n I - ^ I V < . 1 P - r 1 1 ^ on So (2.29) The problem at second order is characterized by the two inhomogeneous free surface boundary conditions (2.28) and (2.29) applied on the still water surface, and the external excitation to the hydrodynamic system is applied through the body surface boundary condition (2.27) as well as 15 the two free surface boundary conditions. The right-hand sides of the two free surface boundary conditions are composed of the known incident wave components at second order and quadratic terms which can be evaluated from the first-order solution. The quadratic terms in the dynamic free surface boundary condition (2.29) physically represent an oscillatory pressure field at double the incident wave frequency applied to the still water surface which would generate second-order forced waves subject to the geometric constraint defined by the quadratic terms in the kinematic free surface boundary condition (2.28). Since the solution at first order includes the sum of the incident and scattered waves, the second-order forced waves generated by the quadratic terms are then composed of three components: (1) a plane wave component associated with the self-interaction of the first-order incident waves; (2) a scattered wave component associated with the self-interaction of the first-order scattered waves; and (3) a third component associated with the cross-interactions between the first-order incident and scattered waves. The plane wave component (1) is part of the incident waves and is excluded from the scattered potential. Its effects on the solution at second order are then accounted for by the body surface boundary condition (2.27). These second-order forced waves are phase-locked with the first-order wave system and do not satisfy any simple dispersion relation. The pressure fields imposed on the still water surface by the forced wave components (2) and (3) are non-uniform in space and are concentrated near the body. The motion of these two components of forced waves on the still water surface would generate second-order free waves at double the incident wave frequency. Imposing the body surface boundary condition, the diffraction of all three components of forced waves gives rise to additional second-order free waves. These second-order free waves 16 propagate at the corresponding celerity independent of the first-order wave system and satisfy the linear dispersion relation in the form of (2co)2 = gk 2 tanh (k2d) (2.30) where k 2 is the wave number of the second-order free waves. In deep water, k 2 approaches 4k while in shallow water k 2 approaches 2k. The scattered potential at second order is then composed of the forced wave components (2) and (3) and free waves generated by all three components of forced waves. Mathematically these forced and free waves represent respectively the inhomogeneous and homogeneous solutions of the boundary value problem. With the problem now defined in terms of scattered components, a radiation condition applied to the scattered potential can now be introduced. Physically, this condition requires that the scattered waves must propagate away from the fluid domain. The Sommerfeld radiation condition, originally derived on the basis of spatially and temporally periodic conditions of a potential function, has been extended by Orlanski (1976) so as to account for unsteady wave motions. In the time domain, the Sommerfeld radiation condition for the scattered potential at k-th order is given by: where c is the celerity of the scattered waves and r is radial distance measured from the body. This has been modified by Orlanski (1976) so as to apply on a control surface S c located at a finite distance from the body (see also Chan, 1977; and Yen & Hall, 1981). The term Vr is no longer necessary, and the condition is applied in the direction normal to the control surface: ( k = l , 2 ) as r -> oo (2.31) = o ( k = l , 2 ) on S c (2.32) dt where c n is the time-dependent celerity of the scattered waves on the control surface S c in the direction of the normal n. 17 However, some attention must be given to the basis for adopting (2.32). In order to examine this, the scattered wave field in three dimensions and at k-th order is expressed in the form: = A k(x,y,z) F k (k x x + k y y - cot) (k = 1, 2) in D (2.33) where k x and k y are local wave numbers in the x and y directions respectively, corresponding to wave propagation in an arbitrary direction. The functions Fk define the spatial and temporal periodicity of the wave field and the functions A k represent the amplitudes of the potential which in general vary with x, y and z. Differentiating (2.33) with respect to the direction n normal to the control surface gives rise to <*-i-*> ° " s < <2-34> which involves an additional term compared to (2.32). However, the amplitude of the scattered waves at first order and the lead amplitude of the scattered waves at second order are both expected to decay as Vr, where r is distance from the body. The third term in (2.34) involving 9Ak/9n may thereby be shown to be of order (1/kr) smaller than the first two terms, and so becomes negligible when the control surface is located at a sufficiendy far distance from the body. (For example, on the basis of the closed-form first-order solution for a vertical circular cylinder, the third term may be up to 5% of either of the other two terms when S c is located about one wavelength away from the body.) Thus, provided that S c is located sufficiendy far from the body, the radiation condition given by (2.34) thereby reduces to the more common form given by (2.32). For diffraction of regular waves at first order, if a large single body such as in the case of a vertical circular cylinder is considered, it is convenient to introduce an additional approximation that the scattered waves at a large but finite distance from the body propagate in the radial direction and that their speed c is known on the basis of the linear dispersion relation. At 18 second order, the scattered potential is made up of components of forced and free waves propagating at different celerities and the radiation condition cannot be satisfied with a unique celerity. Intuitively the scattered potential at second order on the control surface may be separated into forced and free components so as to apply the radiation condition separately to account for the different propagating speeds. However, this approach is only effective for the free wave component and is untenable for the forced waves. The propagation characteristics of the forced waves are quite complicated and, as noted by Mei (1983), the celerity of the forced wave component (3) in the radial direction ranges from c = co/k to infinity, and in the present application it would further vary with time before a steady state condition is obtained. Since the problem appears not be simplified by separating the scattered potential at second order into components, it is more efficient to explore the usefulness of a single radiation condition applied to the entire scattered potential at second order. Based on Orkanski (1976), the Sommerfeld radiation condition is first applied to estimate the time-dependent celerity on the boundary of the flow which in turn is used to calculate the boundary conditions on the control surface for an advanced time. This approach assumes a locally and momentarily periodic condition of the flow on the control surface so that a celerity exists to define the outflow of the scattered waves instandy. By applying the time-dependent celerity in the Sommerfeld radiation condition, an irregular wave train is literally treated and transmitted as a series of periodic waves in time. This approach has been successfully applied to treat the radiation problem of full nonlinear waves by Chan (1977), Yen & Hall (1981) and Jagannathan (1988) and is applied here to treat the diffraction problem to second order. The problems at first and second order may now be defined for a time-independent domain such that the scattered potentials at first and second order in turn satisfy the Laplace equation within the domain and are subject to the boundary conditions on its boundary. The scattered potential at k-th order at a point x located on a smooth part of the boundary is given by the following boundary integral equation involving a Green's function G: 19 «fc ( x ) = j - j[G(x^^{%) - • J f l ^ x , * ) ] dS (k = 1, 2) (2.35) Here £ represents a point (£,y,C) on the surface S' over which the integration is performed and n is measured from the point %. For x located inside the boundary, the factor 1/2TC is simply replaced by 1/4TC. In the present context the surface S' would comprise of the body surface Sb below the still water level, the still water surface S 0 , the control surface S c surrounding the body and the seabed. However, with the body and subsequently the wave field symmetric about the x axis and the seabed horizontal, it is more efficient to simulate one half of the total surface with the seabed excluded from S' and to choose a Green's function that accounts for the double symmetry about the x axis and the seabed. This is GM) = £ ^ (2-36) k=l r k where r k is the distance between the points x and The source points £ k are defined such that£i = (cj, Q is a point on S'; %2 = {%, -\\f, C) is the image of c%\ about the x axis; £3 = (cj, -\|/, -(£+2d)) is the image of %i about the seabed; and 4^ = (£, \\f, -(£+2d)) is the image of £1 about the seabed. Thus, r k is given by r [(S-x)2 + ( y - y ) 2 + (C-z) 2 J 1 / 2 (k = 1) [(S-x) 2 + (v+y)2 + (C-z) 2] " 2 (k = 2) (2.37) [ ( £ - x ) 2 + (\j/+y)2 + (C+2d+z) 2] 1 / 2 (k = 3) [ ( £ - x ) 2 + (\|/-y) 2 + (C+2d+z) 2] 1 / 2 (k = 4) With the subscript k ranging from one to four, the integral equation (2.35) actually represents an integration of the surfaces in the four quadrants of the doubly symmetric geometry. r k = < 20 2.2.3 Two-Dimensional Problem A simplification to the formulation described above may readily be adopted to treat the corresponding two-dimensional problem in the vertical x-z plane. With the y coordinate omitted and the normal vector defined in the x-z plane, all boundary conditions remain unmodified. The scattered waves at first and second order now propagate in either the positive or negative x direction. Although the propagating characteristics of the scattered waves at second order appear to be simpler here, the amplitude and nonlinearities of the scattered potential persist to the control surface and thus require a more accurate treatment of the radiation condition. In adapting the method used here, the surface integral equation deriving from Green's theorem reduces to a line-integral equation in which the influence of the co-ordinate y is absent. When x is located on a smooth part of the boundary, the integral equation for the potential at k-th order ((£ is then given as Here I; represents a point on the surface S' over which the integration is performed and n is measured from the point For x located inside the boundary, the factor — i s simply replaced by -1/2JC. In the present context the surface S' would comprise of the body surface Sb below the still water level, the still water surface S 0 , the control surface S c surrounding the body and the seabed. However, because of the assumption of a horizontal seabed, the seabed can be excluded from S' and a Green's function can be chosen to account for the symmetry about the seabed. This is s ( k = l , 2 ) (2.38) 2 G(x£) = I In (rtf (2.39) 21 where r k is the distance between the points x and £ k . The source points 4k are defined such that %i = (£, Q is a point on S'; and £2 = -(C+2d)) is the image of <%i about the seabed. Thus, r k is given by '[(4-x) 2 +(C-z)2]l/2 ( k = l ) r k = (2.40) . [(£-x)2 + (C+2d+z)2] 1/2 (k = 2) With the subscript k ranging from one to two, the integral equation (2.38) actually represents an integration of the surface S' and its image about the seabed. 2.2.4 Time-Integration and Initial Conditions In the preceding formulation, the boundary conditions at first and second order are defined on the body surface, the control surface and the still water surface, and the solution is expressed in terms of an integral equation relating the scattered potential and its normal derivative on the boundary. For a well-posed problem either the scattered potential or its normal derivative on each portion of the boundary is known and is evaluated from the corresponding boundary condition. The solution to the integral equation can therefore be obtained by a numerical procedure from which the remaining unknowns can be evaluated. In the present application, the normal derivative of the scattered potential on the body surface is known and is given explicitly by the body surface boundary condition, (2.21) or (2.27). Although the normal derivative of the scattered potential is also given in the radiation condition (2.32) and in the kinematic free surface boundary condition (2.22) or (2.28), its value depends on a temporal derivative which is unknown before the solution at an advanced time is obtained. On the other hand, the scattered potential on the control surface and the still water surface at an advanced time can be evaluated by a time-integration of the temporal terms respectively in the radiation condition and the dynamic free surface boundary condition, (2.23) 22 or (2.29). With the flow field known at time t, the scattered potential at k-th order on the control surface and the still water surface at an advanced time (t+At) is given by t+At <t>sk(t+At) = ^(t) + J ^ dt (k = 1,2) on S c and S 0 (2.41) t in which d$Jdt is evaluated from the corresponding boundary condition on".Sc or S 0 . The time-integration of the dynamic free surface boundary condition requires a knowledge of the scattered free surface elevation which can be obtained by a time-integration of the kinematic free surface boundary condition. Likewise the scattered free surface elevation at k-th order at an advanced time (t+At) is given by t+At Tisk(t+At) = Tfk(t) + j^ dt (k = 1,2) on S 0 (2.42) t in which 8rjk/3t is evaluated from the kinematic free surface boundary condition. Various numerical techniques are available to evaluate the integrals in (2.41) and (2.42) and are discussed in the next chapter. To shorten the transient duration in the simulation, a fully developed initial condition similar to that of Isaacson & Zuo (1989) is adopted here. At t = 0, the scattered potential is identically zero and the body has no effects on the flow. The initial conditions are known and correspond to an undisturbed Stokes second-order wave field in the domain. With the body surface boundary condition imposed for t > 0, scattered waves at first and second order are generated on the body surface and propagate into the computational domain. With the flow field at any time instant obtained by solving the integral equation (2.35) or (2.38), the development of the flow in time is obtained by the time-integration of the flow parameters in (2.41) and (2.42). The simulation is continued until a steady state condition is indicated. 23 2.3 Alternative Second-Order Method 2.3.1 Second-Order Expansion of Integral Equation With reference to the full nonlinear problem defined in § 2.1, both the nonlinear free surface boundary conditions and the integral equation are defined on the instantaneous free surface. In the Stokes second-order method discussed in § 2.2, only the free surface boundary conditions are expanded about the still water level and the integral equation is then applied to a stationary boundary which includes the still water surface. On the other hand, the integral equation for the full nonlinear problem may consistently be expanded together with the free surface boundary conditions to be defined on the still water surface. In order to study the effects of this alternative expansion on the predictions of the second-order solution and to simplify the formulation and computation, the special case of a two-dimensional problem in the x-z plane is considered for this approach. In order to highlight the significance of this alternative expansion, the integral equation for the full nonlinear problem in two dimensions is rewritten from (2.7) without approximation as: Here the potential at a point x on the boundary is given by an integral over the fixed boundary defined by S', together with two correction integrals which together account for the fluctuation of the free surface about the still water level and the use of a normal defined with respect to the instantaneous free surface. Literally, the Stokes second-order method expands the nonlinear free surface boundary conditions so as to evaluate the potentials at first and second order on the s (2.43) 24 still water surface. Together with other boundary conditions known on the rest of the boundary, the solutions at first and second order are obtained from the application of the first term of (2.43). In addition to making errors of third order in the context of the perturbation expansion, these solutions neglect the last two terms of (2.43) which may be just as important. In order to reduce the second integral in (2.43) to integrals evaluated at the still water level, the values of $ and G on the instantaneous free surface may consistendy be expanded by a Taylor series about z = 0. Retaining terms up to the second approximation, the following expansions are made: 4>® = [ ¥ £ ) + T l 6 ) f k ) k o (2-44) '3C G(x,4) = [G(x,4) + r,(x)^(x ,4) + rj(4)| |(x ,4)] 2 = 0 ; = 0 (2.45) and dS in the second integral of (2.43) is approximated by dS = 4^ " on S 0 (2.46) where dS is defined on the instantaneous free surface and dS' is defined on the still water surface which is independent of time. The normal derivative of any variable at % can be expressed as a a a ._. 3 - = n?— + n c — (2.47) d n ^ a^c If % is located on the free surface, the terms n^ and n^ can consistendy be expanded in terms of the gradient of the free surface elevation by n? = -e at on So ( 2 - 4 8 ) n^ = 1 - i e 2 ( ^ ) 2 + on SG (2-49) 25 It is understood that in (2.45) rj(x) and r\(%) are taken as zero for x and % respectively not on the free surface. Since G is a function of (z-Q, dG/dz and BG/dC, are equal and have opposite signs. For x = %, the last two terms of (2.45) cancel each other and the Green's function remains unexpanded at this singular point. Similarly, when the two points are both located on the free surface and are in close proximity, the effect of the expansion is small which can be justified on physical grounds. On the other hand, when the two points are far apart, the expansion is insignificant due to the decaying properties of the Green's function. The effect of the expansion is most significant when either one of the two points is not on the free surface and when they are in close proximity. 2.3.2 Alternative Second-Order Integral Equation Thus by substituting (2.11), (2.12) and (2.44) - (2.49) into (2.43) and collecting terms of order e and e 2 in turn, the first- and second-order integral equations evaluated over the time-independent boundary S' are obtained. The integral equation for the first-order problem is given by the same equation as in the Stokes expansion method (2.38), while the second-order integral equation is modified to: <|>2(x) = - i 7C 71 1 S ' 1 Ti i (5 ) [G («^)^<6) - <MS)§oa)]ds (2.50) 26 Comparing (2.43) and (2.50), the two correction integrals in (2.43) taken together are of second order and are approximated by the last three terms in (2.50). More specifically, the second and third integrals in (2.50) account for the second-order effects due the fluctuation of the first-order free surface and the last integral in (2.50) accounts for the normal defined with respect to the first-order free surface. In other words, the instantaneous free surface Sf in the second term of (2.43) is now approximated by the first-order free surface for a second-order approximation. Although the same free surface boundary conditions applied in the Stokes second-order method also apply here, the relation between the potential and its normal derivative at second order is modified due to the expansion terms in (2.50) and a Stokes second-order incident wave field is no longer applicable here. Since a known incident potential at second order is not available for this approach, incident waves to second order are generated in the numerical procedure. With a still water initial condition, an incident wave train can be generated by applying the appropriate excitations on the upwave control surface, analogous to that of a wave maker in a wave flume. The incident waves then propagate into the computational domain and interact with the body. 2.4 Wave Runup and Forces 2.4.1 Free Surface Elevations and Runup One of the primary objectives of this study is to determine the second-order free surface elevation and the runup around a large body. Despite previous efforts in solving the second-order diffraction problem, only a limited number of studies have been directed towards the prediction of second-order runup (e.g. Kriebel, 1897; and K i m & Yue, 1989), and comparisons between theoretical results are scant. 27 With the velocity potential known, the free surface elevation to second order can be given explicitly by the dynamic free surface boundary condition ^ = 4 ^ + 2 L | V * 1 ' 2 + o n S ° ( 2 - 5 1 ) The numerical values of 3<p/3t are obtained by applying a central difference approximation to in time and the numerical values of V<J>i are obtained by a spatial interpolation involving neighbouring facets. For regular wave diffraction to second order, the free surface elevation may be expressed in the usual way as a sum of three components TI = r i ! + T|0 + T|2 (2.52) in which r^, T|0 and TI2 are respectively the first-order oscillatory component at the incident wave frequency, the second-order steady component and the second-order oscillatory component at the twice the frequency. The expressions for Tj^ r i o and TJ2 are given respectively as: T l ^ - ^ C ^ 1 ) on So (2.53) Tlo = -|4lV<t>il 2 + i \ ^ > on S 0 (2.54) Tl 2 = -|(lJr + i l V <W 2 + HiiJ}) - ^0 on So (2.55) where < > denotes a time average. The component r j 0 represents a spatially-varying, corrugated mean water surface due to the non-zero mean of the second-order forced terms in the free surface boundary conditions and is an order smaller than the first-order wave amplitude. The component TI2 represents the free surface elevation at second order which is composed of second-order forced and free waves. 28 With the free surface elevation to second order evaluated, the runup R and rundown R' to second order can readily be obtained by R = Ol i+ i lo + T h W a t w o (2-56) R' = 0 l i+T | o + Tl2)min at w 0 (2.57) where w 0 is the water-line contour associated with the still water condition and the subscripts max and min represent respectively the maximum and minimum values of the free surface elevation in time after a steady state condition is obtained. In the two-dimensional problem, the runup applies only at the two intersection points between the still water surface and the body surface. When the body surface is inclined at the still water level, the motion of the free surface at w 0 is no longer restricted in the vertical direction. The Taylor series expansion in (2.9) and (2.10) which is formally valid for waves interacting with a vertical surface fails to account adequately for the actual fluid motion near w 0 . Although this does not prevent the existence of a viable solution to second order such as that of Kim & Yue (1989), results such as local pressure or free surface elevation near w 0 might not be accurate and must be interpreted more carefully. However, to be consistent with the second-order approximation it is expected that the errors in the solution would be relatively insignificant for bodies with relatively large angles of surface inclination. 2.4.2 Wave Forces and Moments The second application of the second-order solution is the determination of the local pressure acting on the body surface and subsequendy the depth integrated forces. In potential flow, the pressure distribution over the body surface may be determined by the unsteady Bernoulli equation, 29 where p is the fluid density. Once the potentials at first and second order have been obtained, the wave forces and moments acting on the body can be determined by carrying out an integration of the pressure over the instantaneous wetted body surface S w . The force vector is then given as F = J p n' dS (2.59) s w in which n' in the two- and three-dimensional problems are given by (n x , n z , z n x - x n z ) n' = -. (n x , n y , n z , y n z - z n y , z n x - x n z , x n y - y n x ) and F in the two- and three-dimensional problems are given by (2.60) ( F x , F z , M y ) F = i (2.61) (Fx, F y , F z , M x , M y , M z ) In two dimensions, the force vector F represents two components of wave forces in the x and z directions and a moment about the y axis. In three dimensions, the force vector F represents the three components of wave forces and moments in the three translational and rotational directions respectively. To be consistent with the second-order expansion in § 2.2.1, (2.59) defined on the instantaneous wetted body surface can be expanded about the still water level to give an integral defined on the body surface below the still water level together with a correction integral defined at the still water-line contour which accounts for the fluctuation of the free surface about the still water level, 30 F = - J p ( ^ + \ I V O i l V d S _ J j P ( ^ L + gzVdzdw (2.62) where w denotes distance measured along the water-line contour. In the two-dimensional problem, the integration along the still water-line contour is reduced to two terms evaluated at the two intersection points between the body surface and the still water surface. It is also assumed that the body surface intersects the still water surface perpendicularly. If the body surface is inclined at the still water-line contour, an extra factor, (1 - n z ) _ 1 / 2 , should be included in the water-line integral to account for the vertical force that acts on the body surface between the free surface and the still water level. Again this is valid under the assumption of large inclination angles (i.e. n z is small compared to one). Substituting dfyi/dt = - gr\i from the first-order dynamic free surface boundary condition (2.23) into the second term of (2.62) and carrying out the required integration over 0 < z < gives The second integral in (2.63) accounts for the second-order wave force component due to the fluctuation of the first-order free surface about the still water level and represents a resultant between the first-order hydrostatic and dynamic pressures (Lighthill, 1979). For regular wave diffraction to second order, the total force acting on the body may be expressed in the usual way as a sum of three component forces: where F i , FoandF2 are respectively the first-order oscillatory force at the incident wave frequency, the second-order steady drift force and the second-order oscillatory force at twice the wave frequency. The expressions for F i and FQ are given respectively as: (2.63) F = Fi + FQ + F 2 (2.64) 31 F i = - p J^-n' dS (2.65) Sb F 0 = - - p < J IV^! | 2n* dS> + | p g < J t i f n ' d w > (2.66) 2 si w 0 The second-order oscillatory force F 2 is further composed of two components, F ^ and F ^ which are due to <(>i and tyi respectively, = - - p J IV^i I2 n' dS + jpg Jnjn'dw - F 0 (2.67) 2 Sb w ° F ? = - p J ^ n ' dS (2.68) Sb The maximum force and moment to second order which are of practical interest can be evaluated direcdy from (2.64). 32 CHAPTER 3 NUMERICAL FORMULATION 3.1 Field Solution In Chapter 2, the field solution at any time instant is expressed in term of an integral equation relating the potential and its normal derivative on the boundary. A standard method to solve the integral equation is to discretize the boundary into a series of facets over which the potential and its normal derivative are assumed to vary according to an interpolation function. The integrals over each facet are then evaluated by a numerical integration in terms of the potential and its normal derivative at a number of nodes where the boundary conditions are applied. The application of the discretized equation to all the nodes in turn gives rise to a system of linear algebraic equations from which a numerical solution can be obtained. To simplify the numerical formulation and subsequently the computer codes implementing the procedure, the potential and its normal derivative are taken as constant over each facet and applied at the facet centroid. To account for the rapid variation of the Green's function when the source and field points are close together and to speed up the convergence of the numerical solution with increasing number of facets, the integrals involving the Green's function and its normal derivative are evaluated more accurately by a Gaussian integration. 3.1.1 Stokes Second-Order Solution The incident wave field to second order for any waveform is known for the Stokes second-order method and is simply given by a Stokes second-order wave theory (e.g Hasselmann, 1962; and Longuet-Higgins, 1963). For regular incident waves, the incident potential and free surface elevation to second order are given by (2.15) to (2.18) respectively. At t = 0, the 33 scattered potential is identically zero, and the initial conditions are known and correspond to an undisturbed Stokes second-order wave field in the domain. The scattered potential to second order is then developed in time and space through the imposition of the body surface boundary condition and the simulation is continued until a steady state condition is indicated. Three-Dimensional Problem The integral equation (2.35) for the three-dimensional problem is solved by a numerical procedure in which the boundaries Sb, S c and S 0 are discretized into finite numbers of planar quadrilateral facets. The corresponding values of and o^/Bn are then taken as constant over each facet and applied at the facet centroid. The integral equation involving the scattered potential at k-th order may then be approximated as (i = 1, N ; k = 1, 2) (3.1) in which i and j are the indices associated with x and £ respectively; N is the total number of facets in the formulation; and ASj is the area of the j-th facet. It is understood that the surfaces are treated in the order of Sb, S c and S 0 as the indices i and j increase from 1 to N . Equation (3.1) corresponds physically to the influence of the boundary values at all facets from j = 1 to N on the boundary value at the i-th facet. The application of (3.1) to all the facets in turn, with i taken to range from 1 to N , gives rise to the following system of equations: (i = 1, N ; k = 1, 2) (3.2) in which Nb is the number of facets on Sb- The coefficients Bjj and Ajj correspond respectively to integrals of the Green's function and its normal derivative over area of the j-th 34 facet as in (3.1). When i = j , the coefficients are evaluated by a closed form integration and when i * j , the integration is performed by a Gaussian numerical integration. The coefficients Ajj and By are thereby approximated by: Ai i = < 1 4 ™ M (4 k m n-x)-n k — I I I J W m w n 27C k=l m=l n=l rkmn - 1 4 . 1 4 V V T (5kmn-x)-nk + —- I I I J w m w n — 3 27C k=2 n=l m=l rkmn (i * j) ( i= j ) (3.3) Bii = < I 4 M M j — I I 1 / w m w n - — 27C k=l m=l n=l r k m n 1 4 M M l T i + — I I I / w m w n 2K k=2 m=l n=l rkmn (i * j) (i = j) (3.4) where m and n are the indices of the Gauss points in the k-th image of the j-th facet; M is the number of Gauss points to be used in the numerical integration in one direction; and w m and w n are the weight functions at each of the Gauss points. The term Tj in (3.4) represents the integration of the Green's function over the area of the i -th facet when i = j and k = 1, AS; dS ( i = 1, N) (3.5) In principle, a Gaussian integration can consistently be applied to evaluate (3.5). However, the numerical integration is inaccurate near the singular point r = 0 and a substantial number of Gauss points are required to yield satisfactory results. On the other hand, a closed form solution of the integral for any polygonal facet shape has been given by Hogben & Standing (1974) and is adopted here. For an n-sided polygon with vertices at (s^ti), (s n,t n) in a local coordinate system (s,t) and each making angles 6 1 , 0 n in turn with the s-axis at a origin, T is given as 35 2% k=i K S cos(8/2) + T sin(0/2) - R cos(0/2) • -. ek+i ] S cos(9/2) + T sin(8/2) + R cos(6/2) J 6 k (3.6) in which S = s k + i - sk, T = t k + i - t n, Q = s k + 1 t k - SJA+I and R = (S 2+T 2) 1/2. The term / in (3.3) and (3.4) is known as the Jacobian which relates the global and local coordinate systems of a facet in the numerical integration. Its value is given by the determinant of the Jacobian matrix evaluated at the location of each Gauss point. For a planar quadrilateral facet, the Jacobian is given by / = det in which N i = < i=l 3s i 8 N i i=i 3s ti - i=i 3t i=l dt ti | d - s ) ( l - t ) (i = 1) ^ ( l + s X l - t ) (i = 2) ^( l+sXl+t) (i = 3) ^ ( l - s X l + t ) (i = 4) (3.7) (3.8) where Ni is the shape function relating the global and local coordinate systems of a facet. For a rectangular facet, / is constant and is equal to ASj/4. Two-Dimensional Problem The integral equation (2.38) for the two-dimensional problem is solved by a similar numerical procedure in which the boundaries Sb, S c and S Q are discretized into finite numbers of facets made up of straight-line segments. The corresponding values of fa and 3())k/3n are 36 then taken as constant over each facet and applied at the facet centre. The numerical discretization gives rise to the same system of equations as in (3.2) but the coefficients Ay and By are evaluated by different expressions on account of the different integral equation and Green's function that are now applicable. These are now given as A;: = < 1 2 M n t - C k n - x ) - I Z J w m — 2 TC k=l m=l r k m , 1 " J ! n2-(42m-x) 1 - - I / W m — 2 ^ 7C m=l r 2m (i*j) (i=j) (3.9) r i 2 M - I Z / w m l n ( r k m ) 7C k=l m=l Bii = < I M Ti + - I / wm In (r2m) 7C m=l ( i ^ j ) ( i = j ) (3.10) where m is the index of the Gauss points in the k-th image of the j-th facet; M is the number of Gauss points used in the numerical integration; and wm is the weight function at each of the Gauss points. The term in (3.10) represents the integration of the Green's function over the length of the i-th facet when i = j and k =1, and is given by a closed form integration as ASi/2 r , . i J l n | r | d r - ^ [ l n ( ^ i ) - l ] n -ASj/2 (i = 1, ..... N) (3.11) 7C For facets made up of straight-line segments, the Jacobian / for the j-th facet in (3.9) and (3.10) is simply given by AS; J = (3.1?) which is a constant within each facet 37 System of Simultaneous Equations Equation (3.2) represents a system of N equations with 2N variables corresponding to the scattered potential and its normal derivative at each of the N facets. For a well-posed problem, either one of the two variables at each facet is known and a system of simultaneous equations with N equations and N unknowns is obtained. In the present application, the input vector on the right-hand side of the simultaneous equations (3.2) consists of the normal derivative of the potential on the body surface which is known from the body surface boundary condition and the potential on the control surface and the still water surface which is to be evaluated from a time-integration of the radiation condition and the free surface boundary conditions respectively. It is important to note that in (3.2) the coefficients Ay and By are invariant in time. Thus, the solution to the system of simultaneous equations is required only once rather than at each time step compared to the more conventional nonlinear time-domain method in which a new system of simultaneous equations must be generated and solved at each time step. Furthermore, the coefficient matrices are functions of geometry and discretization only and consequently the same system of equations can be applied to different incident wave conditions with modifications only made to the input vector on the right-hand side. This is relatively efficient in comparison to the use of a frequency-dependent Green's function, in which a new system of equations must be generated and solved for each individual incident wave condition. Interior Solutions In some applications, the potential at an interior point is required. Once the potential and its normal derivative are known on the entire boundary, the integral equations, (2.38) and (2.35), provide a straightforward method of obtaining the potential at any interior point for the two-and three-dimensional problems respectively. Based on the same surface discretization as in the development of (3.2), the potential at an interior point x is given by 38 N <)>k(x) = £ [aj(t>^j + b j ( ^ ) j ] ( j = 1 , N ; k = 1, 2) (3.13) T h e c o e f f i c i e n t s b j a n d aj c o r r e s p o n d r e s p e c t i v e l y t o i n t e g r a l s o f t h e G r e e n ' s f u n c t i o n a n d i t s n o r m a l d e r i v a t i v e o v e r a r e a o f t h e j - t h f a c e t a s i n (3.3) a n d (3.4), a n d t h e i n t e g r a t i o n i s p e r f o r m e d b y a G a u s s i a n n u m e r i c a l i n t e g r a t i o n . I n t h r e e d i m e n s i o n s , t h e c o e f f i c i e n t s aj a n d b j a r e g i v e n r e s p e c t i v e l y b y : 1 ^ ^ r (£kmn—x) -n k aj = — I I X-/ w m w n r i (3.14) 47C k=l m=l n=l rkmn I 4 M M i bj = — I X "ZJ w m w n - — (3.15) 47C k=l m=l n=l r k m n a n d i n t w o d i m e n s i o n s , aj a n d b j a r e g i v e n r e s p e c t i v e l y b y 1 v v T (5km-*)-nk 1 „ aj = — I 1 / w m — (3.16) 2K k=i m=i r k m 1 2 M b j = - — I 1 / w m l n ( r k m ) (3.17) 27C k=l m=l w h e r e t h e s y m b o l s a r e d e f i n e d i n t h e s a m e w a y a s b e f o r e . T h e o r e t i c a l l y t h e p o t e n t i a l a t a n y i n t e r i o r p o i n t x c a n b e c a l c u l a t e d f r o m (3.13). H o w e v e r , t h e c o e f f i c i e n t o f t h e i n t e g r a l e q u a t i o n c h a n g e s i n s t a n t l y t o o n e h a l f o f i t s o r i g i n a l v a l u e w h e n x m o v e s a w a y from t h e b o u n d a r y a n d t h e r e e x i s t s a s m a l l d i s t a n c e from t h e b o u n d a r y w h e r e t h e p o t e n t i a l c a l c u l a t e d b y t h e a p p l i c a t i o n o f t h e i n t e g r a l e q u a t i o n i s i n a c c u r a t e . C o m p u t a t i o n a l l y , t h i s c o r r e s p o n d s t o t h e i n a c c u r a c y o f t h e n u m e r i c a l i n t e g r a t i o n d u e t o t h e r a p i d v a r i a t i o n s o f t h e 1/r a n d i n ( r ) t e r m s i n (3.14) - (3.17) w h e n x a n d % a r e c l o s e t o e a c h o t h e r . W h e n t h e G r e e n ' s f u n c t i o n i s a s s u m e d t o b e c o n s t a n t o v e r e a c h f a c e t ( i . e . M = 1 i n t h e n u m e r i c a l i n t e g r a t i o n ) , 39 Liggett & Liu (1983) suggested that x should be located at least one characteristic facet diameter from the boundary. However, the accuracy will certainly be improved if a higher order integration is used. 3.1.2 Alternative Second-Order Solution Unlike the Stokes second-order method, the incident potential to second order for a developed wave train is unknown for this alternative approach. The solution is then developed from an initial condition corresponding to still water in the computational domain and an incident wave train to second order which would interact with the body is generated on the upwave control surface. The simulation is continued until a sufficient duration of steady state condition is indicated. To illustrate this alternative approach and to simplify the formulation and computation, only the special case of a two-dimensional problem is considered here. First-Order Solution At first order, the total potential is separated into an incident and a scattered component as in the Stokes second-order method such that the radiation condition can be applied to each component separately. With a still water initial condition, incident waves are first developed in the absence of the body. The values of 8<t>i73n corresponding to a known incident wave train are specified on the upwave control surface and the values of ^ on the downwave control surface and still water surface are evaluated from a time-integration of the corresponding boundary conditions. The solution for the first-order incident potential can be obtained from the following system of simultaneous equations, Z A^r, + j = N | By = - | ' B i j - 1^ A y (1-1 N) (3.18) where N C 1 is the number of facets on the upwave control surface and N' is the total number of facets on Sc and S0. 40 After the first-order incident potential and its normal derivative are evaluated at a time step, the scattered components can be obtained by the imposition of the body surface boundary condition in the flow. Replacing §\ by (<J>ir+<^ ) in (3.2) and applying the body surface boundary condition (2.21), the solution to the first-order scattered potential can be obtained by the following system of simultaneous equations: (i = 1, N) (3.19) in which <t>Y on the body surface is obtained by the application of (3.13) and on S c and S 0 is evaluated from the corresponding boundary conditions. The first-order solution which is free from reflection can be obtained by the superposition of the incident and scattered potentials evaluated from (3.18) and (3.19) respectively. Second-Order Solution With the same initial condition as in the first-order problem, the total potential at second order is developed on the basis of the first-order solution. Based on.the same numerical discretization as in the Stokes second-order method, the second-order integral equation (2.50) is reduced to the following system of equations: I [Ay <)>2j + By ( A ] = - X Tin [Cij Cg^ -). + Dy OJ j=i J j=i J - j J J > i [ B i J ( a ^ ) i + E l i ' , ' » ] + X ? * [ B B < * ) J + I * * U ] °- 1'-- N ) ( 3- 2 0 > 41 where N ' is the number of facets on Sb and S c , and the right-hand side input vectors are known from the first-order solution. The coefficients Cy, Dy, Ejj and Fy are evaluated numerically in the same manner as Ay and By in (3.9) and (3.10) respectively: 1 3 ?^ r z—Ckm - I 1/ w m — — 7t k=l m=l r k m 1 + - I / w m — — ^ m=l r2m (i * j) ( i = j ) (3.21) I v V f w T - ! ^ 2nk-(^km-x)(Ckm-z)1 = < 7C k=l m=l "Tkm2 r k m 4 1 M - I / W m i 2 ^ 7t m=l r2m C2 2n2-(^2m-x)(C2m-z) r2mq ] Eii = « 1 v v , r n& 2nk-(|k m-x)(c;km-z)1 - I 1/ w mn ? kLr2 - —4 J TC k=l m=l r km r k m 1 V T „ „ T i ^ 2 - 2n2-($2m-x)(C2m-z)-| - 1/ wmn ; 2Lrr2 7 ^ J '2m ( i * j ) ( i = j ) (i * j) (i = j) (3.22) (3.23) Fi; = < ( 1 2 M x _ £ . - I I / W m 2 TC k=l m=l r k m L 0 ( i * j ) (i = j) (3.24) where the symbols are defined as consistent with those of (3.9) and (3.10). With all terms on the right-hand side known from the first-order solution, (3.20) represents a total of N equations with 2N variables corresponding to the values of § 2 ^ d dfa/dn at each of the facets. With 8<j>2/3n = 0 applied on Sb and (J>2 on S 0 evaluated from the free surface boundary conditions, (3.20) is reduced to a system of simultaneous equations if either 9<j)2/3n or <t>2 is known on S c . However, the radiation condition associated with the generation of the 42 incident waves at second order is not properly defined and therefore 9<t>2/3n = 0 is arbitrarily applied on the control surface to simulate a rigid-wall condition at second order. In this approach, the incident waves at second order are generated by the quadratic terms in the free surface boundary conditions at second order and thus are phase-locked with the first-order wave system. Similar to the diffraction problem at second order, these forced waves interact with the control surface and give rise to free waves at twice the incident wave frequency. These free waves are synchronized with the generation of the second-order forced waves and propagate into the computational domain at the corresponding group velocity. Furthermore, the expansion in (2.45) which provides the basis for the development of (2.50) gives bad approximations near the intersections between the still water surface and the control surface. Numerical errors in the form of additional free waves at twice the incident wave frequency are also observed to be generated near the control surface. Since the amplitude and group velocity of these free waves are generally small compared to their first-order counterparts, they will influence the flow near the body at a later time. Therefore, a relatively large domain is used in conjunction with the first-order radiation condition to ensure that a sufficient duration of steady state solution can be obtained before these spurious free waves arrive at the test section. 3.2 Temporal Solution 3.2.1 Radiation Condition In order to be able to simulate a sufficiently long duration in a reasonably sized domain, the application of a radiation condition on the control surface requires particular consideration. In the Stokes second-order method in § 2.2, the Sommerfeld radiation condition with a time-dependent celerity is integrated numerically in time to obtain the required scattered potentials at first and second order on the control surface. 43 Time-Stepping Equation When the control surface is located at a sufficiendy far distance from the body, the Sommerfeld radiation condition (2.31) for the scattered potential at k-th order can be rewritten in term of the normal derivative of the potential and applied on the control surface as indicated in (2.32). In order to develop a numerical solution, (2.32) is expanded by a central difference in time and a leap-frog difference in space (Orlanski, 1976): <fi(x,t+At)- <J>k(x,t-At) + 2At (T L ){^[<^(x, t+At) + (l)k(x,t-At)]-(t)k(x-nAn,t)} = 0 ( k = l , 2 ) (3.25) An * where x is a point on the control surface; At is a time step size; An is a small distance of the order of a characteristic facet diameter; and n is the normal vector at x. From (3.25), the potential on the control surface at an advanced time (t+At) can be expressed in terms of the solutions obtained at previous time-steps: <&(x,t+At) = T=| 44(x,t-At) + <&(x-nAn,t) (k = 1, 2) (3.26) 1+P 1+p in which 3 = ( £ ) c „ (3.27) where |3 is the Courant number on the control surface which should be maintained at less than one for stability reasons. Literally, the boundary value of fa is extrapolated in time and space from the values of fa near the boundary inside the domain at previous time-steps through the knowledge of the celerity c n which is numerically evaluated at each facet on the control surface at time t. 44 Evaluation of Celerity For diffraction of regular waves at first order, it is convenient to introduce an approximation that the scattered waves at a large but finite distance from the body propagate in the radial direction measured from the body centroid and that their speed c is known on the basis of the linear dispersion relation. In three dimensions, the control surface can take on any orientation with respect to the radial direction. In applying (3.26) and (3.27), c n need not be evaluated numerically at each time step, but instead be obtained as c/cosy, where y is the known angle between the normal and radial directions at the control surface. This assumption is generally valid when a large single body in waves is considered and when the control surface is located sufficiently far from the body. When diffraction of irregular waves is considered or when more than one diffracting member are present, this assumption can no longer apply and a numerically evaluated celerity on the control surface should be used. For diffraction of regular waves at second order, the scattered potential is made up of second-order free waves at twice the incident wave frequency as well as second-order forced waves associated with the inhomogeneous terms in the free surface boundary conditions at second order. Therefore, the radiation condition cannot be satisfied with a unique celerity. Instead, the celerity on the control surface at a time t is approximated by its value at the previous time step (t-At) at locations (x-nAn) which are just within the control surface. This is done on the basis of: c n = (3 28) la /^anj ( 3- 2 8 ) This is consistent with the flow of information from within the fluid domain to the boundary. From (3.28), the numerical value of c n at (t-At) and (x-nAn) is given by <|)k(x-nAn,t) - <|)k(x-nAn,t-2At) .(^(x-nAn.t) + (t>k(x-nAn,t-2At) - 2<))k(x-2nAn,t-At). ( k = l , 2 ) (3.29) 45 Unfortunately, even at first order there is a numerical difficulty in using (3.29) in that both d$\/dt and dfa/dn are zero at locations where TJi = 0, and (3.29) becomes undefined. Even at locations near zero-crossing points the numerical values of c n calculated by (3.29) are inaccurate (e.g. Lee & Leonard, 1987). At second order, the numerical procedure is further complicated by the presence of both forced and free waves, in that the value of d<t>2/9t is not necessarily zero when dfa/dn = 0, so that the resulting celerity calculated by (3.29) may be unreasonably large. In order to circumvent these numerical difficulties, the following numerical scheme is carried out. When the absolute values of both 3<j)k/9t and 3()>k/9n are less than a certain prescribed value related to the chosen computational precision, the scattered potential at the advanced time <|)k(x,t+At) in (3.26) is simply given by <)>k(x,t). This is equivalent to taking c n = 0, and is used because fa changes very little in one time step on account of dfa/dt being small. As the Sommerfeld radiation condition requires the scattered waves at large distance from the body to be outgoing, only the outward celerity (i.e. c n ^ 0) is used in describing the boundary condition. Finally, the Courant number on the control surface given by (3.27) is maintained at less than one for stability reasons. The maximum value of the celerity calculated by (3.29) is then limited to have a maximum value of An/At. Thus, the numerical values of c n are given overall by: c n = < 0 for c n < 0, or when 8<|)k/9t, 9<))k/9n are small c n for 0 < c n < An/At (3.30) An/At for c n > An/At where c n is the initially calculated value of the celerity from (3.29). Although this numerical radiation condition was originally developed for hyperbolic flows, its applications in parabolic problems have been shown to be effective by Han et al. (1983). 46 3.2.2 Free Surface Boundary Conditions The free surface boundary conditions are used in a time-stepping procedure to provide the free surface elevation and velocity potential on the still water surface at each time step. Various time-stepping equations are available, and as a simple and reliable procedure the first-order Adam-Bashforth-Moulton method is adopted here. The method also has the advantage of providing an optional iteration scheme to account more accurately for the rapid variation of the solution in time. Initially, the free surface elevation at a new time (t+At) is first evaluated in terms of the known solution up to time t by the first-order Adam-Bashforth equation as T f e t M O - l f e O + f t s ^ ) , - ^ ) ^ (k = l , 2 ) (3.31) The dynamic free surface boundary condition, (2.23) or (2.29), can now be used to provide dfa/dt at (t+At). In turn the potential fa may then be obtained by the first-order Adam-Moulton equation: $ s k(. +A.) = <(«)4[(^)t+(^)t+J (k = l , 2 ) (3.32) With the normal derivative of the potential on Sb known from the incident wave potential and the potential on S c and S Q evaluated from the time-stepping equations, the flow in the domain at the advanced time (t+At) can now be solved from the discretized boundary integral equation (3.2). The output vector provides 3(()k/3n on S 0 at (t+At). The kinematic free surface boundary condition, (2.22) or (2.28), can now be used to provide the corresponding values of dr\ydt at the advanced time (t+At). This completes an initial calculation of all the boundary parameters at time (t+At) and the computation could then proceed to the next time step. However, an improvement to these values may first be obtained by applying an iterative procedure in which improved values of the 47 free surface elevation at (t+At) are obtained by using the first-order Adam-Moulton equation in place of (3.31): r, k(t +At) = nsk(t) + f [(^)t + (^) t +J (k = 1. 2) (3.33) The improved values of T| k can then be used to obtain in turn dfa/dt from (2.23) or (2.29), fa from (3.32), dfa/dn from (3.2) and dv^/dt from (2.22) or (2.28). This iterative procedure is thereby repeated until a sufficient convergence is obtained. In the present application, the time-stepping procedure converges rapidly and normally the iterative procedure is not required. The iterative procedure involving (3.31), (3.32) and (3.33) is thus applied for the first few time steps to account for the transient effects associated with the initial conditions, and the subsequent development of the flow is obtained by the applications of (3.31) and (3.32) without the use of iteration. 48 CHAPTER 4 RESULTS AND DISCUSSION 4.1 Introduction Based on the preceding formulation, two sets of F O R T R A N programs have been developed respectively for the two- and three-dimensional problems. Particular attention is paid to the performance issues of the computer programs running on an I B M 3090/150S with vector facility and virtual memory. In the development of the computer codes, a set of guide-lines published by the Palo Alto Scientific Center of I B M Corporation (Dubrulle et al., 1985) have been followed. In order to demonstrate and verify the present approach, the Stokes second-order method discussed in § 2.2 is applied to examine three classical examples in which numerous theoretical and experimental results have been accumulated in the literature. These are (1) a fully submerged circular cyUnder in two dimensions; (2) a semi-submerged circular cylinder with axis at the still water level in two dimensions; and (3) a bottom-mounted surface-piercing circular cylinder in three dimensions. Only the diffraction of regular waves around rigid bodies is considered here and the extension to irregular wave diffraction can readily be made on the same basis. A number of computational considerations inherent with the formulation is first considered. The proposed method is validated by convergence tests of the numerical solution with respect to the time-step and facet sizes. Results of runup and wave forces are presented and comparisons are made with previous theoretical and experimental studies. In addition, the variation with time of the three-dimensional free surface profiles to first and second order has been transformed to a colour film as an aid to viewing the development of the flow. 49 The alternative second-order method discussed in § 2.3 is also applied to the first example involving a fully submerged circular cylinder in order to examine the effects of the expansion terms in the integral equation on the calculation of second-order wave forces. Wave forces obtained by the two methods are compared and the applicability of the alternative expansion procedure is discussed. 4.2 Computational Considerations 4.2.1 Discretization Schemes Two-Dimensional Problem The two-dimensional numerical model corresponds to that of a narrow wave flume at which sectional models are tested. The body surface, the control surface and the still water surface are discretized into numbers of facets made up of straight-line segments. A F O R T R A N program GEOM2D has been developed to perform the surface discretization and to generate and solve the system of simultaneous equations. About 200 to 300 facets are generally used for the discretization of the surfaces. The coefficients Ay and By in (3.9) and (3.10) are evaluated by a 4-point Gaussian integration and double precision is used throughout the computation. The program is capable of generating circular cylinders of various depths of submergence except for the case when the surface of a fully submerged cylinder is in contact with the still water surface for which the configuration gives rise to a singularity in the coefficient matrices. Figure 2 shows two examples of the numerical model corresponding to a fully submerged and a semi-submerged circular cylinder. In the figure, Nb, N c and N 0 represent respectively the numbers of facets on the body surface, the control surface and the still water surface, while N denotes the total number of facets in the discretization. It is noted that the seabed is excluded from the procedure due to the choice of the Green's function that accounts for the symmetry 50 about the seabed. The facet sizes on the body surface and on the nearby still water surface are smaller than those in the far field and are specified independent of the facet size in the far field. This layout scheme is designed to describe the more rapid variations of the potential near the body and to speed up the convergence on wave force and runup calculations without increasing facet density everywhere. It also provides a more effective treatment for cases when the cylinder size is small or when the clearance between the cylinder and the still water surface is small. A second F O R T R A N program WAVE2D has been developed to time-step the development of the flow and to calculate runup and wave forces at each time step. Two initial conditions can be simulated. One corresponds to still water everywhere in the computational domain and an incident wave train about to be generated on the upwave control surface, while the other corresponds to a fully developed second-order incident wave train in the domain and the scattered potential is developed in time and space through the imposition of the body surface boundary condition. The simulation is continued until a sufficient duration of steady state solution has been obtained Three-Dimensional Problem The three-dimensional numerical model corresponds to that of a rectangular wave basin with a surface-piercing circular cylinder located at the centre. A F O R T R A N program G E 0 M 3 D has been developed to discretize the body surface, the control surface and the still water surface into finite numbers of planar quadrilateral facets, and to generate and solve the system of simultaneous equations. Depending on the geometry and incident wave frequencies of interest, the number of facets used in the discretization may vary from 1500 to 2200. The coefficients Ay and By in (3.3) and (3.4) are evaluated by a 16-point Gaussian integration and double precision is used throughout the computation. Figure 3 shows examples of the surface discretization scheme for the three-dimensional problem. Both examples represent the same geometry but with different numbers of facets on 51 the body surface and on the nearby still water surface. The seabed and one half of the total surface are excluded due to the choice of the Green's function that accounts for double symmetry. Smaller facets are used on the body surface and on the nearby still water surface in order to describe the more rapid variations of the potential near the body and to speed up the convergence on wave force and runup calculations. The number of facets in the z direction is deterrnined such that the facets nearest to the still water surface are square and the facets nearest to the seabed have a specified aspect ratio. For the examples shown in figure 3 the aspect ratio is specified at 1.5. A F O R T R A N program W A V E 3 D has been developed to time-step the development of the flow and to calculate wave forces and runup at each time step. The initial conditions correspond to a Stokes second-order wave field in the domain and the scattered potential is developed in time and space through the imposition of the body surface boundary condition. 4.2.2 Computing Cost When the rank of the matrix equation is large, the computing cost might pose an important issue affecting the applicability of the method. With the rank of the matrix equation denoted by N , the number of floating-point operations for the first part of the geometry programs GEOM2D and GEOM3D which generates the matrix coefficients is proportional to N 2 , while the second part of the programs which solves the matrix equation is proportional to N 3 . Likewise, the number of floating-point operations for the time-stepping programs W A V E 2 D and W A V E 3 D is proportional to N 2 and the number of time steps. Since the computer codes are executed on the I B M 3090/150S computer with vector facility, it is expected that for large values of N the execution wil l be more efficient due to vectorization. On the other hand, the computing cost also depends on a number of factors such as the structure of the computer program, the storage requirement and the availability of efficient library subroutines. The C P U time for a specific job may vary slighdy from time to 52 time depending on the computing environment and the results presented in this section can only be used for reference purposes. Figures 4 and 5 show the CPU time as a function of the rank of the matrix equation for the two- and three-dimensional problems respectively. It is noted that both the CPU time and its slope increase gradually with the rank of the matrix equation. The results also indicate a faster rate of increase in both the CPU time and its slope for the geometry programs than for the time-stepping programs due to their respective relations with N . In the two-dimensional problem, a 4-point integration is used in the evaluation of the matrix coefficients and each coefficient itself involves two facets in the domain and in the image about the seabed, while in the three-dimensional problem a 16-point integration is used and each coefficient involves four facets in the four quadrants of the doubly symmetric geometry. Therefore, for the same value of N G E O M 3 D is computationally more intensive and requires more C P U time compared to GEOM2D. Although GEOM3D is computationally intensive, the matrix equation is a function of geometry and discretization only. The results generated by the geometry programs can be stored in the disk memory and used subsequently in the time-stepping programs for different incident wave conditions. The C P U time for the time-stepping programs is linearly proportional to the number of time steps, and usually about 200 to 300 time steps are required to generate a sufficient duration of the steady state solution. Therefore, the application of the present method with a large number of facets is not economically prohibitive but then the accuracy and stability of the matrix solution require particular attention. 4.2.3 Accuracy and Stability of Matrix Solution Condition Number When the rank of a system of simultaneous equations is large, the round-off error in the matrix solution requires particular attention before it can be applied more generally to the time-53 stepping procedure. A number of different techniques can be used to indicate the magnitude of the round-off error in a matrix solution. One of the most frequently used is the condition number proposed by Turing (1948) which has the property that the higher the condition number, the less accurate is the matrix solution. For orthonormal matrices, the condition number is one which indicates no round-off error in the matrix solution. On the other hand, square matrices with elements chosen at random have a condition number approximately equal to the rank of the matrix. In principle, only when the condition number is significandy greater than the rank of a matrix can the matrix or its solution be regarded as ill-conditioned. However, the order of the condition number is also proportional to the loss of accuracy in the solution process. Therefore, given a computational precision there might exist a limit on the size of a matrix equation to which a valid solution can be obtained For the two-dimensional cases of a submerged and a semi-submerged circular cylinder, the condition number as a function of the rank of the matrix equation is shown in figure 6. In the examples, the facet layouts on the cylinders and on the nearby still water surface are independent of the rank of the matrix equation and the increase in the total number of facets is accomplished by increasing the length of the domain. It is noted that the two curves in the figure are almost identical, and variations of the model geometry and discretization scheme while keeping the total number of facets constant do not affect the condition number significandy. Since the condition number is considerably smaller than the rank of the matrix equation, the matrix solution for the two-dimensional problems is well-conditioned. For the range of matrix sizes of interest and with the use of a double-precision computation, the round-off error in the matrix solution is relatively insignificant Since the rank of the matrix equation is much higher in the three-dimensional problem, the condition number should be investigated more rigourously before the matrix solution can be applied more generally. For the case of a surface-piercing circular cylinder, figure 7 shows the condition number as a function of the rank of the matrix equation. It is observed that the 54 condition number is much smaller in relation to the rank of the matrix equation so that the matrix solution is very accurate. For the rank of the matrix equation as large as 2300, the condition number is less than 4.1 and the loss in accuracy in the matrix solution is insignificant. A n inspection of the coefficients of the matrix equation as given by (3.3) and (3.4) reveals that Ay and By are made up of terms composed of 1/r2 and 1/r respectively, where r is defined as the distant between the facets i and j . The magnitude of the coefficients A y and By is then largest along the matrix diagonal and decays rapidly from the diagonal except at locations where facets i and j are physically close to each other in the discretization scheme. Since the matrix is dominated by elements along and near the diagonal, the round-off error in the solution process is generally small. Therefore, the solution to a very large system of simultaneous equations can confidently be obtained for the three-dimensional problem without inducing much round-off errors in the matrix solution. With the Cray X-MP/48 at the Pittsburgh Supercomputing Center, Korsmeyer et al. (1988) have successfully obtained converging solutions for the first-order diffraction and radiation problems of a tension-leg platform, based on the use of up to 12,608 facets to represent the body surface. Since a large complex matrix was involved in their method, out-of-core storage was required and an accelerated iterative solver of the Gauss-Seidel type was adopted. In the present application, the I B M 3090/150S has a virtual memory of 768 Mb and can, in principle, accommodate a real matrix with double precision of rank up to 9800 without the use of out-of-core storage. However, the performance of the computer is expected to start deteriorating when the virtual memory of the program exceeds about 250 Mb or when out-of-core storage is required. Irregular Frequencies of Surface-Piercing Bodies When a frequency-dependent Green's function is used, the matrix equation obtained from the discretization of an integral equation is a function of the incident wave frequency. When surface-piercing bodies are involved, there exist a series of irregular frequencies at which the 55 matrix equation is singular and the solution is not unique. Even at locations near the irregular frequencies the numerical solution of the matrix equation is inaccurate. Since the irregular frequencies are usually high in relation to the frequency range of interest, their effects on the application of the first-order solution are generally small. At second order, when wave interactions occur at twice the incident wave frequency, numerical problems associated with irregular frequencies are more critical (e.g. Miao & Liu, 1989; and Kim & Yue, 1989). Nevertheless, the Green's functions used in the present application do not depend on incident wave frequency and problems associated with irregular frequencies have been proved not to arise in this case (Angell et al., 1986). The matrix solution generated for a given geometry can be applied to different incident wave frequencies provided the spatial resolution of the wave fields at first and second order is sufficient. 4.2.4 Initial Conditions For the Stokes second-order method, the initial conditions are known and correspond to a Stokes second-order wave field in the domain. The development of the scattered potential in time and space is attained by the imposition of the body surface boundary condition. In the simulation of full nonlinear wave diffraction with similar initial conditions, Stansby & Slaouti (1984) and Isaacson & Zuo (1989) applied the body surface boundary condition instantly at the first time step corresponding to an impulsive-started condition. The transient effects associated with the abrupt initial conditions were found to decay rapidly in time and the established flow is reached after relatively few time steps. On the other hand, by modulating the body surface boundary condition in time it is possible to avoid such an abrupt initial condition and instead allow a gradual development of the scattered potential. Before applying in (3.2), the right-hand sides of (2.21) and (2.27), which correspond respectively to the imposition of the body surface boundary conditions at first and second order are multiplied by the following modulation function, 56 [ l - cos( ;?-) ] for t < T m 2 A m (4.1) I 1 for t > T m where T m is a modulation time. The modulation function is zero at t = 0 and increases gradually to one at t = T m . Consequendy, a smooth transition between the fully developed initial conditions and the subsequendy developed scattered potential is obtained by the gradual imposition of the body surface boundary condition in the computation. In order to illustrate the influence of T m , the two-dimensional case of a circular cylinder with axis at the still water level is considered. Figure 8 shows comparisons of free surface profiles obtained at time t = 5T using different modulation times, T m = 0, T/2 and T. In the figure, A = H/2 is the amplitude of the incident waves, a is the cylinder radius and the cylinder centre is located at x/L = 2. When the modulation time T m = 0, some high frequency dispersive components associated with the impulsive-started condition are just noticeable in the flow at first order. However, their effects are considerably amplified at second order. On the other hand, the same phenomenon is not observed in the profiles corresponding to T m = T/2 and T. Since the amplitudes of these high frequency components are generally small, their effects are mainly confined to the region near the free surface. The influence on the resultant force acting on the body is usually small except for relatively small bodies. However, the runup is affected to a greater extent. Through extensive numerical testing, it has been found that a smooth transition between the initial conditions and the fully developed flow is essential to ensure a stable and steady solution at second order. For ka less than about 0.3, it is difficult to obtain a steady solution at second order without the use of a non-zero modulation time. The impulsive-started condition also affects the performance of the radiation conditions and eventually the duration of simulation. Numerical testing of the effects of T m on the predictions of wave forces and runup to second order has also been carried out for the cases of a fully submerged circular cylinder in two dimensions and a surface-piercing circular cylinder in three dimensions, and this indicates 57 a less noticeable effect of the modulation time on the solution. Nevertheless, in the results subsequently reported in this thesis, a modulation time T m = T is adopted and generally a steady state solution in the vicinity of the body is obtained immediately after the duration of the modulation. Alternatively, it might be possible to adopted a known analytical solution at first order and solve only the solution at second order through the time-stepping procedure. To avoid an abrupt initial condition and allow a gradual development of the solution in time and space, both the right-hand side of the body surface boundary condition and the quadratic forcing terms in the free surface boundary conditions can be modulated by (4.1), and similarly, a steady state solution at second order would be obtained in the vicinity of the body after the modulation. This approach has the advantages of reducing the computing time and providing an accurate evaluation of the quadratic forcing terms in the free surface boundary conditions. However, these known analytical solutions at first order are usually restricted for bodies of simple geometry. Although the present method is illustrated through examples with known analytical solution at first order, it is intended for applications with structures of arbitrary shape. Therefore, the more general approach of solving both the solutions at first and second order numerically by the time-stepping procedure are adopted here. 4.2.5 Stability of Time-Stepping Procedure As with other problems which incorporate time stepping, the stability of the proposed solution requires that the time-step size At is sufficiently small in relation to a characteristic facet diameter AS. This condition corresponds to a minimum value of T/At, for any given value of L/AS, where T and L denote a representative wave period and length. To illustrate the relation between (T/At)nun and L/AS, a two-dimensional numerical model corresponding to a semi-submerged circular cylinder in regular waves is examined, and a constant facet size is used on 58 the entire boundary while in practice smaller facets are used on the body surface and on the nearby still water surface. For incident wave conditions correspond to ka = 1 and deep water, table 1 shows the minimum values of T/At for various values of L/AS which give rise to stable solutions at first and second order. For smaller values of T/At, the solution exhibits an oscillatory instability, while for larger values the solution remains quite stable. In the present application, the time-stepping procedure converges rapidly and the solution appears not to be affected by increasing the value of T/At above the minimum required or by using the iterative approach involving (3.33). However, when the time-stepping procedure with one iteration is used, the minimum value of T/At required to maintain stable solutions at first and second order is found to be less than half that when iteration is not used. The robustness and effectiveness of the present time-stepping procedure can be demonstrated through the observation that a stable second-order solution can be obtained for L/AS and T/At as small as 10 and 7 respectively. The Courant stability condition which specifies that the Courant number given by cAt/AS = (L/AS)/(T/At) is less than one appears to apply to the cases when L/AS is approximately less that 30 and when the iterative approach is not used. For larger values of L/AS or when the iterative approach is applied, a stable solution can be obtained for a time-step size which is significandy larger than that required by the Courant stability condition. This indicates that when a very fine mesh is required in regions where the geometry is complicated or the potential changes rapidly in space, it is not necessary to increase the number of time steps per period proportionally. In the subsequent applications, smaller facets are used on the body surface and on the nearby still water surface in order to render more accurate calculations of the wave force and runup without increasing the number of facets on the entire boundary. The smallest facet in the discretization is then used with table 1 to estimate the required time-step sizes to maintain stable solutions at first and second order. 59 4.2.6 Numerical Accuracy of Discretization Quite apart from the question of stability, the numerical accuracy of the solution must also be assessed in terms of the facet size chosen. The computed runup and wave forces are found to depend on the number of facets used to represent the body surface, and in general, the convergence of the computed runup and wave forces with increasing number of facets depends on a number of factors such as the relative size and water depth. To illustrate these effects and to ensure the accuracy of the numerical solutions, the convergence of the free surface elevations and wave forces at first and second order is examined for each of the three numerical examples. Submerged Circular Cylinder in Two Dimensions For the case of regular incident waves interacting with a submerged circular cylinder in deep water, table 2 shows the convergence of the various wave force components for different values of Nb, the number of facets on the cylinder surface. These results correspond to the conditions ka = 0.4, h/a = 1.75 and 2.00 which are quite typical in the frequency range of interest. Here h is the depth of cylinder axis below the still water level. Each of the first- and second-order oscillatory forces (denoted by F i and F2 respectively) refers both to the x and z components, whose amplitudes are equal for deep water, while the second-order steady force (denoted by F 0 ) refers to the component in the positive z direction only (the x component being zero). Since the cylinder is fully submerged, the contribution from the free surface fluctuation term toward the second-order oscillatory force is zero. For the deep water case considered here, the contribution from the velocity-squared term towards the second-order oscillatory force is zero, and F2 can be evaluated on the basis of fyi only. The first-order force component is found to converge rapidly with the number of facets representing the cylinder and very accurate results can be obtained with Nb as low as 20. Since the second-order steady force is a function of the first-order potential as well as its gradient on the cylinder surface, the convergence should follow that of the first-order force 60 component but with more dependence on the facet size. The convergence of the second-order oscillatory force is relatively rapid except when the depth of submergence is small and when the facet size is large. This is due to the poor resolution of the potentials at first and second order in the region between the cylinder and the still water surface. Since the potential at second order is a function of a number of variables, the convergence of the second-order oscillatory force with increasing Nb is not monotonic. A reduction of time-step size below the maximum value required for numerical stability was found to have no significant effect on the computed force components. Semi-Submerged Circular Cylinder in Two Dimensions For a semi-submerged circular cylinder in two dimensions, the amplitudes of the free surface elevations and wave forces at first and second order are shown in tables 3 to 5 as function of Nb- The results correspond to the case of ka = 0.6 and deep water. The amplitudes of the free surface elevations at first and second order on the cylinder in their respective dimensionless forms are shown in table 3 as functions of Nb. The first-order amplitude Irijl is found to converge rapidly with the number of facets and an accurate result can be obtained with Nb less than 15. Since r | 0 can be calculated from the first-order wave field, its convergence is slower due to accumulation of numerical errors, but generally follows that of the first-order solution. In the deep water case considered here, the wavelengths of the second-order free and forced waves are respectively one-quarter and one-half of the first-order wavelength, and more facets per first-order wavelength are required to obtain a satisfactory resolution of the wave profile at second order. Therefore, the convergence of 1%! is much slower compared to other components and a smaller facet size is required to achieve satisfactory results. Due to the dominance of second-order forced waves, the wave field in the upwave region is more dependent on the first-order solution and the convergence of Ir^l on the upwave side is slower than that on the downwave side. 61 The amplitudes of the wave force components at first and second order in their dimensionless forms are shown in tables 4 and 5. The first-order force components in the x and z directions (denoted by F i x and F i z respectively) are found to converge rapidly with the number of facets representing the cylinder and very accurate results can be obtained with Nb as low as 10. The convergence of the second-order steady force components (denoted by Fo x and Fte) and the second-order oscillatory force component due to <|>i (denoted by F ^ and F2z) is more rapid for the z direction than for the x direction. This is because the force components in the x direction contain an additional contribution from the free surface fluctuation as well as the velocity-squared term in the Bernoulli equation and hence more numerical errors are (2) (2) involved. The convergence of the second-order oscillatory force components, F2X and F 2 z , due to <j>2, is slower compared to that of F2I and F2Z, due to <|>i, while numerical errors are carried over to the resultant second-order oscillatory force components, denoted by F2X and F 2 z . Bottom-mounted Surface-Piercing Circular Cylinder in Three Dimensions For a bottom-mounted surface-piercing circular cylinder in three dimensions, only the force components in the x direction and the moment components about the y axis are non-zero and are considered here. The convergence of the free surface elevations, wave forces and moments at first and second order with respect to the number facets representing half of the cylinder surface is shown in tables 6 to 8. The results correspond to the case of ka = 1 and kd = 1. The amplitudes of the free surface elevations at first and second order on the cylinder are shown in their respective dimensionless forms in table 6 as functions of Nb- In the table, 8 is the azimuthal angle measuring horizontally from the x axis. The first-order amplitude Ir^l has a maximum at 8 = n and is found to converge rapidly with the number of facets. The steady component TJ 0 , which can be calculated from the first-order wave field, also has a maximum at 8 = TC, but its convergence is slower due to accumulation of numerical errors. Since the downwave region of the second-order wave field is dominated by free waves which propagate 62 independently of the first-order waves, the convergence of Iry at 8 = 0 is relatively rapid and is similar to that of the first-order solution. Due to the dominance of second-order forced waves, the wave field in the upwave region is more dependent on the first-order solution and the convergence of Ir^l at 8 = TI is slower than that at 6 = 0. The maximum Ir^l occurs along the side of the cylinder and both the amplitude and the location of the maximum are found to depend on the number of facets representing the body. The wave force components in the x direction and the moment components about the y axis at first and second order are shown in tables 7 and 8 respectively. The first-order force and moment, F i and M i , are found to converge rapidly with the number of facets representing the cylinder, and very accurate results can be obtained with Nb as low as 36. In spite of the accumulation of numerical errors at second order, the convergence of the second-order steady force components Fo and Mo, and the second-order oscillatory force components F ^ and M ^ , both due to $i is relatively rapid. The convergence of the second-order oscillatory force (2) (2) components F 2 and M 2 due to 02 is generally rapid except for very coarse facet sizes. On the other hand, the resultant second-order oscillatory force components F2 and M2 are fairly constant over the range of Nb considered, which might be on account of the cancellation of numerical errors. 4.2.7 Effectiveness of Radiation Condition Application to First-Order Bichromatic Waves Since the radiation condition is satisfied by a time-integration involving a time-dependent celerity, its performance depends on the effectiveness of the numerical scheme in which the celerity is evaluated. Prior to exarnining the application of the present radiation condition to the second-order diffraction problem, the application of this to a linear bichromatic wave train which is known to contain components at two distinct celerities is examined. The flow 63 p o t e n t i a l f o r a l i n e a r D i c h r o m a t i c w a v e t r a i n c o m p o s e d of t w o f u n d a m e n t a l f r e q u e n c i e s m a y be expressed as: where the subscripts 1 and 2 indicate the two harmonic components of the combined wave train. The celerity of the bichromatic incident wave train at x and t can be calculated from the Sommerfeld radiation condition as It is noted that the celerity calculated by (4.3) is no longer a constant in either time or space, and its value is undefined when the denominator is zero. In figures 9 and 10, the variations with time of the celerity and potential at x = 0 calculated respectively by (4.3) and (4.2) are compared with those calculated by the numerical method adopted here, based on (3.30) and (3.26) respectively. Figure 9 corresponds to the case of ai/a 2 = 3, C i / c 2 = 0.5 and deep water and shows two cycles of the biharmonics. From the figure, the numerical values of the celerity calculated from (3.30) is shown to be inaccurate near instants at which dfy/dt = 0. Since both 3<))/8t and dfy/dx are small at these instants, the potential calculated by (3.26) based on the numerical values of the celerity is affected to a lesser extent, and identical agreement with the potential given by (4.2) is observed. It is also noted that the celerity given by (4.3) is fairly constant between the troughs of the potential and is continuous across the crest where an axis of symmetry can be located. This indicates the resemblance between this portion of the biharmonics and a simple harmonic wave train with a constant celerity. Figure 10 corresponds to the case of ai/a 2 = 2, c i /c 2 = 0.8 and deep water and shows half a cycle of the beat phenomenon. Due to the decaying amplitude of the potential, the variation of the celerity in time is less steady compared to that in figure 9. The celerity given by (4.3) is <]>(x,t) = ai cos [ki(x-cif)] + a 2 cos [k2(x-C201 (4.2) c(x,t) = aikjCi sin [ki(x-cit)] + a 2 k 2 c 2 sin [k 2 (x-c 2 t)] a i k i sin [ki(x-cit)] + a 2 k 2 s i n [k 2 (x-c 2 t)] (4.3) 64 discontinuous at the troughs and crests of the potential except at t/T2 = 0 and 2, at which an axis of symmetry can be located. Despite the slight deviation of the celerity calculated by (3.30), the potential calculated by (3.26) is found once more to be identical to that given by (4.2). These two examples indicate that even for a bichromatic wave train which is known to be made up of two harmonic components and propagates unsteadily in time and space, can be accounted for accurately by the present radiation condition. The application of a time-dependent celerity in the Sommerfeld radiation condition to account for the propagation of unsteady waves appears to be valid and sufficient Application to Regular Wave Diffraction to Second Order When second-order diffraction in three dimensions is considered, both the amplitudes and nonlinearities of the scattered waves at second order decrease with the distance away from the body and the effects of the radiation condition are less critical. Although the propagating characteristics of the scattered waves are simpler in the two-dimensional problem, the amplitudes and nonlinearities of the scattered waves persist to the control surface and pose a rigourous test for the radiation condition. To examine the effectiveness of the radiation condition in the present application, a semi-submerged circular cylinder in two dimensions is considered. This is preferable to the case of a fully submerged cylinder in that more waves are reflected by the body and subsequently partial standing waves of non-decaying amplitude are developed in the upwave region. These first-order partial standing waves give rise to a strong second-order forcing and have traditionally created a challenging problem for the second-order radiation condition in two dimensions. Together with the second-order free waves generated by the homogeneous terms, a rigourous test for the radiation condition at second order is posed. In order to examine the effectiveness of the radiation condition at first and second order, the development with time of the oscillatory force components due to §\ and <j>2 predicted by two 65 models of different lengths, Lrj/L = 2 and 4 is compared in figures 11a and l i b respectively. In both cases the cylinder axis is located at the still water level and half way along the domain's length, and the cylinders are subjected to identical incident wave conditions corresponding to ka = 0.6 and deep water. The figures show results for the smaller domain with a rigid wall condition applied at the radiation boundaries, and with both a constant celerity and a time-dependent celerity used in the radiation condition. These are compared to results for the longer domain for which reflected waves are absent. The development with time of the first-order oscillatory force components corresponding to the various treatments at the radiation boundaries is shown in figure 11a. For the present case of regular wave diffraction, the solutions obtained with a constant celerity and with a time-dependent celerity both agree well with the solution obtained with the large domain at which reflected waves are absent. When the radiation condition is not applied, the results indicate significant reflected waves as expected. Since the initial waves in the scattered wave train are unsteady with elongated lengths and thus travel at higher speeds than do the subsequent waves, the solution obtained without the application of the radiation condition diverges from the other solutions much earlier than that predicted on the basis of the first-order group velocity. Without the radiation condition at first order, the reflection of these elongated initial waves would affect the flow near the body well before a steady state solution has developed. When regular wave diffraction at first order is considered, the use of the constant celerity obtained from the linear dispersion relation appears to be sufficient, and the approach is adopted for the subsequent computations. With the first-order radiation condition satisfied by applying the constant celerity of the first-order incident waves, figure l i b shows the development of the second-order oscillatory force components due to fa subject to the various treatments at the radiation boundaries. It is evident from the figure that the solution obtained using the time-dependent celerity gives the best agreement with the solution obtained with the large domain. Although a slight effect of reflection is observed in the wave force prediction, the use of the time-dependent celerity is 66 expected to give even better results in three-dimensional problems in which the scattered wave amplitude decays with distance away from the body. Results using a constant celerity are somewhat worse, on account of the reasons already given, and results using a rigid wall condition indicate significant reflected waves as expected. In frequency domain methods, the validity of a radiation condition has immediate effects on the accuracy of the solution and considerable research has been directed towards the development of a proper radiation condition at second order. In the present method, with the use of a large domain the solution is independent of the radiation condition for a sufficient duration of simulation and the adverse effects of the radiation boundaries on the solution can confidently be ruled out regardless whether or not the radiation condition is applied. 4.2.8 Irregularity of Second-Order Solution For the case of a surface-piercing body, there exists a so-called irregularity in the solution at second order. This is characterized by a discontinuity in the normal fluid velocity at the waterline contour at which the potential gradient calculated from the free surface boundary conditions does not agree with that specified by the body surface boundary condition. This does not necessarily occur only in the solution at second order, but is also present whenever the free surface boundary condition is inhomogeneous. Numerical tests have been carried out to investigate the behaviour of the potential at second order near the waterline contour. The results have demonstrated the discontinuity of potential gradient, although the potential itself is finite and continuous across the waterline contour. The discrepancy in potential gradient is an order smaller compared to the potential gradient away from the body. Since the potential at second order is finite and continuous across the waterline contour, the numerical solution is unique and non-singular. The effects of the irregularity on the free surface elevation and hydrodynamic pressure at second order which themselves do not depend on the potential gradient at second-order are expected to be small. 67 Now that the effects of the rank of the matrix equation, the initial conditions, the facet and time-step sizes, the radiation condition and the irregularity in the second-order solution have been considered, the method is ready to be applied to the three examples mentioned in § 4.1 and comparisons are made with previous theoretical and experimental results where appropriate. 4.3 Submerged Circular Cylinder in Two Dimensions 4.3.1 Free Surface Profiles Figure 12 shows the development with time of the free surface profiles to first and second order for ka = 0.5, h/a = 2, H /L = 0.1 and deep water. The domain is four wavelengths long and the cylinder axis is located at x/L = 2. At t = 0, the free surface profile corresponds to an undisturbed Stokes second-order wave train in the domain. With the imposition of the body surface boundary condition over the first cycle, scattered waves at first and second order are gradually generated and propagate away from the cylinder. A steady state solution is developed near the cylinder after the first cycle. In the figure, the most striking difference between the linear and nonlinear solutions is the formation of second-order free waves at twice the incident wave frequency in addition to the second-order forced waves associated with the first-order wave profile. These second-order free waves satisfy the linear dispersion relation (2.30) and propagate away from the cylinder at the corresponding celerity. This phenomenon has also been observed in laboratory experiments (e.g. Chaplin, 1984) and constitutes one of the most important nonlinear effects. In the vicinity of the body, these second-order free waves interact with the second-order forced waves and both contribute to the second-order wave force. The diffracted free surface profiles to first and second order for the same wave conditions as in figure 12 are compared with the corresponding incident wave profiles in figure 13. For 68 the case of deep water waves considered here, there is no reflection from the cylinder on the basis of both the first- and second-order predictions. At first order, the transmitted waves have exactly the same amplitude and waveform as the incident waves but suffer a phase shift in passing the cylinder, as reported by previous authors (e.g. Ogilvie, 1963). For the solution to second order, the transmitted waves are disturbed to a greater extent due to the generation of second-order free waves in the downwave region while the incident wave field in the upwave region remains undisturbed. This observation confirms the second-order solutions by Vada (1987) and Mclver & Mclver (1990), and the full nonlinear solution by Stansby & Slaouti (1984) in which the same problem was treated. 4.3.2 Wave Forces Wave Forces in Deep Water Figure 14 shows the development with time of the wave force components to first and second order on a submerged circular cylinder for the same wave conditions as in figure 12. At t = 0, the wave forces correspond to those obtained by the Froude-Krylov assumption in which the body is assumed to have no effects on the incident flow and a steady state solution is developed gradually in one wave cycle. Despite the generation of strong second-order free waves as indicated in figures 12 and 13, the wave force at second order is shown to be relatively small compared to the amplitude of the first-order force component. The results also indicate the presence of a steady vertical drift force upwards. However, no horizontal drift force is observed, which is in accordance with the known result that a submerged circular cylinder in deep water does not reflect first-order incident waves. Comparisons of the force components at first and second order have been made with a previous theoretical study (Vada, 1987) and two experimental studies (Chaplin, 1984; and Inoue & Kyozuka, 1988) and are presented in figures 15 to 17. For the deep water case considered here, the contribution from the velocity-squared term towards the second-order 69 oscillatory force is zero and F2 is evaluated from fa only. Each of F i and F2 refers both to the x and z components, whose amplitudes are equal for deep water, while F 0 refers to the component in the positive z direction only (the x component being zero). The computed first-order oscillatory force and second-order steady drift force are compared with Vada's (1987) theoretical results and Chaplin's (1984) experimental results in figures 15a and 15b respectively. Identical agreement is indicated and is expected. Figure 16 presents a comparison between the second-order oscillatory forces obtained by the two second-order methods discussed respectively in § 2.2 and § 2.3 with those of Vada (1987) and Chaplin (1984). The comparison between the results obtained by the Stokes second-order method and Vada indicates excellent agreement for all frequencies which is on account of both methods solving the identical boundary value problem. At low frequencies, the alternative second-order method gives sightly lower predictions in relation to the Stokes second-order method and Vada's numerical model, while the experimental results diverge from the three numerical predictions, as might be expected on account of the increasing importance of viscous effects at small values of ka. At high frequencies, the alternative second-order method tends to give higher predictions in relation to the solutions obtained by the Stokes second-order method and Vada, but is found to give better agreement with Chaplin's experimental results. In order to further ascertain the validity of the present second-order solutions, figure 17 presents a comparison of the results obtained by the two methods with the experimental results of Inoue & Kyozuka (1988) in which similar relations between the two theoretical solutions and the experimental results are observed. In general, the effects of the expansion terms in (2.50) are to give slighdy lower predictions of the second-order oscillatory force at the low frequency region and to give higher predictions at the high frequency range. However, it should also be noted that the second-order oscillatory force decreases much more rapidly in comparison to the first-order force at high frequencies, and the discrepancies between the two second-order solutions at this frequency range have insignificant effects on the total force. 70 Effects of finite water depth and depth of submergence For the case of a submerged circular cylinder in finite water depth, the amplitudes of the x and z components of the first- and second-order oscillatory forces are no longer equal. The contribution from the velocity-squared term towards the second-order oscillatory force is small but non-zero. The x and z components of the second-order oscillatory forces (denoted by F2X and F2Z) obtained by the two methods are shown in figures 18a and 18b respectively. The amplitudes of the x and z force components are almost identical for ka less than about 0.2 and the difference in amplitude between the two components increases with increasing values of ka. The discrepancies between the solutions obtained by the two methods are small and their relations generally follow the same trends as in the deep water case. In general, the solution obtained by the Stokes second-order method is relatively stable for the entire frequency range of interest, while the stability of the solution obtained by the alternative second-order method is found to deteriorate with decreasing ka. Figures 19 and 20 show respectively in deep water and in finite water depth the second-order oscillatory forces obtained by the two methods as functions of the depth of submergence. Both solutions are found to decrease rapidly for small values h/a and the discrepancies between the two solutions are small and further diminish with the depth of submergence. The solutions obtained by both methods are found to be stable for the range of h/a considered. For larger values of h/a, the second-order oscillatory force is small and numerical errors in the form of unsteadiness in the solution are more noticeable in the alternative second-order method since more terms are involved in the integral equation. Although in principle the expansion terms in (2.50) wil l give bad approximations when the cylinder is close to the still water surface, the relation between the two solutions is maintained even for small values of h/a. For h/a —» 1, the clearance between the cylinder surface and the still water surface tends to zero and the perturbation expansion approaches one of its limitations, so that a second-order solution is no longer sufficient to describe the nonlinearities involved. 71 4.3.3 Applicability of Two Methods In the alternative second-order method, the full nonlinear free surface boundary conditions and the integral equation defined on the instantaneous free surface are both expanded about the still water level by a Taylor series expansion whereas in the Stokes second-order method only the free surface boundary conditions are treated in this way. Through the use of the perturbation expansions, the second-order integral equation in the alternative second-order method is shown to contain additional second-order terms in comparison to that of the Stokes second-order method. These are associated with the application of the Taylor series expansion to the integral equation itself. However, the application of the alternative second-order method in more general conditions is problematic. Apart from the ill-defined radiation condition at second order, the expansion terms in (2.50) which require a considerable amount of computational effort and involve subtractions of large quantities constitute a major source of numerical error. Furthermore, a major shortcoming of the alternative second-order method which has not been indicated here is its application to surface-piercing bodies. Similar to the problem with the radiation condition, the expansion in (2.45) which provides the basis for the development of (2.50) gives bad approximations near the intersections between the still water surface and the body surface in the case of a surface-piercing body. The potential at second order and its derivatives near the intersections are found to depend on the nearby facet layout and the computed second-order wave forces and free surface elevation are not reliable. For the flow parameters considered, the difference between the results obtained by the two methods is generally small and has little effects on the total wave force. With its brevity in algebra and its comprehensiveness in treating the free surface boundary conditions and radiation condition to second order, the Stokes second-order method is recommended for general applications. 72 4.4 Semi-Submerged Circular Cylinder in Two Dimensions 4.4.1 Free Surface Profiles and Runup Figure 21 shows the development with time of the free surface profiles to first and second order for ka = 0.6, H / L = 0.1 and deep water. The domain is four wavelengths long and the cylinder axis is located at x/L = 2. The flow immediately adjacent to the body appears to repeat itself after the first cycle, while the flow further from the body takes somewhat longer to reach a steady state. The scattered waves at first and second order propagate steadily away from the body at their corresponding celerities. In the upwave region, the first-order scattered and incident waves propagate in the opposite direction forming a system of partial standing waves which gives rise to a strong second-order forcing. In the downwave region, the first-order scattered and incident waves are out of phase and propagate in the same direction. The resultant wave amplitude is small and the second-order forcing is weak. Figure 22 shows the corresponding variation with time of the free surface elevations to first and second order on the two sides of the body. Despite a short duration of transient effects associated with the initial conditions, the free surface elevation approaches its steady state after only one wave cycle. It is observed that second-order effects contribute significandy to the total free surface elevation on the upwave body surface and that a steady component of setup is also observed on the two sides of the body. The runup to second order and the amplitudes of the free surface elevations at first and second order on the upwave and downwave body surfaces are plotted in dimensionless form against the relative body size ka in figures 23a and 23b respectively. The steady setup r j 0 associated with the increase in mean water level near the body is due to the second-order interactions between the first-order incident and scattered waves. As ka -> 0, the runup to second order and the amplitudes of the free surface elevations at first and second order are expected to approach the values corresponding to regular incident waves, as if the body is not present. Thus, Irjjl/A should approach unity; T\Q/A should approach zero; IT^I/A should 73 approach the second-order wave amplitude; and R/A should approach the maximum free surface elevation of the second-order incident waves. The second-order interactions between the body and the surrounding wave field are observed to be most significant for small values of ka and ITI2' is relatively large, having a peak near ka = 0.1 on both the upwave and downwave body surfaces. The runup to second order in this frequency range is dominated by lrj 2l. For large values of ka, each component of the free surface elevations including the runup to second order approaches a constant value corresponding to the total reflection of incident waves. For the cases considered here, second order effects are important and account for about one-fifth to one-half of the runup to second order on the upwave body surface, depending on the frequency range of interest 4.4.2 Wave Forces Figure 24 shows the development with time of the wave force components to first and second order for the same wave conditions as in figures 21 and 22. At t = 0, the wave forces correspond to those obtained by the Froude-Krylov assumption and a steady state solution is developed gradually in one wave cycle. The results indicate the presence of strong second-order effects which include steady drift forces in the positive x direction and in the negative z direction. In figure 25, comparisons of the first-order oscillatory force components are made with Kyozuka's (1980) theoretical and experimental results and indicate excellent agreement as expected. In figure 26, comparisons of the second-order oscillatory force components are made with Kyozuka's (1980) theoretical and experimental results and also with Miao & Liu's (1989) theoretical results. The present results give good agreement with Kyozuka's theoretical results, whereas the oscillatory variation with ka predicted by Miao & L i u and indicated in figure 26a is not reproduced here. Kyozuka's experimental results exhibit some scatter but have the same trend as other theoretical results. 74 As indicated in (2.66), the second-order steady drift force is made up of two components, one due to the free surface fluctuation on the body and the other due to the velocity-squared term in the Bernoulli equation, both of which can be computed on the basis of the first-order potential. Figure 27 shows the steady drift forces and their component terms in the x and z directions. Since higher surface elevations occur on the upwave side of the body, the horizontal force due to free surface fluctuation is higher on the upwave body surface than on the downwave body surface and therefore the resultant drift force component due to the free surface fluctuation always acts in the direction of wave propagation. On the other hand, the higher surface elevations in the upwave region are associated with higher fluid velocities and hence lower pressure there, and consequendy the horizontal drift force component due to the velocity-squared term always acts opposite to the direction of wave propagation. The two components of the horizontal drift force are counteracting each other. Due to the balance of fluid momentum, the resultant drift force in the x direction always acts in the direction of wave propagation, and for large values of ka approaches a constant value corresponding to a total reflection of incident waves by the body. On the other hand, the free surface fluctuation has no effect on the steady drift force in the z direction, which itself always acts in the negative z direction due to the negative pressure induced by the velocity-squared term in the Bernoulli equation. The second-order oscillatory force consists of a contribution from the potential at second order as well as a contribution from the potential at first order. Figures 28a and 28b show the amplitudes of the second-order oscillatory forces and their components in the x and z directions respectively. In figure 28a, the two components due to the potential at first order in the x direction are of the same phase, while the component due to the potential at second order has the opposite phase over most of the frequency range, and consequendy gives rise to a smaller resultant force. Similar observations are also made in figure 28b for the second-order oscillatory force in the z direction, but the force component due to the potential at second order is found to dominate. It is interesting to note that for large values of ka both the x and z force 75 components due to the potential at second order increase steadily with ka, while the components due to the potential at first order remain fairly constant and do not depend on the size of the body. In terms of the contributions toward the total wave force, the potential at second order plays an important role especially for large values of ka and therefore should be accounted for accurately. 4.5 Surface-Piercing Circular Cylinder in Three Dimensions 4.5.1 Free Surface Profiles Flow Visualization For the numerical model shown in figure 3, approximately 100 kilobytes of numerical results of the free surface elevation are generated at each time step. The most effective way to present such a large volume of data is through a computer animation of the free surface movement. Not only does this provide a convenient visualization and analysis of the predicted wave field, but it also provides a visual assessment of the correctness of the numerical solution. For particular conditions corresponding to ka = 1, kd = 1 and H / L = 0.08, the variation with time of the free surface elevations to first and second order has been transformed to a colour film as an aid to viewing the development of the flow. Although numerical solutions were obtained on one side of the x axis, the free surface elevations were plotted over the entire domain to assist visualization of the flow development Numerical results of the free surface elevation were obtained using the I B M 3090/150S computer at the University of British Columbia and were transferred to a Macintosh IIx computer for post-processing. A spreadsheet program WINGZ Version 1.1 and an image analysis program I M A G E Version 1.33g were used to generate colour-contoured three-dimensional views of the free surface at 20 equal time-steps in sequence spanning a single wave period and shown on a 60 x 60 grid. These images were stored in disk and cycled 76 through a loop in order to generate a moving picture extending over many wave periods. The computer animation was then transferred to videotape using the RasterOps ColorBoard 364 and the RasterOps Video Expander that enable the conversion of a Red-Green-Blue video signal to a National Television Standards Committee (NTSC) interlaced colour signal without significant loss of resolution and colour detail. Similarly, plan views of the flow development can also be generated on the same basis. Free Surface to First Order By studying the first-order solution, some insight may be gained into the nature of the quadratic forcing terms in (2.28) and (2.29), which in turn dictate the second-order wave field. Two of the images corresponding to the free surface at first order are shown in figure 29. The vertical scale of the free surface elevation is exaggerated to make the interactions more pronounced. In the figure, the scattered waves propagate concentrically from the cylinder and their amplitudes decay with distance away from the cylinder. In the upwave region, the scattered waves and the incident waves propagate in the opposite direction forming a system of partial standing waves in front of the cylinder and giving rise to a strong second-order forcing. Along the sides, the first-order scattered waves propagate perpendicularly to the incident waves producing a fairly complicated forcing term near the cylinder. In the downwave region, the scattered waves with decaying amplitude propagate in approximately the same direction as the incident waves and the second-order forcing term is small. Free Surface to Second Order Two of the images corresponding to the free surface to second order are shown in figure 30. Figure 30(a) corresponds to the instant when the surface elevation at the upwave side of the cylinder, 8 = JC, is a maximum, and figure 30(b) corresponds to the instant half a cycle later. For the example shown, the maximum runup to second order is found to be about 40% higher than the first-order solution, and the free surface profiles clearly indicate the importance of second-order effects. 77 The second-order interactions in the downwave region are mainly free waves satisfying the linear dispersion relation (2.30) and propagating away from the cylinder at the corresponding celerity independent of the first-order wave system. From figure 30 and the associated video film, these second-order free waves can easily be identified to originate at the lee and near the lee quarter of the cylinder in which three-dimensional wave pulses are emitted alternatingly. Second-order free waves of smaller amplitudes are also generated along the upwave surface of the cylinder and propagate away from the cylinder concentrically, as may be discerned from the left-hand side of figure 30(b). On the other hand, the second-order forced waves are phased-locked with the first-order wave system to produce an overall wave profile with steeper crests and flatter troughs, which may be clearly seen. Consequendy, the runup to second order at 8 = 7t is effectively increased by the presence of the second-order forced waves. 4.5.2 Wave Runup The computed wave runup and rundown envelopes to second order are compared with the theoretical and experimental results of Kriebel (1987) for three different cases as shown in figures 31a, 31b and 31c respectively. The amplitudes of the free surface elevations at first and second order on the cylinder are also shown in the figures to indicate the wave activity along the cylinder circumference. The comparison between the two theoretical results indicates good agreement, while the comparison with the experimental result is generally good except for the rundown envelope near 8/7C = 0.55 - 0.15. The rundown in this region predicted by the two theoretical methods is significantly lower than the experimental results, which is attributed to the negative mean water level and the large amplitude of the free surface fluctuation at second order as indicated in the figures. As the difference between the theoretical and experimental results appears to decrease with increasing values of ka, the discrepancy is expected to be due to the local flow separation that occurs along the side of the cylinder for small values of ka at which viscous effects are more important 78 As indicated in figures 31a to 31c, the runup to second order and the amplitudes of the free surface elevations at first and second order are functions of the azimuthal angle 6 as well as the incident wave frequency. Due to the redirection of incident wave momentum, the maxima of IT^ and r j 0 always occur at the upwave side of the cylinder, 8 = TC. Since the second-order forced waves are phased-locked with the first-order waves, together with the second-order steady setup, the maximum of the runup to second order usually occurs at 9 = TC. On the other hand, due to the large tangential fluid velocity at first order that occurs along the side of the cylinder, r j 0 is always negative and has a minimum in this region. The amplitudes of r\2 a t 6 = TC and 6 = 0 are dominated respectively by forced and free waves, while along the side of the cylinder both forced and free waves are dominant and usually more than one local maximum might be found for large values of ka as indicated in figure 3 lc . In figure 32, the runup to second order and the amplitudes of the free surface elevations at first and second order at 6 = TC are plotted in dimensionless form as R/A against the incident wave frequency ka for a/d =1 and H/L = 0.1. It should be noted that for a fixed value of a/d, kd decreases as ka decreases, so that as ka —> 0 the perturbation expansion approaches its limit and the solution at second order becomes invalid. For small values of ka, both R and lT]2l are relatively large, on account of the importance of second-order effects in shallow water. For large values of ka, the runup to second order and the amplitudes of the free surface elevations at first and second order each approaches a constant value corresponding to a local reflection of incident waves at 8 = TC. For the two-dimensional case of a semi-submerged circular cylinder with axis at the still water level, the runup to second order and the amplitudes of the free surface elevations on the upwave body surface appear to approach approximately the same values for large cylinders (figure 23a). For the entire frequency range considered here, r\2 is in phase with Tj l 5 so that these taken together with T| 0 give rise to a higher runup to second order. Second-order effects are relatively important and account for at least one-fifth of the runup to second order. 79 As indicated in figures 31a to 31c, the maxima of \r\2\ occur at various locations along the side of the cylinder, depending on the values of ka and kd. In order to illustrate the second-order wave behaviour along the circumference of the cylinder, figure 33 shows the maximum amplitude of r j 2 over 9, in dimensionless form lT | 2 la /A 2 , together with the corresponding location 9 m a x of the maximum, as functions of ka for a/d = 1. The corresponding results of K i m & Yue (1989) are also shown in the figure. Since more than one local maximum might occur along the side of the cylinder, the curve joining the points for 0 m a x is not expected to be continuous, while the curve joining the points for ( lT| 2 la /A 2 ) m a x is not expected to be smooth. The comparison between the two solutions indicates good agreement for the locations of the maxima, while a small discrepancy is observed for the maxima at higher values of ka. With the exception of ka = 0.75, the amplitude of r j 2 for various values of ka is found to be a maximum near 9/7C = 0.25 - 0.5, which is associated with the generation of the three-dimensional wave pulses in the downwave region. 4.5.3 Wave Forces and Moments For the present case of a bottom-mounted vertical circular cylinder, onlythe force components in the x direction and the moment components about the y axis are non-zero and are considered in this section. Figure 34 shows the development with time of the wave forces and moments to first and second order for ka = 1, kd = 1 and H/L = 0.08 in which the steady state solution is developed gradually in one wave cycle and the results indicate the presence of strong second-order effects which include a steady drift force and moment. Figure 35 shows the wave force and moment components at first and second order as functions of ka for a/d = 1 and H / L = 0.1. Also shown in the figure are the force maxima to second order in the positive and negative x directions, denoted by F x and F . x respectively, and the moment maxima in the clockwise and counter-clockwise directions about the y axis, denoted by M y and M . y respectively. To illustrate the relative importance of each component, 80 all force and moment components are made dimensionless with respect to pga 2 A and pga 3A respectively. For small values of ka, the second-order oscillatory force F2 is in phase with the first-order force F i , and thus together with the second-order steady force Fo these give rise to a relatively large resultant force to second order. On the other hand, for larger values of ka the first-order force F i and the second-order oscillatory wave force F2 are out of phase. Despite a second-order steady drift force in the positive x direction, the resultant force F x is found to be smaller than the first-order solution. In order to illustrate the structure of the solution at second order, and to provide a comparison with previous results, the second-order oscillatory force and moment components are plotted in figures 36(a) and 36(b) respectively, as functions of ka for a/d = 1. Eatock Taylor & Hung's (1987) analytic results based on an asymptotic form of the potential at second order which have been confirmed by Abul-Azm & Williams (1988) and K i m & Yue (1989) based on similar approaches are also shown in the figures. The comparison indicates good agreement except for large values of ka at which the two solutions diverge slightly from each (2) (2) other. For small values of ka, the components F 2 and M 2 due to the potential at second order are relatively large and in phase with the components F ^ and due to the potential at first order, so that the resultant second-order oscillatory forces F2 and M2 are even larger. For (2) (2) larger values of ka, the components F 2 and M 2 become out of phase with the components F ^ and M ^ 2 \ such that the resultant second-order oscillatory forces are smaller than the components due to the potential at second order. The force components due to the potential at second order constitute a major factor in determining both the amplitudes and phases of the resultant second-order oscillatory forces, and their amplitudes relative to other second-order quantities increase with increasing frequency. This invalidates many so-called engineering estimates of second-order wave effects (e.g. Herfjord & Nielsen, 1986) in which the potential at second order is ignored in favour of quadratic contributions of the first-order potential. 81 CHAPTER 5 CONCLUSIONS 5.1 Theoretical and Numerical Procedures A time-domain second-order method has been developed to treat the nonlinear diffraction of ocean waves around large fixed bodies of arbitrary shape in two and three dimensions. By solving the second-order boundary value problem in the time domain, the present method has less complicated algebra compared to most conventional second-order methods and is more versatile. It allows the study of deterministic irregular wave diffraction without solving the system of numerical equations every time-step. For the case of regular wave diffraction to second order, the proposed solution is validated by convergence tests of the numerical results and by comparisons with a number of previous analytic, numerical and experimental results. The major difficulties associated with the full nonlinear problem arise primarily from the two nonlinear free surface boundary conditions applied on the instantaneous free surface, which itself is unknown a priori, and from the poorly defined radiation condition in the time domain. Based on the Stokes second-order expansion procedure, the nonlinear free surface boundary conditions defined on the instantaneous free surface are first expanded about the still water level by a Taylor series expansion and terms up to second order are retained through a perturbation procedure. The field solutions at first and second order defined on a time-independent boundary which includes the still water surface are then obtained by the application of a boundary integral equation. Since the boundary is invariant in time, the solution to the system of simultaneous equations obtained through the discretization of the integral equation is required only once rather than at each time step. The flow potential is separated into a known incident potential and a scattered potential which is to be determined. The Sommerfeld radiation condition applied to the scattered potential is modified to include a 82 time-dependent celerity in order to account for the motions of forced and free waves at second order. The initial conditions correspond to a fully developed second-order wave field in the domain, and the scattered potential is developed in time and space through the gradual imposition of the body surface boundary condition over one wave period. The free surface boundary conditions and the radiation condition to second order are then satisfied by a numerical integration in time. Generally, a steady state solution is obtained after the first cycle. Although the computer program for the three-dimensional problem which generates and solves the matrix equation is computationally intensive, the matrix equation involved in the formulation is a function of geometry and discretization only. Results generated by the program can be stored in the disk memory and used subsequendy in the time-stepping program for different incident wave conditions. The application of the present method with a large number of facets is economically feasible. Since the rank of the matrix equation for the three-dimensional problem is as large as 2200, the accuracy of the matrix solution is first assessed in terms of the condition number. It is found out that the matrix equation is very well-conditioned and that the loss in accuracy in solving the system of equations is negligible. The limit on the rank of the matrix equation is basically imposed by the size of the computer memory which is much higher than that required in the present application. The time-stepping procedure is assessed in terms of the maximum time-step sizes allowed for a given facet size in order to maintain stable solutions at first and second order. For time-step sizes above this limit the solutions exhibit an oscillatory instability while for smaller time-step sizes the solutions remain stable. The numerical accuracy of the wave forces and runup is found to be relatively independent of the choice of time-step size, but depends more significandy on the number of facets used to represent the body. The convergence of second-order quantities with increasing number of facets is found to be slow compared to their first-order counterparts. The effectiveness of the proposed radiation condition is tested and verified. At first order, the use of a constant celerity satisfying the linear dispersion relation is found to be sufficient. At second order, the use a constant celerity corresponding to the second-order 83 free waves is shown to be insufficient and a time-dependent celerity is found to give the most satisfactory result A n alternative second-order method based on a different expansion procedure has also been developed in which the nonlinear free surface boundary conditions and the integral equation defined on the instantaneous free surface are both expanded about the still water level by a Taylor series expansion and terms up to second order are retained through a perturbation procedure. The two approaches give rise to identical first-order problems, but give rise to second-order problems which are apparently different. The major discrepancy arises from the second-order integral equation in which more second-order terms are retained through the Taylor series expansion. However, its application in more general conditions is problematic and is restricted for submerged bodies only. Through the various comparisons, it has been shown that the differences between the two second-order solutions are generally small and have insignificant effects on the total wave forces. The second-order interaction is apparently dominated by the nonlinear terms in the free surface boundary conditions rather than the expansion terms in the second-order integral equation. With its brevity in algebra and its comprehensiveness in treating the free surface boundary conditions and radiation condition to second order, the Stokes second-order method is recommended for general applications. 5.2 Numerical Studies The Stokes second-order method is applied to the study of regular wave diffraction to second order around a fully submerged and a semi-submerged circular cylinder in two dimensions, and around a bottom-mounted surface-piercing circular cylinder in three dimensions. Comparisons of the computed wave forces and runup are made with previous theoretical and experimental results where appropriate. For a fully submerged circular cylinder in deep water, there is no reflection of incident waves from the cylinder on the basis of both the first- and second-order predictions. At first 84 order, the transmitted waves have exactly the same amplitude and waveform as the incident waves but suffer a phase shift in passing the cylinder. For the solution to second order, the transmitted waves are disturbed to a greater extent due to the generation of second-order free waves in the downwave region while the incident wave field in the upwave region remains undisturbed. Since the cylinder does not reflect incident waves, the second-order steady drift force in the x direction is zero and a resultant drift force acts in the positive z direction. The amplitudes of the x and z force components are equal at first and second order, and the contribution from the velocity-squared term towards the second-order oscillatory force is zero. Comparisons of the computed wave force components at first and second order are made with a previous theoretical study in which identical agreement is indicated, while comparisons with previous experimental results further reinforce the validity of the present method. The runup to second order and the amplitudes of the free surface elevations at first and second order on a semi-submerged circular cylinder in deep water are presented. The second-order interactions between the body and the surrounding wave field are observed to be most significant for low frequencies and the total runup in this frequency range is dominated by second-order effects. For higher frequencies, each component of the free surface elevations as well as the total runup, approaches a constant value corresponding to the total reflection of incident waves. For the cases considered here, the free surface elevation at second order is important and second-order effects account for as much as 50% of the total runup. The wave force components at first and second order are compared with previous theoretical and experimental results and indicate favourable agreement. The second-order oscillatory force components due to the potentials at first and second order are generally out of phase which gives rise to a smaller resultant force. However, at high frequencies the force component due to the potential at second order is important and increases steadily with incident wave frequency while the component due to the potential at first order remains fairly constant. For the case of a surface-piercing circular cylinder extending from the seabed in three dimensions, the variation with time of the free surface elevations to first and second order is 85 transformed to a colour film as an aid to viewing the development of the flow. Various features of the free surface to second order are identified from the animation. The runup to second order and the amplitudes of the free surface elevations at first and second order are also presented in which second-order effects are shown to be important in evaluating the total runup. At high frequencies a local reflection of incident waves at the upwave side of the cylinder is observed and the results are found to be consistent with those obtained in the two-dimensional problem. Despite a second-order steady drift force in the positive x direction, in the high frequency range the force maxima to second order in that direction are shown to be smaller than the first-order force amplitude. The free surface elevation and wave force components at second order are compared with previous theoretical results and good agreement is obtained except for high frequencies at which the two solutions diverge slightly from each other. 5.3 Recommendations for Further Study In this thesis, a time-domain method is developed to solve the second-order wave diffraction problem. To illustrate and examine the present method, and to compare with known theoretical and experimental results, only the special case of regular wave diffraction around fixed structures of simple geometry is considered. With the present method validated, applications can readily be made to a variety of practical and academic problems in coastal and ocean engineering. Without any modification of the present method, design curves for nonlinear wave forces and runup for a variety of large coastal and offshore structures such as sea-walls, berths, ship sections, gravity platforms and artificial islands can be generated. Nonlinear diffraction of ocean waves around more complicated offshore structures such as semi-submersibles or tension-leg platforms can also be studied. The effects of seabed topography and interferences due to adjacent structures on the second-order wave forces and runup can be investigated. 86 Since the present method solves the second-order diffraction problem in the time domain, the effects of deterministic irregular wave profiles on the wave forces and runup can be examined direcdy without resorting to a stochastic analysis. The diffraction of a random wave train to second order and the corresponding statistical results are also of particular interest. The time series of a second-order incident potential for a specified wave spectrum can be generated by a bi-spectral analysis (e.g Hasselmann, 1962; and Longuet-Higgins, 1963). The results may then be applied to the diffraction model to generate the time series of the wave forces and runup at second order which in turn can be used to evaluate the corresponding spectra and statistical data. On the other hand, the diffraction model can also be used to calculate the wave forces and runup at second order for incident wave trains composed of two fundamental frequencies, and these results can subsequently be used to construct the corresponding quadratic transfer functions. The spectra of the wave forces and runup at second order can then be evaluated by a numerical integration involving the quadratic transfer functions and the incident wave spectrum. The same approaches can be applied to the second-order diffraction of regular or irregular short-crested waves. It should be straightforward to extend the present method to include the effects of small amplitude body motions in waves, in which only the body surface boundary condition needs modifications and all the second-order effects would be accounted for by the numerical integration of the free surface boundary conditions. The effects of arbitrary body motions or manoeuvres in waves which are equivalent to wave-current interactions with fixed bodies can also be incorporated in the method through the use of a coordinate system attached to the mean motion of the body and the free surface boundary conditions should then be modified accordingly. Eventually, it is expected that the present method would be extended and applied to the time-domain simulation of the nonlinear responses of floating bodies moving freely or moored in regular or irregular waves. 87 References Abul-Azm, A . G . & Williams, A . N . 1988. Second-order diffraction loads on truncated cylinders. Journal of Waterway, Port, Coastal & Ocean Engineering, A S C E , 114(4), 436-454. Angell, T. S., Hsiao, G. C. & Kleinman, R. E. 1986. A n integral equation for the floating body problem. Journal of Fluid Mechanics, 166, 161-171. Brorsen, M . & Larsen, J. 1987. Source generation of nonlinear gravity waves with the boundary integral equation method. Coastal Engineering, 11, 93-113. Chakrabarti, S.K. 1987. Hydrodynamics of Offshore Structures. Springer-Verlag, Heidelberg, 440 pp. Chan, R . K . - C . 1977. Finite difference simulation of the planar motion of a ship. Proceedings of the Second International Conference on Numerical Ship Hydrodynamics, Berkeley, C A , 39-52. Chaplin, J.R. 1984. Nonlinear forces on a horizontal cylinder beneath waves. Journal of Fluid Mechanics, 147, 449-464. Chen, M . C . & Hudspeth, R. T. 1982. Nonlinear diffraction by eigenfunction expansions. Journal of Waterway, Port, Coastal & Ocean Division, ASCE, 108(WW3), 306-325. Demirbilek, Z. & Gaston, J.D. 1985. Nonlinear wave loads on a vertical cylinder. Ocean Engineering, 12(5) 375-385. Dommermuth, D.G. & Yue, D.K.P. 1987. Numerical simulations of nonlinear axisymmetric flows with a free surface. Journal of Fluid Mechanics, 178,195-219. Dubrulle, A . A . , Scarborough, R.G. & Kolsky, H.G. 1985. How to write good vectorizable F O R T R A N . Palo Alto Scientific Center Report No. G320-3478, I B M Corporation, 77 pp. Eatock Taylor, R. & Hung, S.M. 1987. Second order diffraction forces on a vertical cylinder in regular waves. Applied Ocean Research, 9(1), 19-30. Eatock Taylor, R., Hung, S .M. & Chau, F.P. 1989. On the distribution of second order pressure on a vertical circular cylinder. Applied Ocean Research, 11(4), 183-193. Garrison, C.J. 1984. Nonlinear wave loads on large structures. Proceedings of the Third International Symposium on Offshore Mechanics & Arctic Engineering, New Orleans, L A , Vol . 1, 128-135. Ghalayini, S.A. & Williams, A . N . 1989. Nonlinear wave forces on vertical cylinders of arbitrary cross section. Journal of Waterway, Port, Coastal and Ocean Engineering, ASCE, 115(6), 809-930. Grosenbaugh, M . A . & Yeung, R.W. 1989. Nonlinear free-surface flow at a two-dimensional bow. Journal of Fluid Mechanics, 209, 57-75. 88 Han, T .Y . , Meng, J.C.S. & Innis, G.E. 1983. A n open boundary condition for incompressible stratified flows. Journal of Computational Physics, 49,276-297. Hasselmann, K . 1962. On the non-linear energy transfer in a gravity wave spectrum. Journal of Fluid Mechanics, 12, 481-500. Herfjord, K . and Nielsen, F .G. 1986. Non-linear wave forces on a fixed vertical cylinder due to the sum frequency of waves in irregular seas. Applied Ocean Research, 8(1), 8-21. Hogben, N . & Standing, R.G. 1974. Wave loads on large bodies. Proceedings of the International Symposium on the Dynamics of Marine Vehicles and Structures in Waves, London, 258-277. Huang, H . , L i , J. & Wang, X . 1989. Asymptotic solutions and radiation conditions for second-order diffraction waves. Journal of Offshore Mechanics and Arctic Engineering, A S M E , 111, 203-207. Hunt, J.N. & Baddour, R.E. 1981. The diffraction of nonlinear progressive waves by a vertical cylinder. Quarterly Journal of Applied Mathematics, 34(1), 69-89. Hunt, J .N. & Williams, A . N . 1982. Nonlinear diffraction of Stokes water waves by a circular cylinder for arbitrary uniform depth. Journal de Micanique Thiorique et Appliqude, 1(3), 429-449. Inoue, R. & Kyozuka, Y . 1988. Nonlinear wave forces acting on submerged horizontal cylinders. Journal of Offshore Mechanics and Arctic Engineering, A S M E , 110, 62-70. Isaacson, M . 1982. Nonlinear wave effects on fixed and floating bodies. Journal of Fluid Mechanics, 120, 267-281. Isaacson, M . & Zuo, Q.H. 1989. Nonlinear wave forces on a circular cylinder. Canadian Journal of Civil Engineering, 16(2), 182-187. Jagannathan, S. 1988. Nonlinear Free Surface Flows and the Application of the Orlanski Boundary Condition. International Journal of Numerical Methods in Fluids, 8,1051-1070. Jamieson, W.W., Mansard, E.P.D. & Mogridge, G.R. 1985. Wave loading on a conical gravity platform. Proceedings of the Forth International Conference of the Behavior of Offshore Structures, BOSS'85, Delft, 673-684. K i m , M . H . & Yue, D.K.P . 1988. The nonlinear sum-frequency wave excitation and response of a tension-leg platform. Proceedings of the Fifth International Conference on the Behavior of Offshore Structures, BOSS'88, Trondheim, Norway, 687-704. Kim, M . H . & Yue, D.K.P . 1989. The complete second-order diffraction solution for an axisymmetric body. Part 1. Monochromatic incident waves. Journal of Fluid Mechanics, 200, 235-264. K i m , M . H . & Yue, D.K.P . 1990. The complete second-order diffraction solution for an axisymmetric body. Part 2. Bichromatic incident waves and body motions. Journal of Fluid Mechanics, 211, 557-593. 89 Korsmeyer, F.T., Lee, C.-H. , Newman, J.N. & Sclavounos, P.D. 1988. The analysis of wave effects on tension-leg platforms. Proceedings of the Seventh International Conference on Offshore Mechanics and Arctic Engineering, Houston, Texas, 2, 1-14. Kriebel, D.L. 1987. A second-order diffraction theory for wave runup and wave forces on a vertical circular cylinder. Ph.D. Thesis, University of Florida, Gainseville, Florida, 236 pp. Kriebel, D .L . 1990. Nonlinear wave interaction with a vertical circular cylinder. Part I: diffraction theory. Ocean Engineering, 17(4), 345-377. Kyozuka, Y . 1980. Non-linear hydrodynamic forces acting on two-dimensional bodies (First report, diffraction problem). Journal of the Society of Naval Architecture of Japan, 148, 45-53. Kyozuka, Y . 1981a. Non-linear hydrodynamic forces acting on two-dimensional bodies (Second report, radiation problem). Journal of the Society of Naval Architecture of Japan, 149, 47-53. Kyozuka, Y . 1981b. Non-linear hydrodynamic forces acting on two-dimensional bodies (Third report, effects of angles of the intersection of the body and the free-surface). Journal of the Society of Naval Architecture of Japan, 150, 166-174. Kyozuka, Y . 1982. Experimental study on second-order forces acting on cylindrical body in waves. Proceedings of the Fourteenth Symposium on Naval Hydrodynamics, Michigan, 73-136. Kyozuka, Y . 1983. Non-linear hydrodynamic forces acting on two-dimensional bodies (Fourth report, motions of a floating body in waves). Journal of the Society of Naval Architecture of Japan, 152, 148-158. Kyozuka, Y . 1988. Second-order wave forces acting on a horizontal circular cylinder in irregular waves. Nonlinear Water Waves, edited by K. Horikawa & H . Maruo. Springer-Verlage Berlin Heidelberg, 261-273. Liggett, J.A. & Liu , P. L . -F . 1983. The Boundary Integral Equation Method for Porous Media Flow. Allen & Unwin, Boston, 255pp. Lighthill, J. 1979. Waves and hydrodynamic loading. Proceedings of the Second International Conference of the Behavior of Offshore Structures, BOSS79, London, 1-40. Lee, J.F. & Leonard, J.W. 1987. A time-dependent radiation condition for transient wave-structure interactions. Ocean Engineering, 14 (6), 469-488. Lin , W.N. , Newman, N J . & Yue, D.K. 1984. Nonlinear forced motions of floating bodies. Proceedings of the Fifteenth Symposium on Naval Hydrodynamics, Hamburg, 33-49. Longuet-Higgins, M.S. & Cokelet, E.D. 1976. The deformation of steep surface waves on water. I. A numerical method of computation. Proceedings of the Royal Society London, Ser. A , 350, 1-25. Longuet-Higgins, M.S. 1963. The effects of non-linearities on statistical distributions in theory of sea waves. Journal of Fluid Mechanics, 17,459-480. 90 Masuda, K . & Nagai, T. 1991. Nonlinear wave forces on a pair of vertical cylinders. Journal of Offshore Mechanics and Arctic Engineering, A S M E , 113,1-8. Matsui, T. 1989. Computation of slowly varying second order hydrodynamic forces on floating structures in irregular waves. Journal of Offshore Mechanics and Arctic Engineering, A S M E , 111, 223-232. Mclver, M . & Mclver, P. 1990. Second-order wave diffraction by a submerged circular cylinder. Journal of Fluid Mechanics, 219, 519-529. Mei , C C . 1983. The Applied Dynamics of Ocean Surface Waves. John Wiley and Sons, New York, 740 pp. Miao, G.P. & Liu , Y . Z . 1989. A theoretical study on the second-order wave forces for two-dimensional bodies. Journal of Offshore Mechanics and Arctic Engineering, A S M E , 111, 37-42. Molin, B . 1979. Second-order diffraction loads upon three-dimensional bodies. Applied Ocean Research, 1(4), 197-202. Molin, B. & Marion, A . 1986. Second-order loads and motions for floating bodies in regular waves. Proceedings of the Fifth International Symposium on Offshore Mechanics & Arctic Engineering, Toyko, Vol . 1, 353-360. Newman, J.N. 1990. Second-harmonic wave diffraction at large depths. Journal of Fluid Mechanics, 213, 59-70. Ogilvie, F.T. 1963. First- and second-order forces on a cylinder submerged under a free surface. Journal of Fluid Mechanics, 16,451-472. Orlanski, I. 1976. A simple boundary condition for unbounded hyperbolic flows. Journal of Computational Physics, 21, 251-269. Papanikolaou, A . 1984. On calculations of nonlinear hydrodynamic effects in ship motion. Schiffstechnik, 31(3), 91-123. Pinkster, J.A. & Wichers, J.E.W. 1987. The statistical properties of low-frequency motions of nonlinearly moored tankers. Proceedings Offshore Technology Conference, Houston, Paper No. OTC 5457, 317-331. Rahman, M . 1984. Wave diffraction by large offshore structures: an exact second-order theory. Applied Ocean Research, 6(2), 90-100. Rahman, M . 1987. A design method of predicting second-order wave diffraction caused by large offshore structures. Ocean Engineering, 14(1), 1-18. Raman, H . & Venkatanarasaiah, P. 1976. Forces due to nonlinear waves on vertical cylinders. Journal of the Waterways, Harbors and Coastal Engineering Division, A S C E , 102(WW3), 301-316. Sarpkaya, T. & Isaacson, M . 1981. Mechanics of Wave Forces on Offshore Structures. Van Nostrand Reinhold, New York, 651 pp. 91 Shimada, K . 1987. Solutions to three-dimensional second-order diffraction problems by means of simple-source integral-equation method. Journal of the Society of Naval Architecture of Japan, 161, 139-145. Stansby, P.K. & Slaouti, A . 1984. On non-linear wave interaction with cylindrical bodies: a vortex sheet approach. Applied Ocean Research, 6(2), 108-115. Turing, A . M . 1948. Rounding off errors in matrix processes. Quarterly Journal of Mechanics and Applied Mathematics, 1,287-308. Vada, T. 1987. A numerical solution of the second-order wave-diffraction problem for a submerged cylinder of arbitrary shape. Journal of Fluid Mechanics, 174,23-37. Vinje, T. & Brevig, P. 1981. Numerical calculations of forces from breaking waves. Proceeding of the International Symposium on Hydrodynamics in Ocean Engineering, Trondheim, Norway, 547-566. Vinje, T., Maogan, X . & Brevig, P. 1982. A numerical approach to nonlinear ship motion. Proceedings of the Fourteenth Symposium on Naval Hydrodynamics, Washington, 245-278. Wang, J.L., Machemehl, J.L. & Ay, L . L . 1983. Design criteria for wave run-up on artificial islands. Proceedings Offshore Technology Conference, Houston, Paper No. OTC 4580, 9-16. Wu, G . X . & Eatock Taylor, R. 1990. The second order diffraction force on a horizontal cylinder in finite water depth. Applied Ocean Research, 12(3), 106- 111. Williams, A . N . , Abul-Azm, A . G . & Ghalayini, S.A. 1990. A comparison of complete and approximate solutions for second-order diffraction loads on arrays of vertical circular cylinders. Ocean Engineering, 15(5), 427-445. Yamaguchi, M . & Tsuchiya, Y . 1974. Nonlinear effect of waves on wave pressure and wave force on a large cylindrical pile. Proceedings of Civil Engineering Society in Japan, 229, 41-53. Yamashita, S. 1977. Calculations of the hydrodynamic forces acting upon thin cylinders oscillating vertically with large amplitude. Journal of the Society of Naval Architecture of Japan, 141, 67-76. Yen, S.M. & Hall , D.R. 1981. Implementation of open boundary conditions for nonlinear free-surface wave problems. Proceedings of the Third International Conference on Numerical Ship Hydrodynamics, Paris, 163-176. 92 L/AS No iteration, (T/At)m First order Second order 1 iteration, (J/At)^ First and second order 10 13 15 7 20 19 21 9 30 28 30 12 40 31 33 13 50 34 36 15 Table 1. Minimum T/At as a function of L/AS for stable solutions of regular wave diffraction at first and second order around a semi-submerged circular cylinder for ka = 1 and deep water. N b h/a = 1.75 h/a = 2.00 Fi/pgaA Fo/pgA 2 F2 /pgA 2 Fi/pgaA Fo/pgA 2 F V p g A 2 20 1.3328 0.3696 0.5465 1.1557 0.2556 0.2929 30 1.3329 0.3869 0.5642 1.1559 0.2673 0.2986 40 1.3330 0.3932 0.5686 1.1560 0.2715 0.2991 50 1.3330 0.3960 0.5682 1.1560 0.2735 0.2988 60 1.3330 0.3976 0.5671 1.1560 0.2747 0.2980 70 1.3330 0.3986 0.5669 1.1560 0.2753 0.2979 80 1.3330 0.3992 0.5668 1.1560 0.2760 0.2979 Table 2. Amplitudes of wave force components at first and second order on a submerged circular cylinder as functions of Nb for ka = 0.4 and deep water. 93 Upwave Surface Downwave Surface N b hlil/A TJoa/A2 ft\2\a/A2 Hill/A Tioa/A2 lri 2la/A 2 10 1.8366 0.4899 0.6058 0.4246 0.0336 0.2346 15 1.8419 0.4942 0.7859 0.4202 0.0365 0.2554 20 1.8474 0.5044 0.9295 0.4208 0.0392 0.2600 25 1.8494 0.5067 1.0024 0.4207 0.0400 0.2634 30 1.8498 0.5092 1.0352 0.4204 0.0412 0.2683 35 1.8520 0.5147 1.0976 0.4209 0.0418 0.2661 40 1.8512 0.5148 1.0890 0.4206 0.0424 0.2734 Table 3. Amplitudes of free surface elevations at first and second order on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water. 94 N b Fi x /pgaA Fox/pgA2 F(2x/pgA2 F S / p g A 2 F ^ p g A 2 10 1.1631 0.3716 1.3482 1.1040 0.3721 15 1.1764 0.3740 1.3574 1.1701 0.3918 20 1.1804 0.3814 1.3612 1.1897 0.4316 25 1.1826 0.3860 1.3606 1.2010 0.4500 30 1.1839 0.3889 1.3587 1.2046 0.4667 35 1.1844 0.3928 1.3593 1.2096 0.4783 40 1.1852 0.3934 1.3571 1.2060 0.4763 Table 4. Amplitudes of wave force components at first and second order in the x direction on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water. N b Fiz/pgaA Foz/pgA2 F^i / pgA 2 Fg / pgA 2 • F 2 z / p g A 2 10 1.0441 -0.3622 0.3433 1.1356 0.8022 15 1.0433 -0.3612 0.3423 1.1776 0.8444 20 1.0430 -0.3602 0.3411 1.2025 0.8689 25 1.0430 -0.3599 0.3406 1.2127 0.8790 30 1.0429 -0.3597 0.3404 1.2135 0.8797 35 1.0429 -0.3595 0.3401 1.2150 0.8810 40 1.0429 -0.3595 0.3401 1.2163 0.8820 Table 5. Amplitudes of wave force components at first and second order in the z direction on a semi-submerged circular cylinder as functions of Nb for ka = 0.6 and deep water. 95 N b T l o a / A 2 iTfcla/A* 0max/ft 0 = 71 0 = 71 0 = 7C 9 = 0max 0 = 0 36 1.6694 0.5220 2.0658 2.9653 1.8783 0.2807 64 1.6831 0.4998 . 2.4289 3.0087 1.9490 0.3109 100 1.6926 0.5162 2.5594 3.1604 1.9850 0.2947 Table 6. Amplitudes of free surface elevations at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1. N b Fj /pga^ Fo/pgaA2 F^VpgaA2 F^VpgaA2 F2/pgaA2 36 3.1864 0.9615 1.0383 1.5225 1.9970 64 3.2339 0.9813 1.0704 1.4592 2.0150 100 3.2509 0.9970 1.0800 1.4623 2.0249 Table 7. Amplitudes of wave force components at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1. N b Mi/pga3A Mo/pga 2A 2 M 2 Vpga 2 A 2 M ( 2 ) /pga 2 A 2 M2/pga2A: 36 1.6928 0.9265 1.1053 0.9732 1.8084 64 1.7283 0.9416 1.1135 0.9542 1.8473 100 1.7407 0.9569 1.1267 0.9674 1.8778 Table 8. Amplitudes of wave moment components at first and second order on a surface-piercing circular cylinder as functions of Nb for ka = 1 and kd = 1. 96 Figure 1. Definition sketch. (a) (b) Figure 2. Examples of surface discretization in two dimensions, (a) fully submerged circular cylinder with N b = 30, N = 30, N o = 76 and N = 136. (b) semi-submerged circular cylinder with N. = 15, N =30 ,N =64andN=109. b ' c ' o 97 Figure 3. Examples of surface discretization of a bottom-mounted surface-piercing circular cylinder in three dimensions, (a) N f a = 144, N c = 312, N o = 1368 and N = 1824. (b) N b = 220, N =312,N = 1402 and N = 1934. c ' o 98 q d I i i i ; i i i i 1 100 200 300 400 500 Rank of Matrix Equation Figure 4. CPU time as a function of the rank of the matrix equation for the two-dimensional problem. , GEOM2D; , WAVE2D (100 time steps). 1400 1600 1800 2000 2200 2400 Rank of Matrix Equation Figure 5. CPU time as a function of the rank of the matrix equation for the three-dimensional problem. , GEOM3D; , WAVE3D (100 time steps). 99 fi 0) _Q O d o E *~ c "a c o o o o I ' i I I I I I i 1 1 0 100 200 300 400 500 600 Rank of Matrix Equation Figure 6. Condition number as a function of the rank of the matrix equation for the two-dimensional problem. , submerged circular cylinder; , semi-submerged circular cylinder. CM I ' 1 p- 1 1 1 ' 1 ' 1 o CD <<- CO C .O TD C o C J co ro 500 1000 1500 2000 Rank of Matrix Equation 2500 Figure 7. Condition number as a function of the rank of the matrix equation for the three-dimensional problem. 1 0 0 o fO | I , I 0.0 1.0 2.0 3.0 4.0 o o • I I I I I L . I I 0.0 1.0 2.0 3.0 4.0 o ro j 1 1 1 r o -^ I I I . I I - J I 0.0 1.0 2.0 3.0 4.0 xA Figure 8. Free surface profiles to first and second order around a semi-submerged circular cylinder obtained by different modulation times T m for ka = 0.6, H/L = 0.1, t/T = 5 and deep water. , solution to first order; , solution to second order. 101 in • o • ! 1 -CM o in i I . J i o r- /I — o * o (a) n i • en 0.0 0.5 1.0 1.5 2 t/T2 O Y i i i i i i 1 0.0 0.5 1.0 1.5 2.0 t/T2 Figure 9. Celerity and flow potential of a bichromatic wave train at x = 0 as functions of time for a^aj = 3, c / c 2 = 0.5 and deep water, (a) celerity, (b) flow potential. -, theoretical solution; , numerical solution. 102 o in d q d J — r n- 1 ' T (a) 1 1 — 1 1 Figure 10. Celerity and flow potential of a Dichromatic wave train at x = 0 as functions of time for a/aj = 2, c/Cj = 0.8 and deep water, (a) celerity, (b) flow potential. , theoretical solution; , numerical solution. 103 Figure 11. Development with time of oscillatory force components at first and second order on a semi-submerged circular cylinder subject to various treatments at radiation boundaries for ka = 0.6 and deep water, (a) components due to tyv (b) components due to §r , hJL = 4; , Lp/L = 2 with time-dependent celerity; , LJL = 2 with constant celerity; , Lp/L = 2 without radiation condition. 104 : o o cs o CS cs I o cs" > ° cs I o cs >: cs q cs* >: cs I 1 1 p . 1 1 0.0 1.6 2.0 3.0 4. *N*A = 2 1 1 0.0 1.0 2.0 3.0 4. * \ t/T = 3 I 1 1 i i 1 1 -0.0 1.0 2.0 3.0 4. = \ t/T = 4 1 1 1 i 1 1 1 i -0.0 1.0 2.0 3.0 4. *v t/T = 5 1 1 1 i 1 1 1 / \ > ^ 1 -0.0 1.0 2.0 3.0 4. x / L Figure 12. Development with time of free surface profiles to first and second order around a submerged circular cylinder for ka=0.5, h/a=2, H/L=0.1 and deep water, i , horizontal location of cylinder axis; , solution to first order; , solution to second order. 105 Figure 13. Incident and developed wave profiles to first and second order around a submerged circular cylinder at t/T = 8 for ka = 0.5, h/a = 2, H /L = 0.1 and deep water, i, horizontal location of cylinder axis; , incident wave profiles; , developed wave profiles. Figure 14. Development with time of wave force components to first and second order on a submerged circular cylinder for ka = 0.5, h/a = 2, H /L = 0.1 and deep water. , solution to first order, , solution to second order. 106 ID U) co d q d . (a) i 1 i - i — i i —°——— />••-•* • — "° ^ ^ \ n ^ ^ ^ ^ f t J i / h it fi ii ft i 1 , 1 , 0.0 0.2 0.4 0.6 ka 0.0 0.2 0.4 0.6 ka 0.8 1.0 0.8 1.0 Figure 15. Amplitude of first-order oscillatory force and second-order steady drift force on a submerged circular cylinder as functions of ka in deep water, (a) first-order force amplitude, (b) second-order vertical drift force. Present study: 0,11/3=1.75; 0,h/a = 2. Chaplin (1984): A , h/a = 2. Vada (1987): , h/a = 1.75; , h/a = 2. 107 00 d 0.0 0.2 0.4 0.6 0.8 1.0 ka Figure 16. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of ka in deep water. Present study: Stokes second-order solution, h/a = 1.75, O, h/a = 2; alternative second-order solution, • , h/a = 1.75, • , h/a = 2. Chaplin (1984): A, h/a = 2. Vada (1987): , h/a = 1.75; - - - -, h/a = 2. 108 Figure 17. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of ka in deep water. Present study: Stokes second-order solution, , h/a = 1.91, ,h/a = 2.55; alternative second-order solution,+, h/a = 1.91, x, h/a=2.55. Inoue & Kyozuka (1988): • , h/a = 1.91, • , h/a = 2.55, x component; •, h/a = 1.91, o, h/a = 2.55, z component. 109 Figure 18. Amplitudes of second-order oscillatory force components on a submerged circular cylinder as functions of ka for kd = 2. (a) x component, (b) z component. Stokes second-order solution: , h/a = 1.75, , h/a = 2. Alternative second-order solution: • , h/a = 1.75; • , h/a = 2. 110 h/a Figure 19. Amplitude of second-order oscillatory force on a submerged circular cylinder as a function of h/aforka = 0.35 and deep water. , Stokes second-order solution; • , alternative second-order solution. 1.4 1.8 2.2 2.6 3.0 h/a Figure 20. Amplitudes of second-order oscillatory force components on a submerged circular cylinder as functions of h/a for ka = 0.35 and kd = 2. Stokes second-order solution: , x component; , z component. Alternative second-order solution: • , x component; • , z component. I l l o CS cs I o cs cs I o cs >: cs I o cs <5t o : o o cs I o cs' <5f o ; O q cs -VT = 1 — i > , i i 1 i I . I , 0.0 1.0 2.0 3.0 4. \/l = 2 *7 I . I . 0.0 1.0 2.0 3.0 4. -t/T = 3 y \\ r l I 1 i I . I , 0.0 1.0 2.0 3.0 4. - \ • f/ 1 ' s I i 1 i 1 i . i , 0.0 1.0 2.0 3.0 4. — \> ^ \ 1 ' ' / • r i i 1 i I . I . 0.0 1.0 2.0 3.0 4. x/L Figure 21. Development with time of free surface profiles to first and second order around a semi-submerged circular cylinder for ka = 0.6, H /L = 0.1 and deep water. , solution to first order, , solution to second order. . 112 Figure 22. Development with time of free surface elevations to first and second order adjacent to a semi-submerged circular cylinder for ka=0.6, H/L=0.1 and deep water, (a) upwave surface, (b) downwave surface. , solution to first order; , solution to second order. 113 o O CN O d 0.0 0.4 0.8 1.2 ka m d o d 0.0 0.4 0.8 1.2 ka Figure 23. Runup to second order and amplitudes of free surface elevations at first and second order on a semi-submerged circular cylinder as functions of ka for H / L = 0.1 and deep water, (a) upwave surface, (b) downwave surface. • , |TJ,|; O, T|0; A , |r|2|; • , T| M ( I X (R). 114 o C N O 0.0 1.0 2.0 3.0 4.0 5.0 o <N I ' 1 ' 1 ' 1 ' 1 1 1 0.0 1.0 2.0 3.0 4.0 5.0 t/T Figure 24. Development with time of wave force components to first and second order on a semi-submerged circular cylinder for ka = 0.6, H /L = 0.1 and deep water. , solution to first order; , solution to second order. 115 Figure 25. Amplitudes of first-order oscillatory force components on a semi-submerged circular cylinder as functions ofka in deep water, (a) x component, (b) z component. • , present study; O , Kyozuka's (1980) experiment; , Kyozuka's (1980) theory. 116 o i r i CD L_ CM o d 1 1 • ( b ) 1 1 ! 1 1 1 o O S . -/ / > o o o o ° ° . • . i 0.0 0.5 1.0 1.5 2.0 ka 2.5 Figure 26. Amplitudes of second-order oscillatory force components on a semi-submerged circular cylinder as functions of kain deep water, (a) x component, (b) z component, •.present study; , Miao & Liu's (1989) theory; , Kyozuka's (1980) theory; O, Kyozuka's (1980) experiment. 117 i n o o y I i i i i i i i i . 0.0 0.5 1.0 1.5 2.0 2.5 ka Figure 27. Second-order steady force components on a semi-submerged circular cylinder as functions of ka in deep water. • , x component due to velocity-squared term; A, x component due to free surface fluctuation; • , x component resultant force; • , z component resultant force. 118 ka Figure 28. Amplitudes of second-order oscillatory force components on a semi-submerged circular cylinder as functions of ka in deep water, (a) x component, (b) z component, x, component due to velocity-squared term; • , component due to free surface fluctuation; A, component due to fy; O, component due to (|>2; • , resultant force amplitude. 119 (a) (b) Figure 29. Views of free surface to first order around a surface-piercing circular cylinder for ka = 1 and kd = 1. (a) 6/2TC = 0.85. (b) 6/2TC = 0.35. 120 Figure 30. Views of free surface to second order around a surface-piercing circular cylinder for ka = 1, kd = 1 and H / L = 0.08. (a) 6/2TC = 0.85. (b) 0/27C = 0.35. 121 o O CN o d q T o CN l 1.0 0.8 0.6 0.4 0.2 0.0 0/7T Figure 31. Runup and rundown to second order and amplitudes of free surface elevations at first and second order on a surface-piercing circular cylinder as functions of azimuthal angle, (a) ka = 0.481, kd= 1.332 and H / L = 0.0697. (b) ka = 0.684, kd= 1.894 and H/L = 0.091. (c) ka = 0.917, kd = 2.536 and H/L = 0.1004. Present study: .fnj; , r i 0 ; ,|r|2|; . T L ^ , r\ (R', R). Kriebel (1987): , second-order solution; • , experimental results. 122 123 o 0.5 1.0 1.5 2.0 2.5 3.0 ka Figure 32. Runup to second order and amplitudes of free surface elevations at first and second order at 6 = it on a surface-piercing circular cylinder as functions of ka for d/a = 1 and H/L = 0.1. •, mj; O, Ti0; • , h 2 |; • , ^ (R). 124 o iri o S o E r o CM CT C M O P C M O c5 0.5 1 r A A A. A A A A A ^ A 1.0 1.5 ka 2.0 2.5 (a) 3.0 0.5 1.0 1.5 2.0 2.5 ka 3.0 Figure 33. Maximum amplitude (over 6) of free surface elevation at second order and location of maximum on a surface-piercing circular cylinder as functions of ka for a/d =1. (a) free surface elevation, (b) location of maximum. A, present study; A, K im & Yue (1989). 125 o i n 0.0 1.0 2.0 3.0 4.0 5.0 o 0.0 1.0 2.0 3.0 4.0 5.0 Figure 34. Development with time of wave force and moment to first and second order on a surface-piercing circular cylinder for ka = 1, kd = 1 and H/L = 0.1. , solution to first order; , solution to second order. 126 o o O CD Q. o CS o d 0.5 o o ro ro c^s L o d 0.5 1.0 1.5 2.0 2.5 3.0 ka 1.0 1.5 2.0 2.5 3.0 ka Figure 35. Wave force and moment to second order and components at first and second order on a surface-piercing circular cylinder as functions of ka ford/a = 1 and H / L = 0.1. (a) wave force components, (b) moment components. . F ^ M ^ O , F 0 , M 0 ; A, F 2 , M 2 ; D,F_X, M_ Y ;B, F x , M Y . 127 0.5 1.0 1.5 2.0 ka 2.5 3.0 to cs < cs a CS o c s o o " 0.5 (b) 3.0 Figure 36. Amplitudes of second-order oscillatory force and moment components on a surface-piercing circular cylinder as functions of ka for d/a = 1. (a) wave force components, (b) moment components. Eatock Taylor & Hung (1989): . F j j M ^ l . F j M ® . F ^ M j . Present study: x, tf£ M (£ A, F ? • , F 2 , M 2 . 128
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Time-domain solution for second-order wave diffraction
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Time-domain solution for second-order wave diffraction Cheung, Kwok Fa 1991
pdf
Page Metadata
Item Metadata
Title | Time-domain solution for second-order wave diffraction |
Creator |
Cheung, Kwok Fa |
Publisher | University of British Columbia |
Date Issued | 1991 |
Description | A numerical method based on potential flow theory is developed for simulating transient, second-order interactions of ocean waves with large fixed bodies of arbitrary shape in two and three dimensions. The physical problem is represented by a mathematical model composed of a fluid domain bounded by the body surface, the still water surface, the seabed, and a control surface truncating the infinite fluid region. The nonlinear free surface boundary conditions defined on the instantaneous free surface are expanded about the still water level by a Stokes expansion procedure. The flow potential to second order is thereby defined with respect to a time-independent boundary which includes the still water surface, and its solution involves a numerical discretization of an integral equation. With the potential separated into incident and scattered components, the Sommerfeld radiation condition applied to the scattered potential is modified to incorporate a time-dependent celerity to account for the transient and second-order effects. The free surface boundary conditions and the radiation condition are then satisfied to second order by a numerical integration in time. An alternative second-order solution is derived based on a different expansion procedure in which the nonlinear free surface boundary conditions and an integral equation defined on the instantaneous free surface are both expanded by a Taylor series, and terms up to second order are retained. The two approaches give rise to identical first-order problems, but give rise to second-order problems which are apparently different. The discrepancy arises from the second-order integral equation in which additional second-order terms are retained. The physical interpretations and limitations of these terms are explored and their effects on the evaluations of wave forces are assessed. Applications of the present method are made to studies of regular wave diffraction around a fully submerged and a semi-submerged circular cylinder in two dimensions, and around a bottom-mounted surface-piercing circular cylinder in three dimensions. The stability and numerical accuracy of the proposed solution and the treatment of the radiation condition to second order are examined. Comparisons of computed wave forces and runup are made with previous theoretical and experimental results and these indicate favourable agreement. |
Subject |
Ocean waves -- Analysis |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0050480 |
URI | http://hdl.handle.net/2429/30987 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1991_A1 C44.pdf [ 8.06MB ]
- Metadata
- JSON: 831-1.0050480.json
- JSON-LD: 831-1.0050480-ld.json
- RDF/XML (Pretty): 831-1.0050480-rdf.xml
- RDF/JSON: 831-1.0050480-rdf.json
- Turtle: 831-1.0050480-turtle.txt
- N-Triples: 831-1.0050480-rdf-ntriples.txt
- Original Record: 831-1.0050480-source.json
- Full Text
- 831-1.0050480-fulltext.txt
- Citation
- 831-1.0050480.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0050480/manifest