Technical Report CHL-98-24 September 1998 US Army Corps of Engineers Waterways Experiment Station Coastal Inlets Research Program Directional Irregular Wave Kinematics by Christopher H. Barker, Rodney J. Sobey, University of California at Berkeley Approved For Public Release; Distribution Is Unlimited am aN a wou Prepared for Headquarters, U.S. Army Corps of Engineers Vo. CHI : Under Scour at Inlet Structures Work Unit 32932 1G = 24 The contents of this report are not to be used for advertising, publication, or promotional purposes. Citation of trade names does not constitute an official endorsement or approval of the use of such commercial products. The findings of this report are not to be construed as an official Department of the Army position, unless so desig- nated by other authorized documents. B) vanres ON RECYCLED PAPER MBL/WHO!I A Coastal Inlets Technical Report CHL-98-24 Research Program September 1998 Directional Irregular Wave Kinematics by Christopher H. Barker, Rodney J. Sobey University of California at Berkeley Berkeley, CA 94720 Final report Approved for public release; distribution is unlimited O 0301 0034586 4 Prepared for U.S. Army Corps of Engineers Washington, DC 20314-1000 Under Scour at Inlet Structures Work Unit 32932 Monitored by U.S. Army Engineer Waterways Experiment Station 3909 Halls Ferry Road Vicksburg, MS 39180-6199 [ia US Army Corps of Engineers Waterways Experiment Station — AG fF 8 a eon oy S fE3 ; Reet SS LABORAT! FOR INFORMATION CONTACT: PUBLIC AFFAIRS OFFICE U.S. ARMY ENGINEER WATERWAYS EXPERIMENT STATION 3908 HALLS FERRY ROAD VICKSBURG, MISSISSIPP! 39180-6199 PHONE: (661) 634-2502 AREA OF RESERVATION ¢ 2.7 sqion Waterways Experiment Station Cataloging-in-Publication Data Barker, Christopher H. Directional irregular wave kinematics / by Christopher H. Barker, Rodney J. Sobey ; prepared for U.S. Army Corps of Engineers. 183 p. : ill. ; 28 cm. — (Technical report ; CHL-98-24) Includes bibliographic references. 1. Ocean waves — Measurement. 2. Wind waves — Measurement. 3. Kinematics. |. Sobey, R. J. Il. United States. Army. Corps of Engineers. Ill. U.S. Army Engineer Waterways Experiment Station. IV. Coastal and Hydraulics Laboratory (U.S. Army Engineer Waterways Experiment Station) V. Title. VI. Series: Technical report (U.S. Army Engineer Waterways Experiment Station) ; CHL-98-24. TA7 W34 no.CHL-98-24 Contents 1 Introduction and Review of Current Methods POTS ackpro ui Geeiysmtek sects 2 th Leb eee Ne EN? oe Nis pnt eet a 2) 22 1.2. Methods for the Interpretation of terecular Waveshecords: 4/72). eons anal Ue Mine EC Roe cg Ia tae POG lobaleMethodsrea 2. emetic. Re ce Re 6 rem chen INGE 2 nae ii 2 onleocalelethiod sitet. Saar Rha) WER Pwr i ge eR ae 2 Two Dimensional LFI Theory 2 Problemebormulations si Seen ieerese sees meee eth Mane nte ns 2) hacer PPL) VAD ATA ATESY 2" She A, MERLOT eck Ch as Re A ae: 1 RAR hs td UL ce 2.2 SE inding the Solution) yr epee. 2 5a ee Ere Mee ne earner cn 3 LFI Method for a Point Pressure Record Sale hormulationsonsolutionsiee nay ae nena cea a a gen Sp ehormulationiomtune) Optimization eee men mee se anne 33a iComputationy Methods stare alee? Ayick eee ae ema ne RICE oe d2oe ly Bre-Rrocessing oivecord’ ts | ein tat = beeen Ms Ge Pee cles Sore OPulnaiZationy lar OCeCUne manne eaten eee St el heoretieal Records: 32.88% jou bau t-te Re eee VAL Clroosmae laren Oi Solin o 55002000006 0000 Stone Laboratory se Measurementsmrars ee eee te ee eee ne MOR IDISCUSSTOMA RAR. Te ie nS cbelc VAIL Le Ayah UR 2) a ie 4 Three Dimensional LFI Theory 41 eeThreedimensionaliSeast Ves 2) o.cu. pees en es Pa as Acs ile Backprounid asyrakt, Bybee saa Reese remen: fai an rie oS ot Aa5. 2 9) yar se PS ahh ere Vee ey Ee REAL Seon eee es ce 4.4 A Local Two-Intersecting-Wave Theory................. iii A ANS Kinematics tae eas © cole Lei ce tc Seca 5 Array of Water Surface Traces Heil ie, 5.3 5.4 5.9 5.6 Rojsanplenareray Ch CNONNOMN 5.659 60 66 5 0,6 oa 9 wa oc Hormulationiof the: Optimization lessee nem ane Computations Methodsiu aye eu aisiey men nena 5duly Bre-ProcessinguomhRecords, yeast hacen ten ee 5rot2 ee Optimizations brocedtke sien alla meee TheoreticalSRecordsiys oye = cps ee. Se ms eee Sra oimpleVavesINecordss aie) c)) aemee chen sasme nnn 5.4.2 Records of Two Intersecting Waves ......... aboratony, Measurements ieaen Cara ace once Heo wllbaboratony, 2 xperiin en tien) ue nen ae amen 505.2) sLaboratory, Nesultsie eee ces cases eae oe ASCUSSTOM OCR. (ec!) sy ee oR ana Ee eee canes ee Sree ee ee 6 Array of Subsurface Pressure Measurements 6.1 6.2 6.3 6.4 6.5 6.6 Rormulationvol Solutionmmesc: (ee oe nee waneun: | ay ey eee Hormillahionsofarher: Optimization men nein meine ey one Computation’ Methods yi -. snenone ae ee ee eee Godel, — WARE Pinoresetynayes Ost Invecoitels 4 5 5 oc ao go oo do & 67322 Optimizations brocedure eatin iene ene heoretical Records 4. i520 Me sete) ahh eine eae cere es G4 Al neSingle Wave Recordswe es ws eae Snares aed: 6.4.2 Records of Two Intersecting Waves ......... BielditMeasurements (ia) sca 4) 2), Seine nears Ral Dire aire bree 6:21, GPreldi Results’ 0 24 QUE ne he te Oe Rieee eee nec DISCUSSION: SAWS 5) eo ss Soe ee See ee een 7 Conclusions Cesk Pinte Workin 3 as cee sie 225) ee ren ae en Re Ameen Bibliography A Intersecting wave theory Al iNe2 A.3 A.4 A.5 Introduction’) i) 6 2 kee eg A ee ee rete Theoretical’ Backeround 24 54. eel lo yi een Analytical, Verifications y..00 pen) cues eh enone keene nen Numerical Verification iy) s hacen 6 ah pace ne VAs ee ichardsonsixtrapolationies piel sesi eee een Aca) aMlore'Complexiseasie)= ie) en. oe eee ee Depthsiof Validity 2:5 pcasapeh ath yteuey ween lan He maeime es. [ante are iv Preface This study was funded by the Operations and Maintenance Division, Headquarters, U.S. Army Corps of Engineers (HQUSACE). It is a product of the U.S. Army Engineer Waterways Experiment Station (WES) Coastal and Hydraulics Laboratory's (CHL) Coastal Inlets Research Program (CIRP) under the “Scour at Inlet Structures” Work Unit 32932. Dr. Steven A. Hughes was the CIRP Principal Investigator. Dr. Nicholas C. Kraus was the CIRP Technical Manager, and Mr. E. Clark McNair was the CHL Program Manager for CIRP. Messrs. John Bianco, Charles Chestnutt, and Barry Holliday were HQUSACE Program Monitors. General supervisory direction was provided by Mr. C. E. Chatham, Chief, Navigation and Harbors Division; Mr. Charles C. Calhoun, Jr., Assistant Director, CHL; and Dr. James R. Houston, Director, CHL. This study was conducted by Dr. Christopher H. Barker and Dr. Rodney J. Sobey, Department of Civil and Environmental Engineering, University of California, Berkeley, California, during the period July 1994 to September 1997. The report is a dissertation in partial fulfillment of the requirements for the degree of Doctor of Philosophy awarded to Dr. Barker by the University of California, Berkeley, California. The authors extend their gratitude to the following people for their assistance. Laboratory results presented in Chapter 3 were provided by Mr. Murray Townsend and Dr. John D. Fenton of the Australian Maritime Engineering Cooperative Research Center at Monash University, Melbourne, Australia. Laboratory data presented in Chapter 5 were collected with the assistance of Dr. Steven A. Hughes, Navigation and Harbors Division, CHL, and the technical staff at the Coastal and Hydraulics Laboratory, WES. Field data presented in Chapter 6 were provided by Mr. Gary L. Howell, Prototype Measurement and Analysis Branch, CHL. The FORTRAN code given in Appendix A was based on source code provided by Takumi Ohyama of the Institute of Technology, vi Shimizu Corporation, Tokyo, Japan. Helpful comments and review were provided by Professor Robert L. Wiegel and William C. Webster, University of California, Berkeley, California. At the time of publication of this report, Director of WES was Dr. Robert W. Whalin. The Commander was COL Robin R. Cababa, EN. The contents of this report are not to be used for advertising, publication, or Promotional purposes. Citation of trade names does not constitute an official endorsement or approval of the use of such commercial products. Chapter 1 Introduction and Review of Current Methods 1.1 Background Wave kinematics are involved in most processes that coastal engineers study. Knowledge of fluid velocities and accelerations are necessary for the study of the wave loading of structures through the use of the O’Brien-Morison equation. Knowl- edge of the kinematics near the sea bed are necessary for studies of sediment transport processes. For this reason, it is important that wave measurements be interpreted as accurately as possible. Despite the need for good data on the kinematics of waves, it is impractical to measure every parameter of interest. In a given situation, one might need to know the velocities, accelerations, and pressure fluctuations throughout the depth at a given location. While, in the laboratory, it might be possible to measure these quantities at many locations, it would be very expensive, and it is totally impractical in the field. A pragmatic approach is to measure just a few quantities, and to use a wave theory to compute the balance of the parameters. This approach can give a complete description of the kinematics from only a few measurements. Unfortunately, the accuracy of this approach is limited by the accuracy of the adopted theory, as well as the accuracy of the measurements themselves. A number of wave theories have been adopted in an effort to describe the kine- matics of waves. The most accessible and frequently used of these is Airy wave theory (also known as linear theory or first order Stokes theory). Airy wave theory is easy to use, but has a number of limitations. The source of these limitations is the sim- plifications of the governing equations of gravity waves that are made to linearize the equations, allowing a straightforward solution. These simplifications are made in the free surface boundary conditions and are justified by an assumption of small wave amplitude. Simplifying the free surface boundary conditions reduces the accuracy of the predicted kinematics. Unfortunately, this compromise is at the free surface, which is the location of the greatest fluid velocities and accelerations in waves, and thus frequently the most important to the forcing of structures. The assumption of small amplitude also renders the theory inadequate for large waves, exactly those of greatest interest to coastal engineers. In order to address these limitations, a number of high order steady wave theo- ries have been developed. Commonly in use are Stokes, Cnoidal, and Fourier wave theories. For a review of these, see Fenton (1990). In general, Stokes methods are successful in deep water, Cnoidal methods in shallow water, and Fourier methods in all depths of water. Within their limitations, all three of these methods provide ex- cellent predictions of the kinematics of steady waves, but are not directly applicable to the irregular waves commonly found in the field. With the possible exception of swell conditions on a very mild-sloped bottom, waves in the sea are neither steady, unidirectional nor monochromatic. Methods for the interpretation of measurements of real sea states rely on Airy wave theory even more heavily than do steady wave methods. With modern computer sys- tems, even the mest complicated of high order steady wave theories is quite accessible. However, such higher order theories are not directly applicable to multi-directional or multi-chromatic waves. The linearity of Airy theory allows superposition, in which any combination of waves of different frequencies or directions can be combined to form a solution. Superposition allows real sea states to be easily characterized by frequency and direction spectra. While accessible, this method results in solutions for the kinematics of the waves that do not satisfy the full free surface boundary "wave" Figure 1.1: Segment of record considered for global approximations conditions, and thus do not take into account the interactions between the individual waves. 1.2 Methods for the Interpretation of Irregular Wave Records Methods used for the interpretation of irregular wave records fall into two general categories: global and local approximations. Global methods seek a solution that matches an entire measured record, or a single complete measured wave, from trough to following trough, or zero crossing to zero crossing (Fig. 1.1). These methods apply the same frequency and wave number (or set of frequencies and wave numbers) for all z (vertical variation) and t (time). Local methods, on the other hand, seek an approximation to each small local segment of a measured wave. In these methods, the frequency and wave number still apply for all z, but are allowed to vary with time, providing a separate solution in each small window in time (Fig. 1.2). Local Segment Figure 1.2: Segment of record considered for local approximations 1.2.1 Global Methods Spectral Methods The most commonly used global method for the analysis of irregular waves is spectral analysis, coupled with superposition of linear waves. Spectral methods based on a single point measurement seek to define the sea state with a variance spectrum, E(w) which specifies the contribution to the variance of the water surface at each frequency, w. In practice, the spectrum is defined in discrete form, with the water surface described by the superposition of many linear waves, M aot) = Se Am COS (Km (x CoS Am + ysin On) — Wmt + Am) (ali) m=1 where 77 is the elevation of the water surface, w,, and k,, are the frequency and wave number of the mth wave, 0,, is the direction of propagation of the mth wave and a,, is the phase of the mth wave at the origin of the coordinate system. The amplitudes are found from the discrete variance spectrum, Om = V 2E (wm )Aw (LZ) where Aw is the frequency spacing of the discrete spectrum. The frequency and wave number of each of the M waves are related by the linear dispersion relation, w? = gkm tanh kph (1.3) where g is the acceleration of gravity, and h is the mean water depth. The dispersion relation defines the phase speed of each component, allowing each separate component to move independently, without interaction with one another. Any Eulerian current is not included in this relation, and is generally ignored in spectral methods. The discrete variance spectrum, EF’, can be computed from a water surface or sub- surface pressure record through the use of variations of the Fast Fourier Transform (Dean and Dalrymple 1991; Newland 1993). A point measurement provides no in- formation about the direction of propagation of the waves, 6,,. It is usually assumed that all the waves are propagating in the same direction. Once the amplitude, fre- quency, wave number and phase of each individual wave are defined, the kinematics and dynamics of the wave field can be computed by superimposing the kinematics of each individual wave, as predicted by linear wave theory. Directional Spectra When measurements are taken by an array of instruments, the directional nature of the sea can be described by a directional spectrum, S(w, @). In this case, the water surface is represented by a large number of linear waves of different frequencies and directions: M N (2508) = es Ss Amn COS (km(x cos, + y sin O,) — Wmt + Amn) (1.4) m=1 n=1 where: Onn N29) Com a NONG (1.5) and Aw and A@ are the sample spacing of the discrete spectrum in frequency and direction space. w,, and k,, are once again related by the linear dispersion relation (Eq. 1.3). The directional spectrum is usually broken down into two parts, the one dimensional variance spectrum, E(w), and the directional spreading function (DSF), D(w, 6): S(w, 0) = E(w)D(w, 8) (1.6) where E(w) is the variance spectrum defined above, and ‘3 iP D(w,0)d0dw = 1 (1.7) There is a great deal of literature about how to best determine the DSF from a variety of arrangements of measurements. Current practice is reviewed in Mansard (1997). Unfortunately, all of these methods are limited in their ability to define accurately the DSF. The time series measurements used to determine the variance spectrum commonly include in excess of 1000 points in time. This much data in the time domain allows a very high resolution computation of the spectrum of a stationary process in the frequency domain. In the spatial domain, by contrast, there are only as many data points as there are instruments in the particular measurement array. This may be as few as three, and perhaps as many as a dozen, but it is too expensive to use many more than that. As a result, there is little information to define the DSF. For example, when data from an array of three instruments is analyzed by the standard linear method, the DSF can be specified by only five independent coefficients (Dean and Dalrymple 1991). While these few coefficients may serve well for determining integral properties, such as the mean direction and radiation stress, it is not enough information to accurately specify the complete kinematics. Another difficulty arises when determining the spectrum from measurements other than wave staffs. Most often, subsurface pressure gauges or a combination of pressure gauges and orthogonal velocity gauges are used. In this case, the measured quantity must be related by a transfer function to the equivalent water surface. The transfer function is most commonly determined from linear wave theory. This use of linear wave theory once again contributes to errors in situations where linear theory is not entirely appropriate. Bishop and Donelan (1987) discuss the difficulties in determin- ing the transfer function from a subsurface pressure measurement to the equivalent water surface. The commonly used methods for determining the DSF are statistical methods, in that they rely on the assumption that the phases of the individual components are randomly distributed. This assumption allows the DSF to be computed by dis- regarding the phase information. Unfortunately, without the phase information, it is impossible to reconstruct the detailed kinematics, only statistical descriptions can be formulated. Even if a frequency or frequency-direction spectrum has been accurately deter- mined, these methods have a number of shortcomings when used to predict the kine- matics of irregular waves. Shortcomings include the inaccuracies inherent in linear wave theory, particularly in shallow water and with large waves. An additional difficulty arises from the superposition of many waves. This problem is known as high frequency contamination of the crest kinematics (Forristall 1985; Lo — and Dean 1986; Bishop and Donelan 1987). Fundamentally, the difficulties arise from the approximations made by linear wave theory in the free surface boundary conditions. In linear theory the free surface boundary conditions are applied at the mean water level, and thus predictions made above that level are strictly out of the solution domain. If the full free surface boundary conditions are not satisfied, the resulting predictions will be inaccurate, particularly near the free surface. In particular, the hyperbolic function quotients that define the vertical variation of the kinematics become very large in the region above the MWL for the high frequency (and high wave number) components. This results in substantial high frequency fluctuations in the predicted kinematics near the crest. In an attempt to reduce the inaccuracies in linear superposition’s predictions of the near surface kinematics, empirical modifications to Airy theory have been adopted (Wheeler 1969; Lo and Dean 1986). This method, known as Wheeler stretching, locally adjusts the vertical dimension to prevent the evaluation of the hyperbolic quotients from being evaluated above the MWL. The result is a hybrid global-local method in that the frequency, wave numbers, and amplitude of each wave are deter- mined globally, but the coordinate system is defined locally, varying with time. The stretching method produces predictions of crest kinematics that seem to match measured data better than simple linear superposition, but it no longer satis- fies either the Laplace equation (mass conservation), or the full free surface boundary conditions. Another attempt to improve on the accuracy of determining the kinematics from directional spectra is the recent work by Prislin et al. (1997), and Prislin and Zhang (1997). This method seeks to reconstruct wave kinematics from a directional spec- trum using a second order Stokes-type interacting wave theory. The spectrum is decomposed into a set of individual free waves, with from one to five directional free waves per frequency. The effects of the second order interactions are subtracted from the measured record to determine the amplitudes and phases of the individual free waves through an iterative procedure. This results in an expression for the potential function and water surface that includes a full set of many free waves, and the corre- sponding second order bound modes. The method succeeds in globally reproducing the measured kinematics in the given deep water field records quite well, although the largest errors are in the vicinity of the crest, where velocities are greatest, and accuracy is most important. Another limitation of the Prislin and Zhang method is that the nonlinear interac- tions are computed through the use of a Stokes-type perturbation expansion in wave steepness. Based on experience with steady waves, a second order expansion of this type is likely to be adequate in deep water, but if the method were to be applied in transitional, and especially in shallow water, a much higher order representation would be necessary (Fenton 1990). While it is theoretically possible to extend the method in this manner, it would increase the magnitude of the computation substan- tially, perhaps prohibitively. Representing an entire record in the global sense requires many interacting free waves. Considering their interactions at high order would pro- duce huge numbers of interacting terms that must be considered. An alternative, local, approach would need to consider far fewer interacting waves. Other global methods rely on zero crossing analysis to identify particular waves that are then analyzed by using steady wave theory for a wave of the same height and period. This approach can provide an order of magnitude estimate for the kinematics, but does not take into account the detail of the record, and thus can not be expected to consistently provide better than order of magnitude accuracy. To include the detail of the record in wave by wave analysis, Dean (1965) adapted his numerical stream function method to irregular waves, seeking a Fourier expansion for the stream function that includes both cos j(kx — wt) and sin j(kx — wt) terms. The method optimizes the Fourier amplitudes to best solve the free surface boundary conditions at the measured water surface elevations from trough to following trough. While this approach takes into account the detail of the record, the method includes z = (water surface) Figure 1.3: Coordinate system for Lambrakos-Baldock-Swan method. only a single free mode, with all other included components being bound modes traveling with the wave at the same phase speed, thus representing an asymmetric wave of permanent form. A solution that does not allow a change in form is unlikely to accurately capture waves in deep water, where frequency dispersion can lead to transient extreme waves, or, indeed, in shallow water where shoaling effects cause a change in form as the waves progress. Seeking to improve on global methods, Lambrakos (1981) developed a method for determining the kinematics of two dimensional irregular waves that includes many free modes, and thus unsteady motion. Baldock and Swan later refined Lambrakos’ method, applying it specifically to large transient waves, first in deep water (Baldock and Swan 1994), and then in shallow water (Baldock and Swan 1996). Baldock and Swan’s method adopts the following potential function that is a double Fourier expansion in space and time: N M NG Ft) = yy ee cosh (nkz) (Anim cos (nkzx — mwt) + Brim sin (nkz — muwt)) n=1 m=1 (1.8) where k, w, Anm, Brym are global constants, and z = 0 at the bed (Fig. 1.3). In this potential function each frequency component (mw) has a corresponding set of wavelengths (nk), allowing each component to travel at distinct phase speeds. 10 This approach can accommodate both steady and unsteady wave forms by including both a set of free waves and a corresponding set of bound modes. Periodicity over time and space is assumed. The record is not likely to be periodic, so the duration of the record is assigned as the fundamental period, (27/w), and the fundamental wavelength, (27/k) is found as part of the solution. The method uses a nonlinear optimization to find the coefficients, A,,., and By, m, the fundamental wavenumber, k, and the Bernoulli constant that produce a solution that globally satisfies the full free surface boundary conditions. In order to accommodate the unsteadiness of the wave profile, the evolution of the wave in space must be known. This is accomplished by applying the free surface boundary conditions at many locations in time and space. Since the elevation of the water surface usually is measured only at a single spatial location, it is predicted at the other locations as part of the solution. The potential function (Eq. 1.8) exactly satisfies the bottom boundary condition and mass conservation. In order to find the solution that best fits the free surface boundary conditions, the method seeks a set of coefficients that minimizes the sum of squares error in both of the free surface boundary conditions over a grid of nodes in time and space. Baldock and Swan identified a difficulty in this basic method, as the squared error was equally considered at all of the nodes, most of which were at locations in which the water surface elevation was also unknown. This resulted in solutions in which the error in the boundary conditions was greatest at the location of the measurement. The difficulty was mitigated by a weighting function, multiplying the dynamic free surface boundary condition errors at the measured location by a factor of 50 in the sum of squares calculation, to force the solution to match well at that point. Comparisons of their results with laboratory data were quite good, but the method has a number of limitations. These include a huge matrix of unknown coefficients that must be found simultaneously (2(/7N +1) unknowns, with typical values of M and N of 18, for a total of 650 unknowns), and the method must solve for the water surface far from the measurement location, removing the focus from the actual measured data. The method makes no assumptions about the steadiness of the wave field; however, it extends the horizontal bottom assumption and introduces unnecessary complication 11 in its attempt to capture an entire segment of an unsteady wave field with a single expression. This complication would become at least an order of magnitude greater if the method were extended to include the second horizontal dimension. 1.2.2 Local Methods Nielsen (1986, 1989) introduced a method that uses a local frequency and linear wave theory to find the location of the water surface from a pressure record. In this method, the wave in the region of each measured point is considered to be a small segment of a linear wave, so that the pressure record is locally represented by a sine function. p = An sin (Wptn — Qn) (1.9) so that: 1 02p n= 4{-—-— 1.10 Bs p Ot? Gey) Estimating the curvature by finite differences, the local frequency becomes: —Pn-1 ate 2Dn — Pn+l1 Sap pi Se ae ial = pAb? ay) where At is the time spacing between points. Alternatively, w can be computed from trigonometric identities, On = 500s" (Batre | (1.12) which is exact for a sinusoid. Once the frequency has been determined, the wave number is computed from the linear dispersion relation, Eq. 1.3. The water surface elevation can then be computed with the linear pressure response function: _ Pn cosh(knh) = 113 pg cosh(ky Zp) Was) Nn where z, is the elevation of the pressure measurement above the bed, and h is the mean water depth (As the method is a local method for the interpretation of a point 12 measurement, the local mean water depth is used, rather than the global still water depth). This may be adapted somewhat by using a variation of Wheeler stretching: cosh (k,, {h + 2 Tn = poe) (1.14) Nielsen also presented an alternative approach, which uses a semi-empirical trans- fer function derived from Fourier wave theory to compute the water surface elevation from the measured pressure and local curvature of the pressure record. While quite accurate at reproducing the water surface, neither method supplies the kinematics, nor do they satisfy either mass conservation and or the free surface boundary con- ditions. Despite these limitations, the efficacy of these methods demonstrates the potential for the local approach. To interpret bottom pressure measurements in the context of the kinematics while preserving the full governing equations, Fenton (1986) presented a method that em- ploys a local polynomial approximation to the complex potential function. In Fenton’s method, the potential function and water surface are represented by separate poly- nomials in each small window in time: M o(x — ct,y) + iv(ae — ct,y) = iene a (1.15) M Nee) = SPE — ct) (1.16) j=0 where z = x+y, y = 0 at the bed, 77 is the water surface, c is the wave celerity, and the a; and 6; are real. The wave is assumed to propagate at speed, c, without change in form. While steadiness is not a valid assumption in the global sense, it is only applied locally, within a small window in time. The method solves for the coefficients, a; and b;, and the wave celerity, c, that satisfy the full nonlinear free surface boundary conditions and fit the measured pres- sure record. This approach provides the complete kinematics and satisfies the full governing equations. Based on a polynomial variation with depth, it works well in 13 shallow water, where Cnoidal methods with polynomial vertical variation are theo- retically appropriate. At the same time, it may not be applicable in transitional or deep water, where Stokes type methods, based on an exponential variation in the vertical, are theoretically more appropriate. In a later paper, Fenton and Christian (1989) presented a simplified version of the method that required less algebraic ma- nipulation, and was still found to be effective in shallow water. In this case, the local wave celerity was assumed to be that given by long wave theory: C= SG (L107) While this is a reasonable approximation for long waves in shallow water, it was not appropriate for shorter waves in deeper water, as might be expected from the polynomial form and the long wave celerity. To find a method that could work in any depth of water, Sobey (1992) developed the Local Fourier Method for Irregular waves (LFI). This approach employs a po- tential function represented by a low order Fourier expansion in a small window in time. It is a method derived for the analysis of a point water surface trace. Local frequency, wave number, and the Fourier coefficients are sought that fit the measured record and the full free surface boundary conditions. The LFI method provides the ' complete kinematics, satisfies the full governing equations, and is successful in all depths of water. A more complete description of this method follows in Chapter 2. Sobey’s method shows a great deal of promise, but was only applied to the analysis of a water surface measurement at a single location. Pressure sensors are easy to deploy, and thus are often used to measure waves, particularly in shallow and transitional depth water. The LFI method is extended in this dissertation to the interpretation of wave measurements from a point pressure sensor in Chapter 3. Point measurements provide no information about the directionality of the mea- sured wave field. Wind seas are not uni-directional, and knowledge of directionality has been shown to be very important in the prediction of the kinematics and forcing of structures by real seas (Dean 1977; Forristall et al. 1978). Arrays of measurements are frequently used in order to capture the directional 14 nature of the wave field. The LFI method has been further extended in this disserta- tion to the interpretation of arrays of water surface measurements in Chapter 5, and arrays of pressure sensors in Chapter 6. 15 Chapter 2 Two Dimensional LFI Theory It is common to deploy a single measurement device to gather information about the wave field at a given location. Measurement at a single point in space does not give any information about the directionality of the sea state, but if the wave field is assumed to be locally two dimensional, a reasonable description of the kinematics can be established. This chapter outlines the LFI method as it is applied to the interpretation of a time series from measurements taken at a single location in space, such as recorded by a pressure gauge or wave staff. 2.1 Problem Formulation The problem formulation for two dimensional irregular waves has much in common with classical steady wave theory. The flow is assumed to be incompressible and irrotational. The kinematics can therefore be represented by a potential function, o(x,z,t), where _ 0 Ox _ 96 U and u and w are the horizontal and vertical velocities, respectively. 16 | Nee surface (7) BU) 10) sea bed Ue Figure 2.1: Coordinate system used for the two dimensional method Field Equation The field equation is mass conservation for irrotational flow, in the form of the Laplace equation: Od 0d Wie = ey ¢ Ox e Oz? : 2) Boundary Conditions The boundary conditions are the bottom boundary condi- tion for a locally horizontal bed (BBC), eee = apie 0 at a) (2.3) (20) the kinematic free surface boundary condition (KFSBC), w—-——-—u— =0 at z=N (2.4) 1 rt —+=-uv+-w?+9n-B= at z= (2.5) where 77 is the elevation of the free surface, g is the acceleration of gravity, and B is the Bernoulli constant. The kinematic free surface boundary condition (Eq. 2.4) includes both the time and spatial gradients of the water surface. When working with a subsurface gauge, 17 the location of the water surface is unknown and part of the solution. When working with data from a water surface probe, the location of the water surface at discrete points in time is known, but the gradients in time and space are not. In order to accommodate this lack of data, the kinematic free surface boundary condition is transformed to eliminate the gradients of the water surface. Following Longuet- Higgins (1962), the kinematic condition is subtracted from the the dynamic condition differentiated following the motion: MKFSBC = —g(KFSBC) + = (DFSBO) (2.6) Resulting in a modified kinematic free surface boundary condition: qP Uap Uo = = (0) at 227) This new form of the boundary condition does not include the gradients of the water surface. Applying both the modified kinematic free surface boundary condition and the dynamic free surface boundary condition completes the formulation. Observational Equations In steady wave theory, periodic boundary conditions are also imposed, forcing the solution to be periodic in both space and time. For irregular waves, however, the periodicity is not known. Rather, the local solution is defined by a local segment of a measured record within a small window in time, together with the field equation and the bottom and free surface boundary conditions. In order to define a solution that fits the measured record, observational equations are identified. These equations will be different depending on which quantity has been measured. In the case of a water surface measurement, they are the free surface boundary conditions, applied at the measured elevations at a number of points in time throughout the window (Sobey 1992). In the case of a pressure measurement, they are the Bernoulli equation, applied at the elevation of the measurement, and also at a number of points in time within the window (see Chapter 3). 18 Solution in Each Local Window A form for the potential function in each local window is based on Fourier steady wave theory. This is the potential function used by Sobey (1992): @(z,2,t) =Ux + yA jE sin i(k — wt) (am) where U and hf are the known cera Eulerian current and mean water depth (see section 3.3.1 for a discussion of these important parameters), J is the truncation order of the Fourier series, A; are the Fourier coefficients, and w and k are the local fundamental frequency and wave number. The above potential function exactly satisfies mass conservation and the BBC. This form for the potential function is periodic in space and time, however the periodicities are not defined apriori, but found to fit the record, defining a local frequency and wave number. The measured record is broken down into individual segments, each in a separate window in time. In each window, a different set of the parameters w, k, kx, and A; are found to fit the segment of the measured record. This represents that segment as a piece of a larger, periodic wave. The entire record is then represented by separate potential functions, each applied to a particular window in time. 2.1.1 Dynamics The potential function provides the complete kinematics, and the dynamics are found through the unsteady Bernoulli equation: Ft 5 +u?) +a24+2 -B=0 (2.8) where p is the mass density of the water, and p is total pressure. Total pressure below the water surface can be broken down into three components. Atmospheric pressure (pa) is the pressure of the atmosphere at the water surface, hydrostatic pressure (p,) is the pressure in the water column due to gravity, in the absence of motion. Dynamic pressure (pq) is the component due to the motion of the fluid. Separating these three components focuses attention on the wave motion. P= Pa + Pa + Ph where Ph = —pgz jor Zhe (2.9) 19 When atmospheric pressure is defined as zero, and dynamic pressure is substituted into Eq. 2.8, the result is: xo) 1 2 2 Pd — = = ) — — == » ag oC ee at (2.10) The Bernoulli Constant In a potential flow, the Bernoulli constant is the same throughout space and time. For the special case of wave motion, the value of the Bernoulli constant can be computed if the kinematics are known (Longuet-Higgins 1975). The Bernoulli equation (2.8) and the bottom boundary condition are applied at the bottom: set stot Bao (2.111) where wu,, pp, and z, are the velocity, pressure, and elevation at the sea bed. If the flow is periodic, a time average over a period results in: T= De a pu toa t > -B=0 (BI) where the over-bars indicate time averaging. In the case of steady periodic wave motion, the total vertical momentum at any horizontal location must be the same at the beginning and end of a period. In order for this to be the case, the vertical momentum averaged over a period must be a constant. In order for the momentum to remain constant, the time averaged net vertical force on the water column at that point must be zero. If the pressure at the water surface (p,) is taken to be zero, then the force of the pressure at the bed must be equal to the gravitational force on the column of water, so that the mean pressure on the bed is hydrostatic, Ps = pg(7 — 2s) (2.13) resulting in a simple and exact expression for the Bernoulli constant. ast i B= gi + su; (2.14) 20 With z defined to be zero at the mean water level (7) and Eq. 2.7 as the potential function, B becomes (Sobey 1992): 2 =i mele hae gkA; 5 i Se su a7 > (5) (2.15) and all the terms in Eq. 2.10 are defined by the potential function, allowing the dynamic pressure to be computed from the kinematics. 2.2 Finding the Solution The bulk of the LFI method is the process of determining appropriate values for each of the parameters, w, k, kx, and A,, in the potential function (Eq. 2.7) for each window in time. The result is a set of potential functions that completely describe the kinematics of the wave field in the region of the measurement, each separate window in time being described by a unique and independent potential function. In the case of subsurface measurements, the location of the water surface is also required. Determination of the parameters of the potential function, as well as the location of the water surface, is accomplished through nonlinear optimization that matches the measured record and minimizes the error in both free surface boundary conditions. Sobey (1992) applied the LFI method to the interpretation of a single water surface trace. The following chapter applies the method to a subsurface pressure record. Chapter 3 LFI Method for a Point Pressure Record Pressure gauges are relatively inexpensive and easy to deploy, and as a result are frequently used in the field to measure waves, particularly in shallow water. The simplest instrument consists of a single subsurface pressure gauge. While a single point measurement provides no directional information, a reasonable estimate for the kinematics of the measured waves can be obtained. This chapter describes a technique _ for applying the two dimensional LFI method presented in the previous chapter to the analysis of a point subsurface pressure trace. 3.1 Formulation of Solution The flow is taken to be two dimensional, with the governing equations described in Chapter 2. These include the potential function, hjk(h +z WG, 21) = er ee 2 sin j(kx — wt) (3.1) the modified kinematic free surface boundary condition (f*), 22 Od Ow Fg oad KD Sic) (eee aes jf = Fp ah ED A Pa) an oe AL Ow dr Oz Ou Ow =F a + w= =0 at B= i (3.2) the dynamic free surface boundary condition (f”), Oo = D 2 2 me puts i = Fe se ath Gj Ba) 2= 7 (3.3) and the unsteady Bernoulli equation ( f8). 0d 1 (de BOON ea o BI ipa Hels rams if Sar al ae B= 0 (3.4) with the Bernoulli constant defined as: fh, 1 1 AlnAle \* B = aCe oe 33, a “ 4 X (5) ce) The unknown parameter, x, appears in the potential function only when coupled with the parameter, k. It is convenient to solve for the non-dimensional parameter, kx, essentially a phase parameter in the potential function. The action of waves is greatest near the surface, so that the fluid velocities and accelerations are largest near the surface and decay rapidly with depth, particularly in deep water. In addition, for pressure sensors placed at or near the bottom, the vertical component of the velocity is damped out completely, as expressed by the bottom boundary condition. The problem, then, is to extrapolate the motion of the fluid up to the surface, where the motion is greatest, using only data extracted from the subsurface, where the motion is much less. 23 Mathematically, this is an ill posed problem, as the flow is governed by an elliptic equation (Laplace, Eq. 2.2), in which the solution is determined by the boundaries, but there are only data at a single location. Some of these difficulties can be sur- mounted through the application of the free surface boundary conditions. While there are no data measured near the surface, the free surface boundary conditions remain appropriate, and necessary to define the solution at that boundary. When the LFI method was applied to a water surface trace (Sobey 1992), the location of the the water surface was known, and the boundary conditions could be directly applied at that location. Unfortunately, when working from a subsurface pressure record, the location of the water surface is not known. While the free surface boundary conditions are well defined, the fact that the location at which they must be applied is not known makes the problem more difficult. This is the complication that leads to the difficulties in finding full nonlinear solutions to all free surface flow problems. In order to apply the free surface boundary conditions when working with a subsurface record, the location of the water surface must be found, together with the potential function, in each window. In order to locate the free surface, the water surface is defined at N surface nodes equally spaced in time throughout the window. The elevation of these nodes “is unknown, and will be sought as part of the solution. Including the water surface nodes as part of the sought solution introduces N additional unknown parameters for a total of 3+ J+ N unknowns in each window (k, kz, w, Ay... Ay, ---7n). The free surface boundary conditions, Eq. 3.3 and 3.2, are applied at each surface node. 3.2 Formulation of the Optimization Finding the unknown parameters in a nonlinear system of algebraic equations is known as nonlinear optimization. The system in this case consists of the nonlinear Bernoulli equation and the nonlinear free surface boundary conditions. The free surface boundary conditions are nonlinear in two respects, involving second order terms in dynamic free surface boundary condition (the u” terms and B in Eq. 3.3), and second and third order terms in the modified kinematic boundary condition (the 24 ae ~~ nodes are sought <——— Bernoulli equation applied a ae —- a Zi at known elevation, zp | al A i a tl si ma mee a > _____ Elevations of water surface 2 ABE fP Figure 3.1: Schematic of system of equations in a window uz and u?2 terms in Eq. 3.2). There is also a significant nonlinearity introduced b Ot Or q g y vy the application of the boundary conditions at the unknown and varying free surface. Observational Equations The given form for the potential function could repre- sent any periodic flow, subject to the bottom boundary condition. The FSBCs define the flow as a gravity-constrained, free surface flow. The observational equations are the equations in the system that force the solution to fit the given record. For a subsurface pressure record, this is the Bernoulli equation, applied at the location of the pressure measurement. The required number of independent equations are es- tablished by applying the Bernoulli equation at a number of times throughout the window considered. The error in the Bernoulli equation is the difference between the measured dynamic pressure and that computed from the kinematics defined by the potential function. The solution is the set of parameters in the potential function and the set of water surface nodes that produces a predicted dynamic pressure that matches the measured record, while simultaneously satisfying the FSBCs. A system of equations is specified if there are as many independent equations as unknown parameters in the system. If there are more equations than unknowns, the solution can be defined as that which results in the smallest squared errors in the equations. This least squares formulation is also appropriate for a uniquely defined 25 system, with an expected error of zero in all the equations. In order to specify the solution to the LF] formulation, the boundary conditions (f? and f*) are applied at each of the water surface nodes, yielding 2V equations, and the Bernoulli equation (f®) is applied at I nodes on the pressure record within the local window (Fig. 3.1). There are a total of 3+ J+ N unknowns (k, kz, w, Ai...Ajz, and m...n). The system is uniquely specified for / + N =3+ J and overspecified for /+ N >3+J, resulting in the following least squares optimization: minimize O(X )= SAP 2G has ten Pat CS Tins t tn)? =P So F(X: lie Bop ta) (3.6) 7=1 XS = (6, the A ooo Alig ih ooo iy) where t, are the locations in time where the water surface nodes, 7,, are sought. t; are the times within the window where the Bernoulli equation is applied, Py; are the measured dynamic pressures at t;, and z, is the elevation of the pressure gauge. Overspecification, particularly with additional nodes on the pressure record, is advantageous for an actual record as a way to minimize the effect of any noise in the measurements. Additional nodes on the pressure record also serve to emphasize the measured data, which directly define the local kinematics. The details of the solution to this system of equations is given in the following section. 3.3 Computation Methods The LFI method can be broken down into the following sequence of steps: 1. Pre-processing of record. (a) Determine estimate for level of noise in the record. (b) Determine MWL and subtract hydrostatic pressure from record. (c) Determine estimate for magnitude and direction of Eulerian current. 26 (d) Define a continuous record from the discrete observations by cubic spline interpolation. (e) Specify spacing of output locations. (f) Compute mean zero crossing frequency. (g) Non-dimensionalize record and all parameters. 2. Primary values of numerical solution parameters are chosen. (a) Window width (70) (b) Order of solution (J) (c) Number of sample locations on record within each window (J) (d) Number of water surface nodes within each window (JN) 3. For each selected output location, a window in the record is defined, and an LFI solution is computed. (a) Initial guess for the optimization is computed. (b) Full nonlinear optimization for all unknown components of the potential function is computed. (c) Results are checked for spurious solution. i. If no solution, or a spurious solution, is found the solution parameters are adjusted, and the optimization repeated. ii. If a good solution is found, progress to the next window. 3.3.1 Pre-Processing of Record Accommodating Measurement Error Measurement error can be a major source of difficulty, as the method relies on the details of a small segment of the record. In Sobey’s work (1992), the primitive kinematic free surface boundary condition (Eq. 2.4) was applied at the measured water surface, requiring an estimate for the time gradient of the surface. The estimation 27 of the gradients from a measured record is very sensitive to error in the record, and thus Sobey found it necessary to smooth field and laboratory measurements with a simple moving average filter. There are no pressure gradient terms in the Bernoulli equation, so the application of the LFI method to a measured pressure trace should be less sensitive to noise. Nevertheless, measurement error remains a problem in applying the LFI method toa pressure record. When there is obvious noise in the record, an alternative to smooth- ing is to substantially overspecify the system of equations. By taking many samples on the pressure record in each window, the least squares optimization may accommo- date the errors in the measurements. This is preferable to smoothing, as no smoothing algorithm can reproduce lost data, but will rather impose some apriori assumption about the nature of the record on the results. Accommodating the measurement error with many samples on the pressure record allows the least squares optimization to find the best fit without imposing any further assumptions on the results. This was the approach taken with the presented laboratory data (Section 3.5). Records gener- ated analytically by Fourier wave theory contained no error, and therefore required no special accommodation. - The Mean Water Level The mean water depth, h, is included in the potential function (Eq. 3.1). This is an unknown value when the only data available is a pressure record. While an approximate depth of water is likely to be known at a given deployment site, the water level fluctuates due to astronomical and storm tides. These processes happen over a much larger time scale than the periods of wind generated waves. An estimate for the local mean water level (MWL) can be determined by averaging the measured pressure (after subtracting out the atmospheric pressure) over a period of time that must be larger than a typical wave period, but smaller than the period of tidal fluctuations, to accommodate changes in the mean water level. This mean pressure is the local hydrostatic pressure, which is used to compute the mean water level, and is subtracted from the pressure record to define the dynamic pressure. While this is 28 not an exact procedure, the MWL is not a critical parameter, as h is only used to locate the origin of coordinates in the potential function. Note that the result is a local MWL, which could be quite different than the global still water level (SWL). The mean water level is the more appropriate reference frame, as the method is local, and only considers a small window in time at once. The elevations of the water surface nodes are found independently, without any requirement that the resulting mean water level be at the origin. Thus it is not necessary that the origin is at the exact MWL, only that the it is near the surface. In shallow water, the bed could be used as the origin (Fenton 1986), but this would not work well in deeper water, as the area of primary interest, the water surface, would be far from the origin. The Eulerian Current The depth uniform Eulerian current, U, appears as a known parameter in the potential function, Eq. 3.1. Unfortunately, a measured pressure record gives no in- dication of the local current. Unlike the mean water level, the Eulerian current is a critical parameter that defines the propagation medium, and varying the value used for the current will have a considerable effect on the results of the computation. In fact, as pointed out by Fenton (1985a), any attempt to apply a wave theory without knowledge of the local current will be inaccurate, even at only first order. The situation is not hopeless, however. The presence of the Eulerian current term in the potential function draws attention to the fact that it is an important parameter, even if it is taken to be zero. This should encourage investigators planning field experiments to establish a method of estimating the local current. This may be as straightforward as placing a current meter nearby, or as simple as using tidal data to estimate the local current. Caution must be taken with the later method, as the accuracy of the estimate may not be consistent with the use of high order interpretation methods. It is also assumed in the LFI method that local current is depth uniform. Although this is unlikely to be strictly true, Kishida and Sobey (1988) demonstrated that for steady waves, a current profile with a realistic level of vorticity was likely to yield very 29 similar results as a depth uniform current profile for the wave generated kinematics, even at high order. It is expected that these results are applicable to irregular waves as well, and thus the assumption of depth uniform current is reasonably appropriate. Another complication to be considered is the direction of any local current with respect to the propagation direction of the waves. The U in Eq. 3.1 is the component of the local current in the direction of wave propagation. In the case of a single point measurement, no information is available to indicate the wave direction. Even if the current magnitude and direction are known exactly, without an estimate for the propagation direction of the waves, an accurate estimate for U is not available. Once again, further information could be used to help resolve this difficulty. Knowledge of the local bathymetry and the wave climate could provide this information. For example, waves near a beach tend to propagate almost perpendicular to the shore, or if there are directional measurements taken nearby available, they could be combined with refraction computations to provide an estimate of the local propagation direction. There must be some information available about both the local Eulerian current and the local wave propagation direction in order to interpret accurately a point measurement. The pressure records used for the examples in this chapter were generated either analytically or in a laboratory flume. In both cases, the Eulerian current and wave propagation direction are known. Spline of Record In order to avoid being restricted by the sampling rate of a particular set of measurements, a continuous record is computed by spline interpolation among the measured points. This allows the window widths and the number and location of samples in each window to be chosen independently of the sampling frequency of the data. Care must be taken, however, to be sure to include an adequate segment of the record in each window. A window that includes only a single measured point is computationally possible, but would result in misleading accuracy in the resultant solution. 30 Output Locations The spacing of the desired output locations must be chosen to determine the placement of the windows. The solution in each window is expected to be most accurate in the center, and thus a separate window is centered at each location where output is desired. As each window is computed independently, there is no restriction on the spacing of the output windows. Non-Dimensionalization The comparisons of the errors in equations of different dimensional quantities may result in spurious solutions that over emphasize certain equations in the formulation. While the familiar dimensional forms of the equations have been presented here for clarity, all parameters and variables are non-dimensionalized by physically identifiable parameters before computation. The mass density of water (p), acceleration of gravity (g), and the mean zero crossing frequency of the measured record (w-) are used to define characteristic length, time and mass scales. Length scale = g/w? Time scale = 1/w, (3.7) 3 Mass scale = 7 5) Zz 3.3.2 Optimization Procedure The primary process in the LFI method is a nonlinear optimization of a system typically involving up to 14 algebraic equations in 12 unknown parameters. A system as complex as this may have numerous local minima that result in spurious solutions. The best way to avoid these solutions, and to ensure efficient optimization, is to start with a good estimate for the unknowns in the system, and, in addition, to have a set of criteria for identifying spurious solutions. 31 Initial Estimate The first step in each window is to establish an initial estimate for the optimization procedure. Linear wave theory can be used to produce estimates for the parameters of correct magnitudes. The linear subsurface dynamic pressure is: pa = acos (kz — wt) cosh k(h + Zp) (3.8) cosh kh where A is the amplitude of the linear potential function. The wave number, k, and a= Apw frequency, w, are related through the linear dispersion relation: (w — kU)? = gk tanh kh (3.9) Frequency Nielsen (1986, 1989) established a method for determining the param- eters of a local linear approximation to waves from a pressure record. His method determined a local frequency and wave number that could be used to find the location of the water surface. A similar method is used here to determine the first estimate for the local frequency and wave number. Frequency of a sinusoidal signal of the form pa = acos(kx — wt) is available from the second derivative: 2 he oe (3.10) This approach requires an estimate for the value of the second time derivative of the record throughout each window. An estimate for the second derivative is directly available from the cubic spline of the measured points. This estimate, however, is very sensitive to random error in the measurements. Nielsen suggests estimating the value of the derivative through the use of divided differences. This solution works well for a smooth, sinusoidal record, but is also very sensitive to noise in the record, particularly in areas of small curvature; estimating a second derivative from a small segment of a noisy record can result in large errors. In order to accommodate the inevitable noise in an actual record, a different approach is taken here. 32 In each window, a third order polynomial is fit to the record in the least squares sense. The second derivative throughout ihe window can then be computed from this polynomial. By using more than four points in each window, any noise in the record is smoothed out by the least squares fit. This approach proved to be consistently reliable for artificially generated noisy records, resulting in reasonable estimates for the value of the second derivative throughout the window. A set of frequencies is computed from the estimate of the second derivative at each of the considered nodes. The mean of these frequencies is used as a first estimate for the local frequency of the record. Amplitude and Phase Once the frequency is known, the amplitude and spatial phase (kz) of a particular record can be found by rearranging the equation as a linear least squares problem by separating the cosine and sine components: Pai = acos(kzr—wt;) = bcoswt;+csinwt; (Soll) where a = Vb?+c? and kx = arctan(c/b). This results in the following matrix equation in the unknown amplitudes, 6 and c. coswt, sinwt, 175 coswt, sinwt> Paz b : = (alZ)} Cc O coswt; sinwty IE ii ? The system is determined if pa(t) is defined at two points. If J > 2, the system is over determined and can be solved in the least squared sense by common algorithms in numerical linear algebra libraries. This method provides a first estimate for the parameters: a, kz, and w. Refining the Linear Estimates These estimates for the parameters can be quite poor, as they are all dependent on the initial estimate for the second time derivative 33 of the record. In order to improve them, the estimates are refined by optimizing for the best least squares fit in the window. I minimize O(X) = (Pai — acos (kr — wt;))* xe ey (3.13) Xe (Gtr) Where P,,; is the measured dynamic pressure at time t;. The result of this optimiza- tion is a linear estimate for the pressure record that fits the measured record most closely in the window. Wave Number and Fourier Amplitudes Once the optimum initial frequency, phase and amplitude are found, the wave number is computed with the linear disper- sion relation (Eq. 3.9), and the first estimates for Fourier amplitudes are computed as follows: a cosh kh i ~ pw cosh k(h + 2,) oe Aj =a Ay (8.115) The value of a is not critical and a = 0.1 was found to be satisfactory. Water Surface The location of the water surface at the N nodes throughout the window is estimated from the linear pressure response function with stretching (Nielsen 1989): Pa(tn) cosh (k (A + pa(tn)/pg)) (3.16) Heats pg cosh k(h + zp) where z, is the vertical location of the pressure gauge, and 7, is the elevation of the water surface node at t,. pa(t) is computed from Eq. 3.8. 34 Nonlinear Optimization Once there is a reasonable first estimate for all the unknown values in the potential function, standard nonlinear optimization routines are adequate for this system. For the results given, the Levenberg-Marquardt algorithm was used as implemented by the MATLAB Optimization Toolbox (Grace 1992). If the optimization routine successfully finds a minimum, the solution is checked to see if a clearly spurious solution is found. Spurious solutions are identified by the following criteria: e Very large or highly variable errors. e First order amplitude smaller than higher order amplitudes. e Unrealistically large or small frequency or wave number. e Large discontinuity between windows in the predicted water surface. It is unusual for the optimization routine to converge to a spurious solution. It is far more common for the routine not to converge at all. If no solution is found, or a spurious solution is found, it is necessary to revise the parameters of the optimization. For the next attempt, the window width is increased by a factor of 1.5 (1.579), and the procedure is repeated. If this is not successful, the window width is increased once more to twice the primary width (27). When increasing the window width is not successful, the order of the potential function is decreased until a solution is found. If none of these adjustments result in an acceptable solution, the window is skipped, and future analysis must be interpolated through that point. These adjustments are most likely to be needed in the long, flat trough of a shallow water wave, where the window needs to be expanded to include some curvature to indicate the frequency. Difficulties may occur in addition near zero crossings, where there is also little curvature in the record. Here the effects of amplitude and frequency may not be independent, as: lim asin (wt) = awt (3.17) t—0 35 At the limit near the zero crossing, a and w have the same effect, and the optimization routine can not distinguish them. Widening the window to include more of the surrounding record avoids this difficulty, and is generally successful in this situation as well as in long, flat troughs. Another complication can be a record that is exactly symmetric about the crest of a wave. In this case, the equations on either side of the crest are not independent. This situation is unlikely to arise in a field record, and can easily be accommodated by using an asymmetric distribution of points in that window. 3.4 Theoretical Records To avoid complications from measurement error in the initial testing of the method, pressure records generated by Fourier Steady wave theory (Sobey 1989) were used. This approach also has the advantage of providing a solution with the complete kine- matics, to compare with results from the LFI method. Field or laboratory data that includes a full set of measured kinematics are not available. Fourier theory provides a near-exact solution for irrotational steady waves and can be applied at any depth (Rienecker and Fenton 1981; Sobey 1989). The use of analytically generated records ‘also allows the method to be tested under a large range of conditions, by varying water depth, wave period, and wave height. This is useful for establishing criteria for choosing the solution parameters, such as window width, number of nodes in each window, and order of the solution. Shallow Water To demonstrate the optimization procedure in a window, figure 3.2 summarizes the results from an initial estimate, before the final optimization. This is a window near the crest of a steep, shallow water wave generated by 18th order Fourier theory. The parameters of the wave are: 5m water depth, 3m wave height, 10s period, and zero Eulerian current, with the pressure record measured at the bottom. The parameters of the LFI solution are: sixth order (J = 6), and window width two seconds (79 = 0.27.) , centered 0.5s before the crest. The top plot shows the measured dynamic pressure as given by the Fourier theory — Actual pg © Estimated pa |— Actual © Estimated 7 ‘| — Bemoulli error DFSBC error — — MKFSBC error 0.2 0.4 0.6 — Actual pg o Estimated pj — Actual 7 o Estimated 7 :| —— Bernoulli error --+ DFSBC error “-:1—— MKFSBC error 0.2 0.4 0.6 0.8 Figure 3.3: Final results in a window 0.8 36 generated record, and the values computed from the Bernoulli equation (Eq. 3.4), with the parameters of the potential function generated by the initial estimate. The second plot is the water surface as predicted by Fourier theory and the elevations of the water surface nodes generated by the initial estimate. Note that the location of the actual surface is given in the plot, but it is not available to help determine the solution. These points were all generated by the method outlined in the previous section, with only the pressure record as a guide. The third plot shows the non-dimensional errors in the objective functions: the Bernoulli equation and the free surface boundary conditions (Eqs. 3.2, 3.3, and 3.4). These are the errors that are minimized by the optimization to find the solution. If the solution were perfect, the error in all equations throughout the window would be Zero. The initial estimates for the dynamic pressure and the water surface have order of magnitude and trend agreement. The errors in the objective functions are on order of .03 and show a systematic pattern, particularly in the Bernoulli equation. This pattern, and the fact that the sharp crest of the wave has not been matched indicate that a better solution can be found. The results after the nonlinear optimization are given in Fig. 3.3. At this point the prediction for the dynamic pressure is essentially exact. This is virtually always the case, as the pressure record is available, and the solution is optimized to that record. The predictions for the water surface are also extremely close. This is an impressive achievement, as location of the water surface was found only by minimizing the errors in the free surface boundary conditions. The non-dimensional errors in the Bernoulli equation and free surface boundary conditions are on order of .001, and show no clear trend. The lack of trend indicates that the remaining error is random, and a good solution has been found. The LFI method was able to capture accurately the crest of a steep shallow water wave. Figure 3.4 shows the results of the method for the complete shallow water steady wave. The LFI method finds the water surface and the kinematics on the surface essentially exactly. While these results show the complete wave, each of the indicated points is in the center of a separate window, and each window was computed com- 38 Dynamic Pressure Fourier fe) LFl Estimate on no on no eS a S| Figure 3.4: LFI predictions (J = 6) and analytical kinematics at the predicted water surface for shallow water wave x Awe /g? © Aqw?/g? + Asw?/g? Figure 3.5: Evolution of the solution for a shallow water wave pletely independently of the other windows. In this case, the parameters of the LFI solution are: primary window width of 2s (7 = 0.2T.), with a sixth order potential ‘function (J = 6), seven samples on the pressure record, and seven water surface nodes (I = N =7), resulting in 21 equations in 16 unknowns in each window. The window width of two seconds is one fifth of the period of the wave, and is a reasonable length of time to extend the locally steady approximation. Evolution of the Solution Fig. 3.5 shows the evolution of the solution as the window is moved along the wave. The top figure shows the non-dimensional values of the local fundamental wave number, k, and the local fundamental frequency, w. The importance of the local nature of the solution is apparent in this figure, as the solution varies substantially from window to window. The wave number and frequency both increase to maximum magnitude near the crest (4 < tw, < 7), and have lower values in the trough region (near tw, = 3 and tw, = 9). This pattern suggests that the kinematics in the crest region are similar to those of a much higher frequency lower 40 order wave, and the trough kinematics similar to a lower frequency wave. Deep Water Figure 3.6 shows the results of the LFI method for an entire deep water wave. The wave was generated by 10th order Fourier wave theory, with param- eters: 100m water depth, 10m wave height, 10s period, and zero Eulerian current, with the pressure record measured 10m below the mean water level. The parameters of the LFI solution are: Primary window width of 1s (7 = 0.1T-), with a fourth order potential function (J = 4), five samples on the pressure record, and five water surface nodes(J = N = 5), resulting in 15 equations in 12 unknowns in each window. Once again, the LFI solution matches the actual solution essentially exactly. 3.4.1 Choosing Parameters of Solution Unlike steady wave theory where only the order requires prior specification, this LFI application requires prior specification of four parameters: e Order of the solution (J) e Window width (79) e Number of nodes on the water surface (V) e Number of nodes on the pressure record (J) Experimentation has demonstrated that the resulting solutions are not highly sensi- tive to the values of these parameters; however, a reasonable solution will not result if these parameters are not selected judiciously. The experience acquired while devel- oping the LFI method on analytically generated records has provided some guidelines to help select appropriate values. Despite these guidelines, individual judgment must be used when applying the LFI method to a particular measured record. Number of Nodes on the Water Surface and Pressure Record Mathemat- ically, the system of equations is specified if there are at least as many equations as unknown parameters. In this case, /+ N > 4+ J. While meeting this criterion is 4] Dynamic Pressure Fourier LFl Estimate 12 Figure 3.6: LFI predictions (J = 4) and analytical kinematics at the predicted water surface for deep water wave 42 the minimum mathematical requirement, there are other factors that must be taken into account to assure a reasonable physical solution. The local water surface is rep- resented by N nodes in each window. Defining the surface with N = J+ 1 points achieves a representation for the water surface of the same order as the potential function being used. Care must be taken with very low order solutions; for example, 2 points would define the water surface to first order, but would provide no indication of the curvature of the surface. As a rule of thumb, at least three points should be used, regardless of order. Defining the water surface at J +1 points provides 2(J + 1) equations (f? and f*, at the J+1 points). To keep the order of approximation consistent, the Bernoulli equation is also applied at J+1 points within each window. Using J+1 water surface and pressure record nodes overspecifies the solution at all orders. This arrangement of equations proved effective for all the examples given in this chapter on analytically derived records. While that approach was effective on these few examples, it is important that the final solution matches the measured record closely. It’s possible that this suggested arrangement of points would allow the optimization routine to be biased towards the more numerous free surface boundary conditions, giving a solution that does not match the pressure record well. The Bernoulli equation could be applied at more points on the pressure record than the number of nodes defining the water surface. Increasing the number of points at which the Bernoulli equation is applied will shift the emphasis of the optimization to the measured record. Additional points on the measured record can also be useful for accommodating measurement noise that may be present in field or laboratory records. Order of Solution Similarly to high order steady wave theory, the order chosen for the solution is influenced by a number of factors including the height of the waves, the depth of the water, and the accuracy desired. As with steady wave theory, higher order results in greater accuracy at the expense of computational simplicity, and is necessary for larger waves and for shallow water. The following examples will help to provide guidelines for the order chosen. 43 Paws pg? Fourier Q ILI nw g B — Bernoulli 6 5x ‘—- DFSBC - - KFSBC -0.8 -0.6 -0.4 -0.2 (0) 0.2 0.4 0.6 tw Figure 3.7: Crest of a shallow water wave at order 1 Fourier Oo LFI — Bernoulli 0.6 ‘=. DFSBC - - kKFSBC -0.8 -0.6 -0.4 -0.2 (0) 0.2 0.4 0.6 tw, Figure 3.8: Crest of a shallow water wave at order 2 0.06 Paw? qze 0104 Fourier Oo. LFI nw g — Bernoulli ‘=. DFSBC - - KFSBC -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 tw, Figure 3.9: Crest of a shallow water wave at order 3 — Bernoulli ‘-- DFSBC - - KFSBC _. _. _ -0.8 -0.6 -0.4 -0.2 (0) 0.2 0.4 tw, Figure 3.10: Crest of a shallow water wave at order 4 44 45 0.06 Paw? i Fourier es : o LFI ato nw g a — Bernoulli 0.6 5 Xx -—: DFSBC - - kFSBC OF? Eagle Sa Fo SBIR STS Hy 8 -5 -0.8 —0.6 -0.4 -0.2 0 0.2 0.4 0.6 tw, Figure 3.11: Crest of a shallow water wave at order 5 0.06 Paw? pg? 0.04 Fourier Oo. 6LFI Nw g — Bernoulli ‘=. DFSBC - - KFSBC tw, Figure 3.12: Crest of a shallow water wave at order 6 46 Figures 3.7 through 3.12 show the crest of the same shallow water wave shown in Fig. 3.4, computed to orders J = 1 through 6. At first order, the non-dimensional errors in the Boundary conditions and the Bernoulli equation are of order 10~°. These are small errors, but the sharp curvature at the crest has not been captured. As the order is increased, the magnitude of the errors increases. The increase in the magnitude of the errors is due to the increase in the number of equations included. There are three additional equations included in the optimization for each increase in order (the two FSBCs at each additional water surface node, and one additional pressure record node), but there is only one more parameter in the solution vector. In order to best satisfy this system, the optimization finds a solution that results in slightly more error in each equation. The magnitude of the errors does increase, but the solution slowly converges to very precisely match the sharp crest at sixth order. An asymmetric distribution of points was used in this window to accommodate the symmetry about the crest. Figures 3.13 through 3.16 show the crest of the same deep water wave as Fig. 3.6, computed to orders 1 through 4. Even at first order, Fig. 3.13, the solution is very good. Deep water waves generally do not require a very high order solution, linear wave theory often being reasonably adequate in these conditions. It’s important to keep in mind that, although this solution is first order, it is still nonlinear, having found a minimum in the errors of the full nonlinear governing equations, and the frequency and wave number are free to vary, not being bound by the linear dispersion relationship. The local nature of the solution would allow it to change with time, accommodating an irregular profile at low order better than an equally low order global solution. In this case, first order provides an acceptable solution; however, the water surface is more accurately matched as the order is increased, and the solution converged well at higher order, so there is little penalty in using up to fourth order. For an irregular record, higher order is more likely to be successful in matching the irregularity in the record. As the examples show, deep water waves can be well represented at low order, while higher order in necessary in shallow water. —— Fourier -0.4 ; : : Oo LF 0.3 — Bernoulli ‘=. DFSBC - - KFSBC Fourier °o. 8 6LFI — Bernoulli ‘= DFSBC - - KFSBC Figure 3.14: Crest of a deep water wave at order 2 47 48 Paw? pg? Fourier (e} Nw. 9 — Bernoulli --. DFSBC - - KFSBC Figure 3.15: Crest of a deep water wave at order 3 Pyw? pg? Fourier I nw? 9 — Bernoulli --: DFSBC - - KFSBC Figure 3.16: Crest of a deep water wave at order 4 49 Window Width Choice of window width is a balance between including sufficient curvature in the record to identify the local wave response while minimizing the extent of the locally steady approximation. Another factor to be considered is the sampling rate of the data. The LFI method employs a cubic spline for interpolation among the measured points. This approach allows the selection window widths and sampling locations without being restricted by the data sampling locations. While this freedom is very useful, it can be deceiving, as the data are not truly continuous. As the data are not continuous, it is important to be sure that the window width chosen includes sufficient measured data points to justify the order being used. Unfortunately, field data are often sampled at a frequency of as low as 1Hz. In this case, a window width of one tenth the mean zero crossing might be on order of 1 second. This would provide a maximum of two actual data points in each window. It may not be appropriate to attempt a high order solution with such little hard data. Ideally, the LFI method should be used with data sampled at a higher rate. If that isn’t possible, wider windows, and perhaps lower order solutions, should be used. In the case of the theoretical records used in the previous section, the data sam- pling rate could be arbitrarily high. Without being restricted by measurement fre- quency, some guidelines for window width have been determined. A primary window ' width of one tenth the zero crossing period provides a good starting point. Using a smaller window often did not allow convergence of the solution. It was necessary to occasionally increase window width at some locations on the record, as discussed in section 3.3.2. When a primary window width is determined that works for most of the record, while needing to be increased at a only a few locations, the optimum width has been found. In the case of the shallow water wave discussed in the previous section, a window width of 0.17. was successful for up to a fourth order solution. Unfortunately, the fourth order solution did not fully capture the sharp crest of the wave. When a higher order solution was attempted, the optimization would not converge with the given window width. A wider window had to be used in order to obtain a solution of high enough order to capture the sharp crest of the wave, and thus a window width of 0.27, was used in that example. In general, it has been necessary to use wider windows for 50 higher order solutions. 3.5 Laboratory Measurements The following results use laboratory data collected by Murray Townsend and John Fenton at the Australian Maritime Engineering Cooperative Research Center at Monash University, Melbourne, Australia in June, 1996. The dynamic pressure, fluid velocities and water surface were measured at the same horizontal location along the flume; the pressure at a variety of elevations, and the horizontal and vertical velocities at an elevation of -0.8 meters, with the origin at the mean water level. The still water level was 1.55 meters, with the data sampled at 60 Hz. The wave flume is 52m long by 2.2m wide with two working sections 4m and 2.5m long connected by a ramp. A false floor had been built into the shallower working section to bring the depth to 1.55m. At one end of the flume is a hydraulically oper- ated bottom hinged wave paddle and at the other a beach that absorbs a minimum of 96% of wave energy across the range of operating frequencies. The measurements were taken approximately 22m from the beach end of the flume. The beach length is 6m. The wave surface was measured with a capacitance wave probe with typical cali- brations providing a correlation coefficient R? greater than 0.999. The pressure was measured with Druck PDCR35/D transducers located in the flume near the center of a long (2.5m long by 2m high) plywood board aligned with the direction of wave propagation. The measurement face of the probe was flush with the board surface to eliminate local flow separation. The R? values from calibrations were in all cases above 0.999. Fluid velocities were measured with a factory calibrated Alec electronics ACM250-A electromagnetic current meter which is capable of velocity measurements in two dimensions. One of the difficulties encountered with pressure measurements is the question of what has actually been measured. If the transducer itself is absolutely accurate, the pressure recorded will be the actual pressure at the location of the transducer. Unfortunately, the presence of the instrument is likely to alter the flow in its vicinity, 51 creating substantial differences between the pressure at the transducer, and the pres- sure that would exist if the transducer were not there to disturb the flow. Another source of error is the frequency response of the transducer (Raichlen et al. 1990; Ippen and Raichlen 1957). The frequency response of pressure transducers is not flat, and as a result, the recorded signal could be quite different from the actual pressure. In the experimental setup described above, the pressure transducers were mounted on a large plywood panel oriented in line with the flume. The measurement face of the instrument was flush with the surface of the plywood to minimize dynamic effects near the gauge. The panel was large enough for the boundary layer to be fully developed near the surface of the plywood, to prevent edge effects from the edge of the plywood affecting the measurements. This arrangement is expected to have resulted in minimal dynamic effects on the measured pressure. There is no information available about the frequency response of the pressure transducers used for this experiment. However, a spectral analysis of the records can help identify any potential problems. The top plot of figure 3.17 gives the measured water surface and subsurface pressure for a short segment of a record. There are some clear higher frequency fluctuations in the water surface that do not appear in the pressure record. In this particular segment, there is a sharp secondary crest in the first trough (near the 2s point) in the water surface record. This loss of high frequency information in the pressure records is confirmed by an examination of the discrete Fourier transform of the records. The second plot in figure 3.17 is the smoothed variance spectra of the two records from which the above segments were taken. The water surface record has a great deal more variance at the higher frequencies. Note that the spectrum of the pressure record has decayed to almost zero at w = 4s~! point, while there is still considerable variance in the water surface at that frequency. Linear wave theory predicts that the motion of high frequency waves decays with depth much faster than low frequencies. In this example, the non-dimensional water depth at the peak frequency is w*h/g ~ 1. This is generally considered intermediate depth water, and the decay with depth is expected to be moderate. At the frequency where there is essentially no energy in the pressure record, the non-dimensional water Water Surface and Dynamic Pressure 0 1 2 3 4 5 6 7 8 fs) 10 Variance Spectrum of Dynamic Pressure — Measured — — Predicted 0 5 6 7 10 Variance Spectrum of Water Surface — Measured — — Predicted Figure 3.17: Water surface, pressure record, and variance spectra of records from flume experiment 53 depth is w?h/g =~ 2.5, which is considered to be deep water. In deep water, the action of the waves decays very quickly with depth, so that little energy is expected to remain at half the water depth, the depth at which the pressure measurements were taken. The third plot in figure 3.17 is the spectrum of the subsurface pressure record. The solid line was computed from the measured record, and the dashed line was computed from the water surface record, using the linear pressure response factor: cosh k(h + z,) ap(w) = a,(w) cosh( Eh) (3.18) where a,(w) and a,(w) are the magnitude of the pressure and water surface amplitude spectra at a given frequency, and the wave number, k, is computed from the linear dispersion relationship: (w — kU)? = gk tanh(kh) (3.19) The predicted spectrum matches the measured spectrum very closely. This indicates that the loss of high frequency information in the subsurface pressure record is the result of the predicted decay with depth of the energy at the higher frequencies, rather than a result of the frequency response of the pressure gauge used. The final plot in figure 3.17 is the spectrum of the water surface record. The solid line was computed from the measured water surface, and the dashed line was computed from the pressure record, also using the linear pressure response factor. The predicted spectrum of the water surface matches the measured spectrum well near the peak frequency, but strongly over-predicts the high frequencies. If this method were used to predict the water surface, a somewhat arbitrary high frequency cut off would have to be established. An examination of the top plot in figure 3.17 reveals what appears to be two dominant free modes in the water surface. There is a large low frequency mode with a period of approximately 2s, and a higher frequency mode with a period of approximately 1s. If the LFI method were applied to the water surface record (Sobey 1992), the moving widow could capture a different dominant mode in each window. Unfortunately, the higher frequency mode has decayed with depth too much to be 54 Water Surface 0.25 - : : Original Record =e | : — LFl Estimate 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 time Figure 3.18: Measured and predicted water surface near the crest of a sharp wave in the flume clear in the pressure record; applying the LFI method to the pressure record results in predictions that capture the dominant low frequency mode, and miss the higher frequency mode. If the higher frequency information is not in the local segment of the record, it is not possible to reconstruct it. In the previous section, the LFI method was applied to pressure records generated by Fourier steady wave theory. These records had only a single free mode, and thus the LFI method was able to identify that mode. In addition, in shallow water, there is little decay with depth of the wave action, and the LFI method could be expected to be effective, as it was on the Fourier record. In deep water, the method can be effective if the pressure gauges are located near enough to the surface. For the deep water steady wave presented previously (figure 3.6) the pressure record was taken from a depth of one wave height below the water surface, one tenth of the total depth. At this depth, there is adequate information to identify the action at a wide range of frequencies. Given the limitations of the data, a segment of the record was chosen in which there did not appear to be a substantial high frequency component to the water Dynamic Pressure — Measured Record Waiter Surface © LFI Estimate Figure 3.19: Segment of record from flume measurements 59 56 surface record, and thus the effect of the high frequency decay should cause few problems, and the LFI method can be applied. While this particular segment was chosen to minimize the effects of the high frequency decay, it was not possible to eliminate it entirely. Figure 3.18 shows the measured water surface and the LFI predictions of the water surface in three windows in the vicinity of a sharp crest in the record. The solution parameters are: J = 3, ™ = 0.3T., and J = N = 4. Both the pressure and the velocities were measured at an elevation of 0.8m below the MWL, with the mean water depth 1.55m. The LFI method located the water surface fairly accurately, although it clearly missed the very sharp peak of the crest. Varying the solution parameters did not improve the accuracy at the crest. In fact, a smaller window width resulted in a discontinuity in the water surface between windows. The inability of the LFI method to capture the sharp crest is likely due to the missing high frequency information. If the high frequency part of the signal is missing from the examined pressure record, there is no way to recapture that information. In order for the method to be more effective in this situation, the pressure would have to be measured at a depth closer to the water surface. Figure 3.19 shows the results of the LFI method applied to the entire segment of the record from which the crest was taken. The water surface is captured fairly well, as are both the vertical and horizontal velocities. 3.6 Discussion The given results demonstrate the potential of the LFI method in the interpre- tation of submerged pressure traces in a variety of conditions. In the case of steady waves, the LFI method accurately computed the detail of the wave, using only data from a small window in time. In particular, the method was able to capture the pronounced sharp crest of a steep, shallow water wave. The method did not perform as well on laboratory records, failing to capture some of high frequency detail in the water surface. This is due to the limitations of subsurface pressure records, where much of the high frequency information is missing of from the record. The method is likely to perform better in shallow water, or with records that are measured closer to the water surface. Despite this limitation, the method was able to capture much of the detail of a irregular laboratory record, and provide fairly accurate estimations of the local kinematics. The analysis of regular waves provides guidelines for the parameters to be used in the analysis of irregular waves. In shallow water, higher order solutions and wider windows must be used than in deep water. Window widths of one fifth of the zero crossing period and a sixth order potential function are adequate for the shallowest waves, and window widths as small as one tenth of the zero crossing period and a third order potential function are adequate for deep water. The laboratory results indicate that the method requires good precision and care in the measurements. Any high order method demands very accurate data, but this need is exacerbated by the ill posed problem of determining the near surface kinematics from a subsurface record. While not particularly sensitive to random noise in the record, the decay with depth of the high frequency information makes it very difficult to capture the high frequency modes. Fundamentally, the LFI method is designed to capture the dominant free mode in each given window. As the higher frequency modes decay faster with depth, and if the measurements are taken far below the surface, the dominant mode will always be one of lower frequency modes. This difficulty would be exacerbated by any limitations in the frequency response of the gauges. The computer resources required for the method are substantial, but not pro- hibitive. As computers continue to get faster, computation time is not the considera- tion that it once was. The method was developed and all computation done using the MATLAB computational package. MATLAB is an interpreted language that provides a very easy to use interface to a robust and complete library of computational and visualization routines, allowing for rapid development of new methods. Being an in- terpreted language, the resulting code is not as fast as it might be if the routines were programmed in a compiled language, such as Fortran or C. The computational speed is also affected by the degree to which the program pauses to provide visual output. Perhaps a better measure of the computational intensity of the method is a count of the total number floating point operations (flops) used to compute the solution. 58 The examples in this dissertation were all computed on a Intel 90MHz Pentium processor based PC with 32MB of RAM, running the Linux operating system. As an order of magnitude estimate, the shortest computation time for an example in this chapter is 4.8 minutes and 7.8 x 10° flops for the complete laboratory record presented in figure 3.19. The longest computation time required was 71 minutes and 582 x 10° flops for the shallow water wave presented in figure 3.4. There are a number of reasons for this large variation in computation time. The first is simply the number of window solutions computed. As each window solution is computed separately, the computation time increases linearly with the number of windows computed. If computational time is a concern, this can be taken into account when choosing the output spacing. The other reason that the shallow water wave takes much longer to compute is that there are a number of windows that did not converge with the first set of computational parameters. The optimization is run for many iterations to ensure that it won’t converge. As the parameters are varied, the computation is repeated. This process takes a great deal of time. With further development, it may be possible to determine a set of criteria for the computational parameters by examining the segment of the record in a given window. This would be far more efficient than simply attempting a solution with a variety of values until convergence is achieved. 59 Chapter 4 Three Dimensional LFI Theory The previous chapters were concerned with the determination of wave kinematics from the analysis of a time series of a single physical quantity at a single location. By assuming that the flow is two dimensional, a reasonable approximation of the wave kinematics can be found. Unfortunately, this analysis does not give any information about the directionality of the sea state. There are some processes in which the wave directions are directly important, such as sediment transport. Even in situations where the wave directions are not directly important, it has been shown that omitting _ the directional nature of the sea results in substantial errors in the prediction of scalar quantities, such as the maximum velocities and accelerations in a measured wave crest, or the resultant forcing on structures (Dean 1977; Forristall et al. 1978). In order to capture the directional nature of the sea, an array of instruments must be used. The result is a set of time series of a single physical quantity at a number of different locations, or a set of different physical quantities at the same or different locations. This chapter outlines a method for determining the directional kinematics of irregular seas that can be adapted to accommodate virtually any combination of such time series. 60 ~ - a surface (7) ae y PUUn Hees s\ ose eeeosss Tea MUG (@ = 0) sea bed Figure 4.1: Coordinate system used for 3-d method 4.1 Three Dimensional Seas The development of the two dimensional LF] method was closely anchored to the very complete understanding of two dimensional steady waves. In contrast, the understanding of three dimensional wave fields is not nearly as complete. Much of the literature on three dimensional seas attempts to describe the motion through the use of a directional energy spectrum. Far less attention has been directed to the determination of the detailed kinematics of directional seas. 4.2 Problem Formulation The governing equations for three dimensional gravity waves are a straightforward extension of those in two dimensions to include the third dimension. The flow is taken to be irrotational and incompressible, and thus the kinematics can be represented by a potential function, d(x, y, z,t), in a Cartesian coordinate system (Fig. 4.1), where: _ Of _ Od (aka) u and v are the horizontal velocities in the z and y directions, and w is the vertical velocity. 61 Field Equation The field equation is mass conservation for irrotational flow, rep- resented by the Laplace equation: Od Oo O76 Vv? = a Me Ox? i Oy? x Oz? 0 (4.2) Boundary Conditions The boundary conditions are the bottom boundary condi- tion (BBC) for a horizontal bed, as => = W 0 at z=—h (4.3) the dynamic free surface boundary condition (DFSBC), 1 sree re to )tqp—e=0 a 2eR (4.4) SS oP US POR KS OS at =7 (4.5) u v Zz : . where 77 is the elevation of the water surface. As with the two dimensional formulation, the kinematic free surface boundary condition requires the gradients of the water surface (22, on and a): Knowledge - of these gradients is likely to be limited. In the case of an array of water surface measurements, estimates for the gradients could be computed by interpolating among the measured elevations, but this would result in a low order estimate, and compound the inevitable error in the measurements. In the case of subsurface measurements, there are no data as to the location or the gradients of the water surface. To eliminate these gradients, a modified kinematic free surface boundary condition is defined by differentiating the DFSBC following the motion (Longuet-Higgins 1962), as was done in two dimensions in Section 2.1. MKFSBC = —g(KFSBC) + =. (DFSBC) = Od Ov Ow OP tow +20 +05 + Oa v Ow (4.6) Ou sp CO ap CO a0 (I) at 2 =I This form does not include the gradients of the water surface, and all the terms are defined by the potential function. Applying both the modified kinematic free surface boundary condition and the dynamic free surface boundary condition completes the formulation without the need for knowledge of the gradients of the water surface. Observational Equations The field equation and the boundary conditions de- scribe a free surface potential flow. In steady two dimensional wave theory, a wave height and periodicity in space and time are specified to complete the formulation. In the two dimensional LFI theory (Chapter 2), a form for the potential function is chosen, with the parameters determined to fit the free surface boundary conditions and the measured record in a small window in time. In three dimensions, the pro- cess is essentially the same. A three dimensional form for the potential function will be presented, with parameters that are found to fit the measured records and the boundary conditions. In order to define a solution that fits the measured record, observational equations are established. These equations are defined to make use of the particular quantities that have been measured. In the case of an array of water surface measurements, they are the free surface boundary conditions, applied at the measured elevations and hor- izontal locations at a number of points in time throughout the window (Chapter 5). In the case of an array of pressure measurements, they are the Bernoulli equation, applied at the elevation and horizontal locations of the measurements, and also at 63 a number of points in time within the window (Chapter 6). The method could be adapted to virtually any combination of measured physical quantities by establishing appropriate observational equations. In addition to the aforementioned water surface and pressure measurements, these could include water surface gradient or accelera- tions measured by wave buoys, subsurface velocity measurements, or any combination of these. 4.3 Formulation of the Solution in each Window 4.3.1 Background Two dimensional steady wave theory provided the inspiration for the development of the two dimensional LFI method (Chapter 2). Unfortunately, the literature does not provide as solid a basis for nonlinear interpretations of directional seas. There has, however, been some work that can be used as a basis for a directional LFI method. In an early attempt to explore the nonlinear nature of directional seas, Longuet-Higgins (1962) computed the interaction of two intersecting steady waves in deep water through the use of a double perturbation expansion in the steepness _of the waves up to third order. The result was a potential function that contained terms representing the higher order interaction between the phases of the intersecting waves: bo” = fi(z, 51) + fo(z, So) po) = fs(z, 51 + 52) + faz, 51 — So) $°) = f5(z, 251 + S2) + fe(z, 2:51 — S2) + fr(z, Si + 252) + f(z, 51 — 252) Sn = (Kn -X — Wnt + An) (4.7) where n = [1,2], d'”) is the mth order potential function, S; and S> are the phase functions of the two waves, k,, is the vector wave number of wave n, x is the horizontal position vector, w, and a, are the angular frequency and initial phase of the nth wave. The functions, f;, were determined algebraically by expanding the free surface boundary conditions in a Taylor series about the mean water level, and solving for 64 the potential function at each order. Nonlinear frequency modulation was not taken into account, so the frequencies and wave numbers of each wave independently satisfy the linear dispersion relation for water of infinite depth: w, = glk,| (4.8) fi and fz are independent of each other, and are the familiar linear wave solution: ja S erie sin $1 i mre sin So (4.9) The higher order terms are all functions of both waves, and thus include the interac- tion of the two waves. A number of investigators subsequently expanded upon this work. Hsu and Chen (1992) presented a detailed examination of Longuet-Higgins (1962), pointing out diffi- culties that arise from the assumption of the linear dispersion relation. They presented a more mature analysis, including higher order modulation of the wave frequencies, and higher order self interactions of the individual waves. This resulted in a com- plete theory up to third order for two intersecting waves in deep water. Hsu and Chen also proposed a systematic ordering in the phase relationships generated by the interactions of two steady waves to arbitrary order. Expanding upon this work, Ohyama, Jeng and Hsu (1995a) extended the per- turbation expansion method in a number of ways. The most recent version of the method accommodates any number of waves, allows for water of finite depth, and is accurate to fourth order. This last method can compute the water surface and full kinematics of a highly irregular sea, as produced by a large number of intersecting waves. A more detailed discussion of the method is given in Appendix A. Ohyama, Jeng and Hsu’s work suggests a form for the potential function that can be used for a local method for the recreation of three dimensional kinematics from arrays of wave measurement devices, including water surface arrays, pressure arrays, directional current meters, or any combination of these. The general form for the potential function representing N intersecting steady 65 waves at order J is: P(z,y,2,t) = Urx st Uyy ate J=i1| J=TNTt lim| J s cosh Ajj(h+2) . A=—Jie=—J+li| — iy=—J+ 3 N=1 lim LI = (GCip@ayce2 stay) Of] = (4151 + 2252 +--- + in Sy) Sn = (kent + ky ny — wnt + On) “= ae 2 if aan 2 where U, and U, are the components of the depth uniform Eulerian current, k,,, and ky, are the vector wave number components of the nth wave in the z and y directions. There is a separate summation for each wave considered. As the notation in Eq. 4.10 is quite confusing, examples with two and three waves are given below. With two waves (N = 2): P(2,Y,2,t) = Ure + Uyy + J J—|t1| cosh 1&2, Blo + z) : : - oe SS As) Tepe han aue Ts (215) + i252) (4.11) LEG a And with three waves (N = 3): P(z, Y, ais) = Uf ate Uyy AF J—|i1| J—|t1|—22| d INA, By eal Se 2) oh pe : ye ST ye Ae se i fasis(h + 2) ss sin (2154 + 1252 + 2353) Goaln KE A Bll y=-J in=—J4 iy] ig=—J +e | +2 CDE (4.12) 66 Sn = (Kent + ky ny — Wnt + On) fe = [1,23] (t1k2.i: + takzin + t3kz,ig)? + (21k yi: + 22h yin + takty,is )? This form for the potential function exactly satisfies mass conservation, in the form of the Laplace equation, and the bottom boundary condition for a locally horizontal bottom. It allows for high order representation of each of the steady waves, as well as the interaction of each wave with every other wave. The balance of the solution is specified by the full three dimensional free surface boundary conditions, Eqs. 4.4 and 4.6. 4.3.2 Dynamics The dynamics are available through the Bernoulli equation in three dimensions: ag ot 1 =e shag ts aa a =1() (4.13) The dynamic pressure is the difference between the total and hydrostatic pressure (Pa = P — Ph)- The Bernoulli Constant In Chapter 2, an explicit and exact expression for the Bernoulli constant, B, is given for steady two dimensional waves. (Eqs. 2.14 and 2.15). A similar expression can be established for unsteady three dimensional waves. The Bernoulli equation is applied at the bed: Od, 1 i >) a 7 CO) eos 2 (4.14) where the subscript, 6 indicates the value at the bed. The two dimensional approach is based on analysis first presented by Longuet- Higgins (1975). In that work, the analysis was applied to steady waves. With steady waves, the average over either time or space of the pressure at the bed is the hydro- static pressure. With unsteady waves, the pressure at the bed averaged over space at any given time, or averaged over time at any given location, is not necessarily the 67 hydrostatic pressure. If, however, an average is taken over both time and horizontal space, the result must be the hydrostatic pressure. b = pg(7] — 2) (4.15) 3 with the double overline indicating an average over time and horizontal space. When the Bernoulli equation at the bed is averaged over time and horizontal space, the average time gradient of the potential function is zero, resulting in a simple expression for the Bernoulli constant: (u2+ v7) +97 (4.16) For the potential function given by Eq. 4.10, and z = 0 at the mean water level (7), B becomes 1 ree rat Ay ky \? Psa (awmy2k So. esa 4. 2 UE + U;) a 4 x ~ (= oF) ey) y=—J in=—J+ NI} lim| giving a complete expression for B as a function of the parameters of the potential function. 4.4 A Local Two-Intersecting-Wave Theory While a large number of intersecting waves could capture a sea state of virtually any complexity, it would be difficult to distinguish between the effects of each in- dividual wave. It is important to remember that the familiar directional spectrum description of a sea state is an integral description. While there are many different frequency and direction modes represented in the spectrum, they do not necessarily all play a significant rule in the kinematics the entire time. In fact, when observing irregular seas, it is often the case that at any given time, there appears to be a single dominant wave of a particular frequency and direction. Over time, there is a series of such dominant waves, each of a different frequency and direction. The integral effect of this process is a broad directional spectrum, but if time is separated into individual 68 small windows, in each window there may be only one or two free waves dominating the motion. Taking advantage of the local nature of the approach, only one or two intersect- ing free waves are considered in each window. While unlikely to capture the full complexity of a broadly directional sea, this approach should be able to capture the dominant modes in each window. The direction, frequency, and amplitude of the dominant modes will vary substantially from window to window as they do in the two dimensional approach, having an integral effect that includes a large variation in directionality. The two wave method is expected to be particularly effective for the case of standing waves or short crested waves, as would be found near a reflecting surface when the incident wave field is almost unidirectional, as it often is in shallow water. Expanding Eq. 4.11 to fourth order, the potential function for two intersecting waves is: 20 P(x, Y, B..h) =U,2+ Uyy + Sy A; (KA, Kz) sino; (4.18) cl where: : : cosh K;(h + z) : 00;7 Oo,” [Mo L632) = —$—— ——— CV ie-we eran er M2) cosh hh B Ox i Oy Sy = (kear + kyiy — wit + a) So = (krat + ky2y — wot + az) At first order: Oo1= (5) On) = (S2) At second order: 03 = (25;) 04 = (252) o5 = (Si + S2) a6 = (51 — 52) 69 At third order: o7 = (351) og = (352) Og = (25; + 52) G19 = (25; — So) 011 = (Sy + 252) 012 = (5) = 2S2) 13 = (45)) O14 = (452) O15 = (35; + $9) Cig = (35; — $2) O17 = (25 + 252) o1g = (25, — 22) O19 = (5; + 352) 729 = (51 — 3S) 01,02,03,04,07,08, 013, and O14 are self-wave interactions. The balance of the o; are Wwave-wave interactions. Because the phases are arguments of the sine in the potential function, a sim- ple expansion of Eq. 4.11 results in redundant terms that have been removed. For example, if there are two components: By sin(S;—S2) and Bposin(S2—5;) (4.19) Because sine is an odd function: By sin(S2 — $,) = —Bz sin($; — S2) (4.20) the two components can be combined into a single component: (B, — B2)sin(.S; — Sz) = A;sin( Si — S2) (4.21) If both components were included, the effects of B; and B2 would be indistinguishable in the optimization; they must be combined into the single coefficient, A;. This form for the potential function results in 28 unknowns at fourth order (20A; plus kz, ky, w, and a for each of the two intersecting waves) that must be found 70 for each window in time. At least 28 independent equations must be defined by applying a combination of each of the free surface boundary conditions (Eq. 4.4 and 4.6) and the observational equations at a number of nodes throughout each window. The specific combination of these equations will be determined by the layout of the instruments. 4.4.1 Kinematics In order to apply the free surface boundary conditions, the full kinematics must be known. While the kinematics are completely specified by the potential function, Eq. 4.18, some of the algebra may not be obvious. The full expressions are provided here. ; MT (GI Zh) = = =| ULt Sy Ac ((K;h, K;z) cos o; (4.22) jail ad ac v(z,y,z,t) = fe te A @(K;h, K;z) cos 0; (4.23) 7=1 ad Zl w(2, Yy, Zab) = BE = S- AK; SC(KG, Kz) cos 0; (4.24) i 2! ; I : oP (eus zt) = 2, Ase (Kh, Bjz) cos 9; (4.25) se Eb Ae ope on Z Asie \ pS een (= aa) (E28) 71 71 where: sees CORBI ae cosh h;h oa ey , Sul AG((lo ) ( S2) og = (352) 010 = (2S; ei 52) 012 = (5; a 252) O14 = (4582) O16 = (35; — $2) 01g = (25, — 25>) 020 = (5; — 352) It might be possible to include a larger number of intersecting free waves in order to capture a broadly directional sea, but it would introduce additional complication in distinguishing the effects of the individual waves, and will not be considered here. Both of these potential functions satisfy the field equation (Laplace, Eq. 4.2) and the bottom boundary condition (Eq. 4.3) exactly. The balance of the solution is determined by the free surface boundary conditions: the modified kinematic free (63) surface boundary condition (f*), ro Ou Ov Ow A — ) 2 —— ) if a 1 Ie (un tua, twa) oF eee oF ape! ar renee 6) Ox Ox (5.3) rey te EI ra 6) Oy Oy ie cps aa Sih 8 w— )— —__ = 3 = UD a oh iO we Ae a n and the dynamic free surface boundary condition (f”), 0d 1 ae J et a Gig il) at B= i (5.4) with the Bernoulli constant (B) defined as: 1 Lea Ake \? TH Zire D\ oS i P 2 eeu, ep (<= 7) oe) The problem of determining the kinematics of irregular waves from a set of mea- sured water surface traces is a mathematically better posed problem than interpreting a subsurface record. The flow is governed by the elliptic Laplace equation (Eq. 4.2), so that the solution is determined by the boundaries. While the complete boundaries of the solution domain are not known, the boundary conditions and location of the boundary are known at both the top and bottom of the solution domain. The bottom boundary condition is well defined, and the location of the water surface is measured at a few locations in space and many points in time, allowing the direct application of the free surface boundary conditions. This is in contrast to working with subsurface records, in which the location of the free surface must be determined in order to apply the free surface boundary conditions. The need for horizontal boundary conditions is eliminated by the assumed periodicity of the chosen potential function. However the fundamental wavelength(s) and period(s) must be found as part of the solution. When working with a point measurement, the fundamental frequency is fairly well defined by the time evolution in the window chosen, as long as the window is wide enough. There is no direct information available about the spatial evolution of the signal, however, so the wave number is determined only through the application of 76 the free surface boundary conditions. An array of measurements provides information about the spatial evolution of the signal, helping to better define the fundamental wave number. This usually results in faster and more robust convergence of the optimization. 5.2 Formulation of the Optimization The formulation of the optimization for the LFI method as applied to the in- terpretation of an array of water surface measurements has a great deal in common with the formulation for a subsurface pressure record (chapter 3) and a point water surface measurement (Sobey 1992). Less detail is presented here than in the previous chapter, but the framework will be presented, with an emphasis on the additional information necessary for applying the method to an array of measurements. Observational Equations The governing equations presented in the previous sec- tion represent a free surface potential flow, with one or two components propagating in an arbitrary direction. The observational equations are the equations in the system that force the solution to fit the given measured record. As the location of the water surface has been measured, these are the free surface boundary conditions (Eqs. 5.3 and 5.4), applied at the horizontal location of each of the nodes in the array. Sufficient independent equations are defined by applying the boundary conditions at a number of times along the measured records, within the window in time considered (Fig. 5.1). The solution is the set of parameters in the potential function that result in the least error in the FSBCs. In order to specify the solution, there must be at least as many independent equations as there are unknown parameters in the system of equations. The free K m,n surface boundary conditions ( and f?_) are applied at each of the N measured locations and at M time samples in the window, resulting in 2MN independent equations. In the single wave case, there are 4 + J unknown parameters sought in Kq. 5.1 (kz, ky, w, a, and A,...A,z) in each window, so that if 2MN > 4+ J, the solution is specified. In the two wave case, there are 2)> 7 +8 unknowns in Eq. 5.2 Ue t Figure 5.1: Schematic of system of equations in a window (2k, 2ky, 2w, 2a, and A,...Azy: 10 at first order, 14 at second order, 20 at third order, and 28 at fourth order). The solution is specified when 2MN > 2507 +8. This system results in the following least squares optimization. D 2 minimize O(X ass SE ( 1G a Ma Uncen SPS Oa Ontbnenta) (GS) m=1 n=1 where: X= (ke Sp Ds Qa, vAre D009 Aj) for the one wave case, and: X = (koi, ky1, 1, 01, ke,2, ky,2,W2, 02, Ar,--. , Az) for the two wave case. z, and y, are the coordinates of the nth gauge, and mn is the measured elevation of the water surface at time t,, and gauge n. As with the analysis of a pressure record, overspecification can be helpful in accom- modating the measurement errors in field records. More than the minimum number of time samples in the window may be required to define the shape of the water surface in each window. This is less likely to be necessary with two waves than with one, as the number of unknown parameters is much larger. It should be kept in mind, however, as it would be a factor for low order solutions with a large number of measurement locations. 78 5.3 Computation Methods The LFI method for an array of water surface measurements can be broken down into a similar sequence of steps as with the analysis of a pressure record: 1. Pre-processing of record. (a) Determine estimate for level of noise in the record. (b) Determine estimate for magnitude and direction of Eulerian current. (c) Define a set of continuous records from each gauge from the discrete ob- servations by cubic spline interpolation. (d) Specify spacing of output locations. (e) Compute mean zero crossing frequency. (f) Non-dimensionalize record and all parameters. 2. Primary values of numerical solution parameters are chosen. (a) Window width (70) (b) Order of solution (J) (c) Number of time samples of the water surface records within each window (M) 3. Global solution is computed on an entire wave to provide first estimates for local optimization 4. For each selected output location, a window in the record is defined, and the an LFI solution is computed. (a) Initial guess for the optimization is determined from the global solution. (b) Full nonlinear optimization for all unknown components of the potential function is computed. (c) Results are checked for spurious solution. 79 i. If no solution, or a spurious solution, is found the solution parameters are adjusted, and the optimization repeated. ii. If a good solution is found, progress to the next window. 5.3.1 Pre-Processing of Record Accommodating Measurement Error Measurement error can be a major source of difficulty with any high order data interpretation method. Local methods can be especially sensitive, as each window solution relies on the detail contained in a small segment of the record. In applying the LFI method to a single pressure trace, Sobey (1992) found it necessary to apply a simple moving average filter to field and laboratory data. In his work, the primitive kinematic free surface boundary condition was used, requiring an estimate for the local gradients in the water surface. In the current work, a modified version of the boundary condition (Eq. 5.3) is used which does not require these gradients. This makes the method less sensitive to noise, so it was not necessary to apply any filtering for the results presented here. If there is substantial noise in the measured record, the system of equations can be - substantially overspecified, allowing the least squares optimization to accommodate the noise in the record. When this is possible, it is preferable to applying a smoothing filter to the record, as it does not impose any assumptions on the nature of the record. However, if the error bands are very large on the data, it may still be necessary to apply filtering to the raw measurements. The Mean Water Level The mean water depth, h, must be specified as part of the potential function. As a time series of the water surface is provided, it is a simple matter to compute a mean of the measured records. The mean should be taken over a period much longer than a typical wave period, but short enough to accommodate changes in the local water level due to astronomical and storm tides. In keeping with the local nature 80 of the approach, this is the local mean water level, rather than a global still water level, which might be different, and would be less appropriate. As with subsurface measurements, the precision of this value is not critical, as h is only used to locate the origin of the coordinate system for the potential function. The Eulerian Current The Eulerian current, U, and U,, completes the description of the propagation medium of the waves, and an accurate estimate of its magnitude and direction is important. The measured water surface traces provide no information about the current field, so the information needs to be provided from other sources. See Section 3.3.1 for a detailed discussion of this important parameter. Spline of the Records As with the previous chapter, a set of continuous records is computed by cubic spline interpolation among the measured points at each gauge, allowing complete flexibility in the choice of window widths and location of samples in time of the records. Output Locations The spacing of the desired output locations must be chosen to determine the placement of the windows. Each window is computed independently so there is no restriction on the spacing of the output windows. In addition to selecting the output spacing in time, in can be useful when using an array of instruments to determine a single central location within the array at which to define the solution. Non-dimensionalization In order to prevent spurious solutions due to the comparisons of errors of different dimensional quantities, all parameters and variables are non-dimensionalized before computation with scales defined by the same physically identifiable parameters used 81 in Chapter 3: the mass density of water (g), acceleration of gravity (g), and the mean zero crossing frequency of the measured record (w-). Length scale = g/w? Time scale = 1/w, (5.7) 3 Mass scale = ee Ww. 5.3.2 Optimization Procedure The nonlinear optimization in the LFI method for the analysis of arrays of water surface measurements is very similar to that used for the analysis of a subsurface pressure trace. In the single wave approach, the solution is somewhat easier. The potential function used by the single wave approach is almost the same as that used with a point measurement (both Chapter 3 and Sobey (1992)), with the addition of a directional component to the wave number. This results in a single additional unknown, but the array of measurements provides at least three times as much data, with gauges at at least three spatial locations to specify uniquely the directional structure of the sea. The optimization tends to converge more rapidly and robustly than with a single point measurement. When applying the method with the two wave potential function, there are many more unknowns, and the optimization once again becomes somewhat tenuous. As with the previous chapter, in is essential to identify a good initial estimate to reduce the chances that the optimization will converge to a spurious minimum. Global Solution The strength of local methods is that they seek a single solution for only a short segment of a record at a time. The downside is that a single short segment often may not contain sufficient information to identify the directional nature of the wave field. In the method presented in the previous chapter the initial estimate could be computed from the data in the current window. In contrast, when working with spatial arrays, it is necessary to examine a larger segment of the record, for example 82 an entire wave from crest to crest or trough to trough, to get a general sense of the directional trend of the wave field. This global solution can then be refined to fit a small segment very closely. This approach is used to establish the initial estimate for the optimization procedure in each window. For the initial estimate to the global solution, it is assumed that the water surface records can be approximated with linear wave theory: n = acos(kztn + kyyn —wtm +a) (5.8) for the single wave method, and n = ay cos (een. + kyaYn cae Wit iy a1) (5.9) +2 COS (Kr22n =F ky 2Yn i Wotm oP a2) for the two wave method, where a, = A,w/g, and A, is the amplitude of the linear potential function of the nth wave. Directional Trend The first step is to determine the directional trend of the wave field. This determination is accomplished by examining the gradients of the water surface throughout the segment considered. Estimates for the spatial gradients and the elevation of the water surface at the center of the array are computed by finite difference approximations. The water surface is expanded in a two dimensional Taylor series about the center of the array: = One One UN(25 2) = ih a Bs) ae (@ = Ba) ae (y — Yc) By +... (5.10) where xz, and y, are the coordinates of the center of the array. This expansion is written for each of the N measured gauges, at each of the M points in time, resulting 83 in a system of linear equations that can be expressed as the following matrix equation: Ll B= 2) We— Mn iol 2 3 o o ina If fe=2y Ve = Op 2,1 Net Ne2 - - - Nem One. 9Ne.2 OneM | — Ox Ox AW RE Ox ONne,1 ONc,2 One,M oy Cun so ae ay LT Be ay We — hy WN 0 = « « ii (5.11) O 6 fo) ead where zx, and y, are the coordinates of the nth gauge in the array, 7m, se and Stent are the water surface and gradients of the water surface at the center of the array at time, t,,, and Nr, 1s the water surface elevation at gauge n and timet,,. If there are three gauges, the gradients are uniquely specified. If there are more, the system is solved in the least squared sense. The water surface is traveling either toward or away from the direction of the water surface slope depending on whether it is moving up or down at the time. The direction is thus determined by the spatial gradient of the water surface, and the sign of the time gradient, yielding a set complex direction vectors: «(Om \ ( ONm Om By = sign (Fe) (Fe + i) (5.12) A cubic spline of the water surface at the center of the array (n.) would provide an explicit piecewise polynomial expression that could be directly differentiated to obtain the gradients in time of the water surface. These estimates would be very sensitive to measurement errors in the record. To obtain more robust estimates for the time gradients, a smoothing cubic smoothing spline (de Boor 1978) is employed instead. This algorithm provides a smoothing parameter, p, that can be set at any value between 0 and 1, where p = 0 results in the linear regression fit, and p = 1 results in the “natural” cubic spline. The smoothing parameter may be varied to accommodate varying levels of noise in the record. For the examples in this chapter, p = 0.9 was found to be satisfactory. 84 Figure 5.2: Direction vectors for a short crested wave For the single wave method, it is expected that the individual directions will vary around a central dominant direction, and the propagation direction estimate is the orientation of the mean of the complex direction vectors: 6 = angle (3; SD 3. (5.13) For the two wave method, it is expected that the propagation directions at each point in time will vary around two dominant directions. In the case of a standing wave, these two directions are clearly defined. In a more complex sea, it is not so straightforward. To separate the two dominant propagation directions, the set of direction vectors, Der is divided into two sets. The two dominant directions are defined as the angles of the means of the sets of direction vectors within 7 greater 85 than and less than the mean direction (see Fig. 5.2): M, 1 = e2 6, = angle (a: ) 3. where: 0< angle( D m) <0+7 (5.14) 1 = = 62 = angle (si SS B. where: 6—-7< angle( Dm) © ¢ : wave ray water surface gages Figure 5.3: Layout of water surface array 5.4 Theoretical Records 5.4.1 Single Wave Records The following results are from “measured” data generated by Fourier steady wave theory (Sobey 1989). Fourier theory provides a near-exact solution for steady waves that provides the complete kinematics. This approach provides a complete set of data to compare with the results of the LFI method, without the complications introduced by the inevitable errors of data collected in the field or the laboratory. Theoretical records were used also because field or laboratory data that included a full set of measured kinematics to compare with the results are not available. A triangular array of water surface measurements was used, as indicated in Fig. 5.3. The array is an equilateral triangle with the same dimensions as the DWG-1 pressure array (Howell 1992). Three measurements were chosen, as three is the min- imum number required to provide directional information. Additional measurements would provide overspecification, and can easily be accommodated in the formulation. With actual field data, additional measurements are recommended, as increasing the number of instruments would provide redundancy in the case of instrument failure, as well as helping to accommodate measurement error. 91 Kinematics at the predicted water surface —— Fourier LFl Figure 5.4: LFI predictions (J = 3) and exact kinematics at the center of the array at the predicted water surface for a steady deep water wave 92 Fig. 5.4 is a steady deep water wave generated by 10th order Fourier theory with the following parameters: wave height = 10 m, period = 10 seconds, water depth = 100 m, zero Eulerian current, and direction of travel 10 degrees from the x-axis. This is a fairly large wave in deep water. The window width is 1 sec, 1/10 the zero crossing period, and the LFI solution is computed to third order (7 = 3). The LFI method has captured the location of the water surface and the kinematics at the surface essentially exactly. While linear wave theory might do an adequate job of approximating much of the kinematics of a deep water steady wave like this, it is important to remember that each of the points in Fig. 5.4 was computed from a small segment of the record surrounding that point. In this case, the window width is 1/10 of the zero crossing period, or ls. The local nature of this method extends its applicability to irregular wave records. Fig. 5.5 is a steady shallow wave generated by 18th order Fourier theory with the following parameters: wave height = 3 m, period = 10 seconds, water depth = 5 m, zero Eulerian current, and direction of travel 10 degrees from the x-axis. This is a fairly extreme wave in shallow water. The window width is 1 sec., 1/10 the zero crossing period, and the LFI solution is computed to third order. As with the deep water wave, the LFI method has captured the location of the water surface and the kinematics at the surface essentially exactly, including the pronounced sharp crest. Choice of Order In order to determine the order necessary to accurately capture the kinematics of measured waves, it is particularly useful to examine a window near the crest of a wave. The crest is usually the region that requires the highest order solution. This is particularly true for shallow water waves, but higher order wave theory in all depths of water indicates that, as the wave height increases, the crest tends to get sharper, and the trough flatter. Capturing this sharp crest requires a high order solution. Deep Water A window near a crest of the deep water steady wave given in Fig 5.4 has been computed at orders 1 through 4. The results in that window are given in Fig. 5.6 through 5.13. Figures 5.6, 5.8, 5.10, 5.12 are the non-dimensional errors in 93 Kinematics at the predicted water surface Figure 5.5: LFI predictions (J = 3) and exact kinematics at the center of the array at the predicted water surface for a steady shallow water wave 94 Crest of a deep water wave x10~ gauge | 5 0 -5 - -0.4 50-3 -0.2 -0.1 0 2 0.1 0.2 0.3 0.4 x 107 gauge 4 -0.4 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.6: Free surface boundary condition errors for a deep water wave at order 1 Crest of a deep water wave 0.24 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 0.1 — actual g 1 : 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.7: Free surface and velocities at the center of the array for a deep water wave at order 1 95 Crest of a deep water wave x10~ gauge 1 an : ; : : 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.8: Free surface boundary condition errors for a deep water wave at order 2 Crest of a deep water wave 04 -03 -02 -0.1 0 0.1 02 #038 04 0.1 — actual WW, 5 © predicted g 04 -03 -0.2 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.9: Free surface and velocities at the center of the array for a deep water wave at order 2 96 Crest of a deep water wave x10” gauge 1 04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.10: Free surface boundary condition errors for a deep water wave at order 3 Crest of a deep water wave 22 -04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 — actual © predicted 1 : : 04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.11: Free surface and velocities at the center of the array for a deep water wave at order 3 Crest of a deep water wave x10” gauge 1 Figure 5.12: Free surface boundary condition errors for a deep water wave at order 4 Crest of a deep water wave 04 -03 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.1 — actual WW, 6 © predicted g 04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.13: Free surface and velocities at the center of the array for a deep water wave at order 4 oi 98 the free surface boundary conditions (Eqs. 5.3 and 5.4) at each of the measured nodes. These are the data that would be analyzed in a practical situation. Figures 5.7, 5.9, 5.11, 5.13 are a comparison between the predicted and actual values for the water surface and the velocities at the surface at the center of the array. The actual values were computed using Fourier wave theory. The actual values would not be available for comparison when analyzing field records. The first order computation results in errors in the free surface boundary condi- tions of only order 10~° (Fig. 5.6) as well as very accurately predicting the velocities at the water surface at the center of the array (Fig. 5.7). It is a surprisingly accu- rate first order solution. This is because of the local nature of the method. When a single first order solution is used to capture the entire wave, the error is larger, of order 107°. It also should be noted that this first order solution is not the same as a linear wave theory solution, even locally. The full nonlinear boundary conditions are preserved, and the frequency and wave number are free to vary, and are not bound by a dispersion relationship. For steady waves in deep water, linear wave theory is fairly accurate. Linear theory is not, however, a local solution, and is not directly applicable to irregular waves. At second order, the free surface boundary condition errors are smaller, of order 10-®, and the velocities at the surface match the Fourier solution visually perfectly. At third and fourth order, the errors in the free surface boundary conditions continue to decline, and the water surface velocities continue to match the Fourier solution well. In deep water, for waves of this height, second order is more than adequate to capture the surface kinematics of this wave. Higher order solutions are likely to be necessary to capture the irregularity of field records, even in deep water. Choice of order is dictated by the desired accuracy of the solution, and by the ease of convergence to the solution. As there are more free parameters at higher order, a higher order solution will always have smaller errors in the free surface boundary conditions. In the case of this example, the solution converged at all orders very quickly, so there is little penalty in using third or fourth order. With an irregular record, in contrast, convergence can be more difficult, and it is occasionally necessary to resort to lower order. 99 As the order is increased, there are more free parameters, and care must be taken to include sufficient points in each window. As discussed above, for fairly low orders and an array of measurements, the limiting factor is the minimum number of points needed to define the curvature in the window, rather than to specify the system of equations. The above examples were computed from an array of three points, and so required a minimum of only one sample to specify the system at first and second order, and two points for up to order 8. When attempted, the solution would not converge with only one sample. With two samples, the optimization algorithm found a reasonable solution, but this small number of samples would not define the shape of the water surface adequately if it were part of an irregular record, and so three points were used in each window for all orders. In general, a minimum of three points should be used, and more may be necessary to accommodate a highly irregular profile, or the inevitable measurement error in a field record. In the case of theoretically generated records the need for more sampling of points poses no problem, but with field records, there are limitations as to the spacing of the sampled points. In order to free up the spacing of points for the LFI method, points are sampled from a cubic spline interpolation of the actual record. This allows the points to be sampled anywhere within the window. While computationally it is ' possible to sample as many points as necessary in a small window, if that window is, in fact, defined by only a couple of actual data points, it is not appropriate to try to fit a high order solution to a segment defined by only a few observational points. In order to include sufficient actual data points to justify the increased order, the window must be increased in size. While increasing the size of the window permits a higher order solution, it also compromises the local nature of the method. The goal of the LFI method is for the solution to be as local as possible, which is achieved by selecting as small a window as possible at fairly low order. Shallow Water A window near a crest of the shallow water steady wave given in 5.5 has been computed at orders 1 through 5. The results are given in Fig. 5.14 through 5.23. These figures are analogous to those previously discussed, but on a shallow water wave. Crest of a shallow water wave x107 gauge 1 04 303 -02 -0.1 (0) 0.1 x10 gauge 3 -04 -03 -02 -0.1 ie) 0.1 0.2 0.3 0.4 Figure 5.14: Free surface boundary condition errors for a shallow water wave at order 1 Crest of a shallow water wave — actual © predicted -0.2 A 4 a -04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.15: Free surface and velocities at the center of the array for a shallow water wave at order 1 100 101 Crest of a shallow water wave 3 gauge 1 x10 —0.01 5 1 -04 -03 -02 -01 0 0.1 0.2 0.3 0.4 Figure 5.16: Free surface boundary condition errors for a shallow water wave at order 2 Crest of a shallow water wave — actual -0. 02 © predicted WW» 0 4 eae) SAT Ae -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.17: Free surface and velocities at the center of the array for a shallow water wave at order 2 Crest of a shallow water wave x10~ gauge | -5 i -0.4 70.3 -0.2 -0.1 {0} 2 0.1 0.2 0.3 0.4 x10 gauge -2 i - 04 03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.18: Free surface boundary condition errors for a shallow water wave at order 3 Crest of a shallow water wave (0) -0. — actual 0.2 © predicted WW 0 t5 g -0.4 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.19: Free surface and velocities at the center of the array for a shallow water wave at order 3 103 Crest of a shallow water wave x10~ gauge 1 Z ; : 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.20: Free surface boundary condition errors for a shallow water wave at order 4 Crest of a shallow water wave 0.1 Tu ee: Ae ° Nw — actual © predicted -0.2 i wip 3 © of G2 @ ow Figure 5.21: Free surface and velocities at the center of the array for a shallow water wave at order 4 Crest of a shallow water wave x10" gauge 1 -1 4 i 04 -03 -02 -01 0 01 O2 O38 04 Figure 5.22: Free surface boundary condition errors for a shallow water wave at order 5 Crest of a shallow water wave -04 -03 -02 -0.1 (0) 0.1 02 03 04 Figure 5.23: Free surface and velocities at the center of the array for a shallow water wave at order 5 104 105 In this case, the first order computation results in substantial errors in the free surface boundary conditions of order 10~? (Fig. 5.14) and does not predict the ve- locities at the water surface at the center of the array as well as it did in deep water (Fig. 5.15), noticeably under-predicting the horizontal velocities. The local nature of the method allows a first order solution to get fairly close, but steady waves in shallow water are highly nonlinear, and a first order solution simply can not capture the sharp crest. At second order, the free surface boundary condition errors are slightly smaller, but still of order 10~?, and still under predicting the horizontal velocities at the surface. At third and fourth order, the errors in the free surface boundary conditions continue to decline, but the water surface velocities do not match the Fourier solution visually perfectly until fifth order (Fig. 5.23). In shallow water, it may be necessary to use up to fifth order to capture accurately the surface kinematics of a fairly large wave. In the case of this example, the solution converged at all orders very quickly, and there is little penalty in using fourth or fifth order. With an irregular record, convergence will be more difficult, and it may sometimes be necessary to resort to lower order. As previously discussed, the use of higher order solutions requires that more points ’ be sampled in each window. This is not likely to be a limitation with the single wave form with an array of measurements, as it is likely to be necessary to overspecify the system in order to define the curvature within the window. 5.4.2 Records of Two Intersecting Waves In order to develop the method using two intersecting waves, theoretical water surface records were generated by Ohyama’s fourth order intersecting wave theory (Ohyama, Jeng, and Hsu 1995a). Ohyama’s method is a Stokes-style solution for irrotational intersecting waves that is accurate to fourth order, and applicable in deep water (see Appendix A). As the resulting water surface from intersecting waves can be quite complex, a desired wave height can not be specified. Rather, the amplitude of the first order component of each of the intersecting free waves is specified. The 106 first order components are the largest components of the resulting wave field, and thus give an approximation of the size of the resulting waves. This wave theory assumes a zero Eulerian current. Standing Wave Figure 5.24 shows the results of the method for a standing wave generated by two identical intersecting waves traveling in opposing directions. The parameters of the waves are: Period = 10 sec., first order amplitudes = 3 m, prop- agation directions, 15 and 195 degrees from the x axis, in deep water. The LFI method to third order (J = 3) finds the kinematics at the water surface almost ex- actly. While these results show the complete wave, it is important to keep in mind that as with all the previous results, each of the indicated points is in the center of a separate window, and was computed independently of the other windows. In this case, the standard window width was 1 sec., or one tenth the zero-crossing period of the wave, with a third order potential function (J = 3), and four water surface nodes distributed equidistantly in time in each window (M = 4), at each gauge, for a slight overspecification, with 24 equations in 20 unknowns Short Crested Wave Figure 5.25 is the water surface at a point in time of the short crested wave that would result in the reflection of a steady wave from a sea wall, as indicated by figure 5.26. The parameters of the wave are: Period = 10 sec., first order amplitudes = 3 m, propagation directions, 30 and 150 degrees from the x axis, in deep water. Figure 5.27 shows the results of the method for this short crested wave. The LFI method to third order again finds the kinematics for this wave almost exactly. In this case, the standard window width was also 1 sec., or one tenth the zero-crossing period of the wave, with a third order potential function (J = 3) and four water surface nodes (M = 4) in each window, for a slight overspecification, with 24 equations in 20 unknowns. Choice of Order As with the previous examples with steady waves, it is useful to examine a window near the crest in order to determine the order necessary to accurately capture the 107 Kinematics at the predicted water surface — Ohyama et al. | © (al Figure 5.24: LFI predictions (J = 3) and exact kinematics at the center of the array at the predicted water surface for a standing wave 108 ip Ye Ss Vises tL -5 CYL PEERS SOS SGP SSS ~10 =“ USGS 300 300 0 y o x Figure 5.26: Schematic of reflection that would generate the short crested wave given in figure 5.25 kinematics of a short crested wave. A window near a crest of the short crested wave given in Fig. 5.27 has been com- puted at orders 1 through 4. The results in that window are given in Fig. 5.28 through 5.35. Figures 5.28, 5.30, 5.32, 5.34 show the errors in the free surface boundary con- ditions at each of the measured nodes. These are the data that would be analyzed in a practical situation. Figs. 5.29, 5.31, 5.33, 5.35 give a comparison between the pre- dicted and actual values for the water surface and the velocities at the surface at the 109 Kinematics at the predicted water surface —— Ohyama et al. Oo LFI Figure 5.27: LFI predictions (J = 3) and analytical kinematics at the predicted water surface for a short crested wave 110 Crest of a short crested wave x 107 gauge | 04 -03 -02 -0.1 to) 0.1 0.2 0.3 0.4 Figure 5.28: Free surface boundary condition errors for a short crested wave at order 1 Crest of a short crested wave 0.4 -03 -02 -01 0 0.1 02 O38 04 f — actual WW, 3 © predicted 04 -03 -02 -0.1 te) 0.1 0.2 0.3 0.4 Figure 5.29: Free surface and velocities at the center of the array for a short crested wave at order 1 Crest of a short crested wave x107 gauge 1 -1 04 ~<03 -02 -01 x10 a 2 ° 04 <03 -02 -0.1 x10 ~ : : 04 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.30: Free surface boundary condition errors for a short crested wave at order 2 Crest of a short crested wave 04 -03 -02 -0.1 0 0.1 02 O38 04 OF -0.1 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.31: Free surface and velocities at the center of the array for a short crested wave at order 2 © predicted ILL Crest of a short crested wave x10” gauge | 10 5 (0) —-5 i -0.4 50.3 -0.2 -0.1 (0) 2 0.1 0.2 0.3 0.4 x 107 gauge -1 -0.4 -03 -02 -0.1 (0) 0.1 0.2 0.3 0.4 Figure 5.32: Free surface boundary condition errors for a short crested wave at order 3 Crest of a short crested wave 04 03 -02 -01 0 0.1 0.2 0.3 0.4 — actual O° predicted 04 -03 -02 -01 0 0.1 02 #03 04 Figure 5.33: Free surface and velocities at the center of the array for a short crested wave at order 3 113 Crest of a short crested wave x10" gauge 1 -04 -03 -02 -0.1 ie) 0.1 0.2 0.3 0.4 Figure 5.34: Free surface boundary condition errors for a short crested wave at order 4 Crest of a short crested wave -0.1 i -04 -0.3 -0.2 -0.1 0 O41 O02 O38 04 Figure 5.35: Free surface and velocities at the center of the array for a short crested wave at order 4 114 center of the array. The theoretical values used Ohyama’s fourth order intersecting wave theory. The actual values would not be available for comparison when analyzing field records. It is clear that the first order computation results in substantial errors in the free surface boundary conditions of order 107° (Fig. 5.28), as well as significantly underestimating the horizontal velocities (Fig. 5.29). As with the steady deep water wave previously discussed, the first order solution is quite reasonable. This is because of the local nature of the method, and the fact that it is not a linear solution. The full nonlinear boundary conditions are imposed, and the frequencies and wave numbers are free to vary, and are not bound by a dispersion relationship. At second order, the free surface boundary condition errors are better, of order 10-°, and the velocities at the surface match the analytical solution visually perfectly. At third and fourth order, the errors in the free surface boundary conditions continue to decline, and the water surface velocities continue to match the theoretical solution well. As with the single wave method, choice of order is dictated by the desired ac- curacy of the solution, and by the ease of convergence of the optimization. For the above examples, four water surface nodes (M = 4) distributed equidistantly in time in each window, at each gauge, at orders 1 through 3. These values providing an over- specification and sufficient points to define the shape of the water surface throughout the window. Five points (M = 5) were used at fourth order, as four points provide only 24 equations, and there are 28 parameters to be found. Five points provides 30 equations for a slight over specification. At fifth order, there are 38 unknowns, and seven points in time (M = 7) would have to be used in each window. Sampling this many points in a single window would require very closely spaced data or a wider window. Another solution would be to use an array with additional measurement locations. An array of four gauges, for example, would provide eight equations per point in time, and would allow a fifth order solution to be computed from five points (M = 5) per window. It would be advantageous to use as many gauges as possible in shallow water, where higher order solutions are necessary. Unfortunately, the theoretical solution used for this analysis, while the best avail- 115 able for intersecting waves, is only accurate to fourth order, and only being applicable in deep or transitional water. This precluded an analysis of the behavior of the so- lution in shallow water. When the method is applied in shallow water, the lessons learned in the two dimensional case should be applied, and a higher order solution, perhaps to fifth order, should be used. It is also the case that the Ohyama solution is a single solution that is accurate to fourth order globally. The local nature of the LFI method allows each separate window in time to be represented by a unique solution, representing the entire record with far more accuracy than a global solution of the same order. Choice of Window Width As mentioned in the previous section, it is important to keep the window width at a minimum. For the short crested wave shown in Figure 5.27, the standard window width was one tenth the zero-crossing period (0.17,). It was necessary to widen the window to 0.157. at three locations in order to find a solution. This indicates that the standard window width is set close to the smallest size that would be effective. In general, it is necessary for the window to contain enough of the record to define the local curvature, and include sufficient data points to reasonably expect to define a high order solution. Most often, the limiting factor will be the sample rate of the data. For example, the above examples are computed for a short crested wave with a period of ten seconds. In those examples, the window width of 1/10 the period (1s) was adequate to capture most of the wave. Unfortunately, field records are often sampled at only 1Hz. At this sampling frequency, and a one second window, there would be a maximum of only two data points in each window, and it may be necessary to widen the window to include more data. Effects of Array Size The examples above are all computed from data sampled by a fairly small array (Fig. 5.3). The 1.6m size of the array is about 0.01 of the wavelength of 10s wave in deep water. That particular size was chosen because it has the same dimensions as 116 the DWG-1 pressure array that is used for the analysis of pressure records in Ch. 6. A small array has the advantage of having a small distance to interpolate to the center of the array, and is local in the sense that each of the sensors is in essentially the same part of the wave. The disadvantage of a small array is that the differences between the measured values at each of the nodes are small, and thus more sensitive to measurement error. With the theoretical records used for the previous analysis, the method was not sensitive to array size, giving equally accurate results with arrays up to ten times the size of the DWG-1. The spacing between the gauges of an array will determine the wavelengths to which the gauge is most sensitive. Window widths of approximately one tenth a typical period have been found to be effective. This indicates that a gauge spacing of about one tenth of the expected dominant wave length might be an optimum spacing. 5.5 Laboratory Measurements The results from the analysis of analytically derived records given in the previous sections demonstrate the method’s capabilities in a variety of conditions. They have also been very useful in helping to determine the range of solution parameters that must be adopted, including order of solution, window width, and the number of time samples taken in each window. The question still remains, however, as to how well the method works with actual irregular waves. 5.5.1 Laboratory Experiment The following results are taken from a laboratory experiment performed together with Dr. Steven A. Hughes at the Coastal Engineering Research Center, Army Corps of Engineers Waterways Experiment Station in Vicksburg, MS, in June of 1995. The waves were generated in a 0.46m wide flume by a programmable, hydraulically driven, piston type wave board (see Fig. 5.36). At the other end of the flume was a beach with 1:30 slope, and a single layer of wave-absorbing horsehair matting. This beach has been shown to have bulk reflection coefficients that vary from about 5% for IOLA WG Y y y y y Z oo 5 ; At Plan Beach with horsehair matting Wave board 1 a Uj WW Elevation (not to scale) Figure 5.36: Schematic of wave flume Hp»0/h = 0.1 up to about 15% for Hmo/h = 0.4 (where Hmo is the incident spectral significant wave height and h is the water depth) (Sultan 1992; Sultan and Hughes 1993). The water surface was measured by capacitance wave rods, calibrated with a cubic calibration function. The velocity data were collected using a Dantek laser Doppler velocimeter (LDV) system operated in the backscatter mode. The LDV system fea- tured a 2-watt argon-ion laser equipped with a fiber-optic probe that measures two orthogonal water velocity components (horizontal and vertical). Velocity data were converted in real time to engineering units (m/s) and written to a computer file simultaneously with the wave rod data. The wave rods and LDV were placed near the middle of the flume, with the wave rods arranged in an equilateral triangle with leg length of 0.17m (see Fig. 5.37). The LDV was situated to measure the horizontal and vertical velocities at the center of the array, at a variety of vertical elevations. The flume was allowed to reach quiescence in between runs, and the waves were measured for only a short time after starting the 118 e Wave maker ———~> — LDV location 0.17m gauge 3 gauge 1 Plan Figure 5.37: Layout of array in flume wave maker. At a depth of 0.35m, the approximate wave celerity (\/gh) is 1.8m/s. As the measuring location is located 9m from the beach, it takes about 10s for the waves to travel from the measurement location to the beach, and for the reflected waves to travel back to the study location. Thus, the first 10s of the record are assured of being uncontaminated by reflection from the beach. Using only the beginning of the record also prevents a Stokes drift induced return current to be established, and so the Eulerian current is taken to be zero. 5.5.2 Laboratory Results Figure 5.38 shows a sample wave taken from an irregular wave record measured in the flume. The waves were generated with a target TMA spectrum (Bouws et al. 1985) with a peak period of 2.5s and mean wave height of approximately 12cm. The water depth was 0.35m, with the velocities measured by the LDV at an elevation of 0.05m below the mean water level. The solid lines on figure 5.38 are the horizontal and vertical velocities as measured by the LDV. The circles are the velocities predicted by the LFI method at the center of each computed window. The solid line in the water surface plot is the measured elevation at wave rod 2, at the same x coordinate as the center of the array. The circles are the water surface elevations at the center of the Figure 5.38: Water surface Results on a sample wave in a flume 119 -3 node # 1 70.4 -03 -02 -01 01 O2 O38 04 node #3 —P ; 04 -03 -02 -0.1 0 0.1 0.2 0.3 0.4 Figure 5.39: Free surface boundary condition errors at the crest of a sample wave in a flume Figure 5.40: Computed and measured velocities near the crest of a sample wave in a flume 120 array as computed by the LFI method. This solution was computed to third order, with a window width of 1/10 the mean zero crossing period of the record (approximately 0.2s). Four time samples were taken in each window, providing substantial overspecification of the system, and allowing the shape of the water surface to be well defined in each window. The first trough is fairly flat, and the optimization did not converge to a solution quickly with the small window. After the window was widened, the optimization converged quickly to the given results. The LFI method has captured the detail of the kinematics of this irregular wave very well. The comparisons to the measured data are given just below the troughs of the wave. Note that the LFI method accurately captured the secondary hump in the vertical velocity (w) near the 2.5s point. The data used to compute the kinematics were the measured water surface. The computed water surface was found from the predicted kinematics, and it matches the measured water surface almost exactly. The accurate computation of the water surface indicates that the predicted kinematics near the surface may, in fact, be more accurate than the predictions at the elevation where the kinematics were measured. This is to be expected, as the data used to compute the solution was measured at the surface, and thus the results should be - most accurate there. Achieving the greatest accuracy near the surface can be useful, as the surface is the location where the water velocities and accelerations are greatest. Figures 5.39 and 5.40 show the results of the computation in a window near the crest of the wave given in Fig. 5.38. The errors in the DFSBC are of order 107%, and substantially smaller in the MKFSBC. The computed and measured velocities do not match exactly, but the magnitude and trend are very close. The results are similar in the other windows. 5.6 Discussion The LFI method is an effective and efficient way to re-create the kinematics of irregular waves from the measurements of an array of wave rods. It is able to capture the kinematics from primarily uni-directional seas as well as the seas near a reflect- 122 ing surface. The local nature of the approach allows a nonlinear solution without prohibitive computational costs, using a fairly simple form for the potential function, and allowing the parameters of the potential function to change with time. The examples in this Chapter were all computed using MATLAB running under the Linux operating system on an Intel 9)MHz Pentium processor based PC. The shallow water steady wave (Fig 5.5) took the longest to compute, about 30min, and 85x 10° floating point operations (flops). The steady deep water wave (Fig. 5.4) took the least time, 10min. and 15x10® flops, the standing wave (Fig. 5.24) took 20min. and 244x10® flops, the short crested wave (Fig. 5.27) took 11min. and 93x10® flops, and the laboratory record (Fig. 5.38) took 21min. and 83x10® flops. The examples given in this chapter provide guidance as to the parameters of a solution to be applied to field records. Far from reflecting surfaces, using a potential function representing a single wave is effective. Near a reflecting surface, a potential function representing two nonlinear intersecting waves is capable of capturing the standing or short crested waves that are likely to develop. Window widths of 1/10 of the zero crossing period are small enough to maintain the local nature of the solution, and capture the detail of the record, while being large enough to include the local trend of the wave field. On certain segments of the record, the window width must be increased in order to include sufficient curvature in the record to find a solution. The widening of the window is most often needed in long, flat troughs in shallow water, or near zero crossings of the record. In either case, the window rarely needs to be larger than 1/5 of the zero crossing period. Low order solutions are quite adequate for deep water waves. For most waves, second order is adequate. For the very steepest waves, slightly higher order may be appropriate, and there is little computational penalty in including the higher order terms. In shallow water, higher order solutions are necessary. Third order is adequate in many cases, but including up to fifth order is recommended for extreme waves. When attempting a high order solution such as this with the two wave method, it may be necessary to use arrays with more than three gauges in order to provide sufficient equations to specify the solution. 123 Chapter 6 Array of Subsurface Pressure Measurements Arrays of subsurface pressure gauges are installed in the sea in order to capture directional wave field data. These instruments have an advantage over arrays of water surface wave rods, as they are mounted below the surface, where they are less susceptible to damage, by either the strong wave forces experienced during storms, as well as vandalism or accidents with vessels. Unfortunately, the reason they are less susceptible to damage from wave forces is because they are placed under the surface, often at the sea bed, where the action of waves is the smallest. This leads to the mathematically ill-posed problem discussed in chapter 3. This difficulty is somewhat mitigated by the added information made available by the multiple gauges in the array, but still limits the ability to determine the detail of the kinematics, particularly when there is a lot of high frequency energy in the sea state. 6.1 Formulation of Solution The lessons learned in the previous chapters provide for a fairly straightforward determination of the formulation for adapting the LFI method to the analysis of arrays of pressure gauges. The formulation for a single subsurface pressure gauge was 124 presented in chapter 3, and for arrays of water surface measurements in chapter 5. The formulation for arrays of subsurface pressure gauges is very similar to that for arrays of water surface measurements, with the addition of the need to determine the location of the water surface, as was done with the analysis of a subsurface pressure trace. The following description includes the required information; greater detail has been given in previous chapters. The flow is assumed to be irrotational and incompressible, with a potential func- tion that represents either a single directional wave, or two separate intersecting waves. The potential function locally representing a single wave is Eq. 4.10 reduced to a single wave: cosh 7A (h + z) cosh j Ah k= eae where U, and U, are the components of the known depth uniform Eulerian current i (x, z,t) =Urzr + Uyy + > A; j=l sinj(kzr +kyy—wt+a) (6.1) in the z and y directions, h is the mean water depth, J is the truncation order of the Fourier series, A; are the Fourier coefficients, w is the local fundamental frequency, k, and k, are the components of the local fundamental wave number in the z and y directions, and A is the magnitude of the local wave number. The potential function for two intersecting waves to fourth order is Eq. 4.18, repeated here: 20 oz, Y, 2; t) = U;,z + Uyy ale S- A; C(K;h, K;z) sin 0; (6.2) t=] where: 3 : cosh h;(h + z) 0c;° 00; K;h, Kz = SST Gg= =—— pee (K;h, K;z) eon an K ae ait a Si = (sane + ky iy — wt + a) So = (kr 2t + ky sy — wot + a2) At first order: a1 = (51) At second order: G3 = (2Sa)) 05 = (5; + S2) At third order: 07 = (35) 09 = (2S, + S2) oy, = (Si + 282) And at fourth order: o13 = (451) o15 = (35; + S2) O17 = (25; + 25>) o19 = (51 + 352) 125 oD = (S2) 04 = (252) 06 = (Sy — 52) 0g = (352) 010 = (2S, = 52) 012 = (S) = 2S) Both of these potential functions solve the field equation (Laplace, Eq. 4.2) and the bottom boundary condition (Eq. 4.3) exactly. The rest of the solution is determined by the free surface boundary conditions: the modified free surface boundary condition Ga ao fF ==> + 9wt 2% Ot Ow (6.3) 126 the dynamic free surface boundary condition (f?), Oe aa J Sata ea era F=0 BM 2S i 2) the Bernoulli equation (f”), Od 1 Pa = fF =—+-(v’ tv? +w’?)+—-B=0 6.5 ot | ) Pp \ I BWR = nal 1 2 2 A;K; Sati D i) + (= “| () w=1 In order to apply the free surface boundary conditions, the location of the free surface must be found, together with the potential function, in each window. In order to locate the free surface, the water surface is defined at M nodes, located at the horizontal location of the center of the array, and equidistantly spaced in time throughout the window. The elevation of these nodes is unknown, and will be sought as part of the solution. 6.2 Formulation of the Optimization The formulation of the optimization for the LFI method as applied to the inter- pretation of an array of pressure measurements has a great deal in common with the formulation for an array of water surface measurements and the formulation for a sub- surface pressure record (chapters 3 and 5). Only the framework is presented, together with any additional information unique to the application to a pressure array. Observational Equations The observational equations are the equations that con- fine the solution to fit the measured records. The predicted subsurface pressure is available from the Bernoulli equation, Eq. 6.5, which is applied at the location of the gauges in the array to define the observational equations. The error in the Bernoulli equation is the difference between the measured dynamic pressure at the given lo- cation and time and the dynamic pressure predicted from the kinematics defined by TG the potential function. Sufficient independent equations are defined by applying the Bernoulli equation at a number of points in time throughout the considered window. In order to specify the solution to a system of equations, there must be at least as many independent equations as unknown parameters in the system. In order to specify the solution to the LFI formulation for a pressure array, the boundary conditions (f? and f*) are applied at each of the water surface nodes, yielding 2M B.) is applied at J times on each of the N n,t equations, and the Bernoulli equation ( measured pressure records within the local window, yielding NJ equations. In the single wave case, there are 4+ J + M unknown parameters sought in Eq. 6.1 (kz, ky, w, a, Ay... Aj, and 7...) in each window, so that if (2M@+N/)>(4+J+M) the solution is specified. In the two wave case, there are 2)> 7 +8+M unknowns in Eq. 6.2 (2k,, 2ky, 2w, 2a, Ay... Az, and m...nvw: 10+ M at first order, 14+ M at second order, 20+ M at third order, and 28+ M at fourth order). The solution is specified when (2M + NI) > (25>7 +8+M). This system results in the following nonlinear least squares optimization. minimize O(X yr DO aa gin ae n=1 t=1 (6.7) SPX: Ley Ye, "Mm; aa) ap fn (Xs Le, Ye, Mm, bia) m=1 where: 28 = (Re, py Oy Oh Bigg coe g Agi o0- Mn) for the one wave case, and: DC (eet pil she Cig [eobin | Way Ob, As oo guava Mil os olfhve) for the two wave case. ,, Yn, and z,,p are the coordinates of the nth gauge, x. and y. are the coordinates of the center of the array, and 77,, is the elevation of water surface node m. As with the previous analysis, overspecification is helpful in accommodating the measurement errors in field records, as well as being required to define the shape of measured records and water surface in each window. 128 6.3 Computation Methods The LFI method for an array of pressure measurements can be broken down into the same sequence of steps as with the analysis of a water surface array: 1. Pre-processing of record. a) Determine estimate for level of noise in the record. b Determine MWL and subtract hydrostatic pressure from the records. d (a) (b) (c) Determine estimate for magnitude and direction of Eulerian current. (d) Define a set of continuous records from each gauge from the discrete ob- servations by cubic spline interpolation. (e) Specify spacing of output locations. (f) Compute mean zero crossing frequency. (g) Non-dimensionalize record and all parameters. 2. Primary values of numerical solution parameters are chosen. (a) Window width (79) (b) Order of solution (J) (c) Number of time samples of the pressure records within each window (1) (d) Number of water surface nodes within each window (M) 3. Global solution is computed on an entire wave to provide first estimates for local optimization 4. For each selected output location, a window in the record is defined, and an LFI solution is computed. (a) Initial guess for the optimization determined from the global solution. (b) Full nonlinear optimization for all unknown components of the potential function is computed. (c) Results are checked for spurious solution. i. If no solution, or a spurious solution, is found, the solution parameters are adjusted, and the optimization repeated. ii. If a good solution is found, progress to the next window. 6.3.1 Pre-Processing of Record The pre-processing of the measured records is essentially the same as presented in the previous chapters. Estimates for the mean water level (h) and Eulerian current (U, and U,) must be computed to define the propagation medium. Cubic splines of the measured records are computed to provide continuous records for computation, and the desired output locations are chosen. The records and all parameters and variables are non-dimensionalized by the following parameters: the mass density of water (p), acceleration of gravity (g), and the mean zero crossing frequency of the measured records (w,). 6.3.2 Optimization Procedure The nonlinear optimization in the LFI method for the analysis of arrays of pressure measurements is very similar to that used for the analysis of a single subsurface pressure trace. In the single wave approach, the solution is somewhat easier. The potential function used by the single wave approach is almost the same as that used with a point measurement (both chapter 3 and Sobey (1992)), with the addition of a directional component to the wave number. This results in a single additional unknown, but the array of measurements provides at least three times as much data, with gauges at at least three spatial locations to specify uniquely the directional structure of the sea. The result is that the optimization tends to converge more rapidly and robustly than with a single point measurement. When applying the method with the two wave potential function, there are far more unknowns, and the optimization once again becomes somewhat tenuous. As with the previous chapters, it is essential to identify a good initial estimate to reduce 130 the chances that the optimization will converge to a spurious minimum. Global Solution As with arrays of water surface measurements, a single short segment of the mea- sured pressure records often may not contain sufficient information to identify the directional trend of the wave field. As a result, it is necessary to examine a larger segment of the record, for example an entire wave from crest to crest or trough to trough, to get a general sense of the directional trend of the wave field. This global solution can then be refined to fit a small segment very closely. This approach is used to establish the initial estimate for the optimization procedure in each window. The procedure for computing the initial estimate to the global solution for an array of pressure measurements is virtually identical to that used in the previous chapter for arrays of water surface measurements. It will be outlined here. For more detail, see section 5.3.2. For the initial estimate, it is assumed that the pressure records can be approxi- mated with linear wave theory: Pa = acos(kzt, + kyYn — wtm + a) (6.8) for the single wave method, and Pa = 4 C0S(ke itn + ky Yn — Witm + A) (6.9) +a2 C08 (kz,22n + ky,2Yn — Watm + Q2) for the two wave method, where: an = AnpWn coy Za as *») cosh kh (6.10) Ky, = \/k?2 + k2 A, is the amplitude of the linear potential function of the nth wave, and 6. = angle (x SS 3. where: Oa=mr< angle( Dm) <6 Frequency The frequency is computed from the second derivative of the pressure at the center of the array. 1 O-pitem ee 6.14 Pd,cjm ot? ( ) Wm = 8? Da.c.m at? wave method, it is assumed that the bimodal sea is the result of reflection, so that where is computed from the same smoothed spline used above. For the two the same frequency is used as the first estimate for both waves. Wave Numbers The wave numbers are estimated from the linear dispersion rela- tion. (w — K,U,)? = gK, tanh K,h, [Big cy IS COS Oa. OG sim ae (Gell) 132 where U,, is the component of the Eulerian current in the direction of the nth wave. U, = ,/U2 + U2 cos (can (=) - é) (6.16) Amplitudes and Phases The amplitudes and phases the record are found by separating the cosine and sine components so that the equations can be solved as a linear least squares problem. Dac = acos(k,x + kyy — wt + a) = bcos (k,z + kyy — wt) + csin (k,x + kyy — wt) (6.17) G=ViPEe a@ = arctan (—c/b) for the single wave method, and - Dac = 4,008(ke it + kyiy — wt + a1) + a2 008 (ker + ky oy — wt + a2) = 6, cos(kzix + kyiy — wt) +c, sin (kz iz + ky iy — wt) +-b2 cos (kz,22 + ky,2y — wt) + co sin (kz2r + ky2y — wt) Th = b? + c? (6.18) a= /+c Q, = arctan (—c,/b;) Qa, = arctan (—c2/b2) Refining the Linear Estimates These rough estimates are refined by optimizing for a best linear wave theory fit to the record: N M minumize O(X) = DS Ni (ives 3 22 OS Monta) (6.19) (=) ToSil where: fP4 (X35 Ln, Yn, tm) = acos (kptn + kyYn — wtm + @) X = (kz, ky, a, a) K=,/k2+k 133 for the single wave method, and: fP4 (X35 Ln, Yns tm) = a1 008 (kein + ky i Yn — Witm + 1) + a2 608 (kz22n + ky,2Yn — Watm + Q2) xX = (ke, kya, 1, kz2, ky2, a2, 41, a2) Kn = [Rn + hey for the two wave method. The frequencies are determined from the linear dispersion relation: wn = /gK, tanh K,h + K,U, (6.20) where U,, is defined by Eq. 6.16 This optimization results in a linear estimate for the dynamic pressure that fits the measured records most closely in the segment considered. The final step in computing the initial estimates for the final window-by-window optimization is to compute a full order global solution to this large segment of the record. The initial parameters for this full order global optimization are provided by the computed linear fit to the water surface records, with the amplitudes adjusted for ' the potential function: An cosh h,,h 7 poh coshlka(h-= 25) (6.21) The higher order Fourier amplitudes are all initially set to zero. Water Surface The location of the water surface at M nodes in the center of the array throughout the segment is estimated from the linear pressure response function with stretching (Nielsen 1989): a, (tm) cosh (k (h + pa(tm)/pg)) Wiphtre = RETR Aaa = (6.22) pg cosh k(h + zp) where k is the mean of the two wave numbers, Zp is the elevation of the pressure array, and 7, is the elevation of the water surface node at the center of the array at tm- Pa(tm) is computed from Eq. 6.8 or Eq. 6.9. 134 Using these linear wave theory estimates as the first guess, the full order optimiza- tion, Eq. 6.7, is computed to determine the best full order fit to a global segment of the record. The parameters of the potential function computed by this optimization are then used as the initial estimate for the final optimization in each defined small window in time. Phase Shift The phase parameters, a; and a2 are adjusted to accommodate the change in the time reference frame from the global to the local solution: On = —w,At+ Qn,global (6.23) where At is the difference between the time in the two reference frames. Nonlinear Optimization in Windows The established global solution is used as the first guess for the parameters in each window, and these parameters are refined in the same manner as the previous chapters by solving Eq. 6.7 with any standard nonlinear optimization routine. The resulting solution is then checked to see if is clearly spurious. Spurious solu- tions can be identified by the same criteria used in the previous chapters: e Very large or systematically variable errors. e First order amplitude smaller than one of the higher order amplitudes. e Unrealistically large or small frequency or wave number. e Large discontinuity between the windows in the predicted kinematics. If no solution or a spurious solution is found, the parameters of the solution are revised for another attempt. These nevisllons include increasing the width of the window, and decreasing the order of the solution. If neither of these adjustments result in an acceptable solution, the window is skipped, and future analysis must be interpolated through that point. As with the analysis of a point pressure measurement, these adjustments are most likely to be needed in the long, flat troughs of shallow water waves, or near zero 135 crossings in the record. The additional data provided by the array of measurements provides for a more robust optimization than with a single subsurface pressure mea- surement. As well as providing additional data, the spatial array provides information about the evolution of the wave field in space, helping to define the local wave number more clearly. This is in contrast to the analysis of a point measurement, where the spatial evolution of the wave field must be determined entirely from the measured temporal evolution, coupled with the governing equations. As a result, adjusting the solution parameters is necessary less often than with a single measurement. 6.4 Theoretical Records 6.4.1 Single Wave Records The following results are from “measured” data generated by Fourier steady wave theory (Sobey 1989), providing a near-exact solution for steady waves that provides the complete kinematics. The data used are a simulation of data that might be col- lected by a DWG-1 pressure array (Howell 1992). The DWG-1 is a reliable, easy to deploy pressure array. The unit is capable of including up to six independent pressure transducers, but is frequently used with three transducers, arranged in an 1.6 m equi- lateral triangle (fig. 6.1). Three measurements were chosen for the theoretical data, as three is the minimum number required to provide directional information. Additional measurements would provide overspecification, and can easily be accommodated in the formulation. Shallow Water Figure 6.2 is a steady shallow wave generated by 18th order Fourier theory with the following parameters: wave height = 3m, period = 10s, water depth = 5m, direction of travel 10 degrees from the x-axis, and an opposing Eulerian current of 2m/s. This is a near limit wave in this depth water, with this opposing current. The pressure array is located at the bottom. The parameters of the LFI solution are: window width = ls (7 = 0.1T,), fifth order (J = 5), with 6 samples on the water surface, and 6 samples on each of the pressure records (M = J = 6), resulting in 30 136 wave ray Pressure Gauges Figure 6.1: Layout of DWG-1 pressure array equations in 15 unknowns. The solid lines are the dynamic pressure, water surface, and kinematics as com- puted from the Fourier solution. The circles are the predictions of the LFI method, with each circle being located at the center of a separate window. The kinematics are measured at a depth of 1.5m below the surface, just under the trough. The LFI method has captured the location of the water surface and the kinematics essentially exactly, including the pronounced sharp crest. As with the single pressure gauge (Sec. 3.4), a fairly high order solution (J = 5) was necessary to capture the sharp crest of this shallow water wave. The additional data provided by the three gauges allowed for a narrower window than with the single measurement, even at this high order. Deep Water Figure 6.3 is a steady deep water wave generated by 10th order Fourier theory with the following parameters: wave height = 10m, period = 10s, water depth = 100m, direction of travel 10 degrees from the x-axis, and zero Eulerian current. The pressure array is located 10m below the surface. The parameters of the LFI solution are: window width = 1s (7 = 0.1T.), second order(J = 2), with 3 samples on the water surface, and 3 samples on each of the pressure records (M = J = 3), resulting Kinematics at the center of a pressure array — Fourier Oo LFI tw. Figure 6.2: LFI predictions (J = 5) and exact kinematics at the center of the array at z = —1.5 m for a steady shallow water wave. 138 Kinematics at the center of a pressure array Figure 6.3: LFI predictions (J = 2) and exact kinematics at the center of the array at z = —5 m for a steady deep water wave. 139 in 12 equations in 9 unknowns. The LFI method has captured the location of the water surface and the kinematics essentially exactly, at only second order. Higher order may be necessary to define the kinematics of a more irregular record, and can be included with little computational penalty. 6.4.2 Records of Two Intersecting Waves Theoretical pressure records were generated by the same fourth order stokes-type intersecting wave theory used in the previous chapter (Ohyama, Jeng, and Hsu 1995a). This method provides the complete kinematics, is applicable in deep water, and as- sumes a zero Eulerian current (see appendix A). Standing Wave Figure 6.4 shows the results of the method for a standing wave. The parameters of the Ohyama solution are: period = 10s, first order amplitudes = 3m (resulting in a total wave height of slightly over 6m), propagation directions = 15 and 195 degrees from the x axis, in deep water. The pressure array is located 10m below the mean water level. The LFI method applied to third order (J = 3) finds the the water surface and the kinematics at the elevation of 3m below the surface, just below the trough of the wave. While these results show the complete wave, it is important to keep in mind that as with all the previous results, each of the indicated points is in the center of a separate window, and was computed independently of the other windows. In this case, the standard window width was 1s (7 = 0.1T.), with a third order potential function (J = 3), four water surface nodes (M = 4) and six samples on the pressure records (J = 6) distributed equally in time in each window, resulting in 26 equations in 24 unknowns. Short Crested Wave Figure 6.6 shows the results of the method for the short crested wave shown in Figure 6.5 The parameters of the wave are: period = 10s, first order amplitudes = 3m (resulting in a wave height of just over 6m), propagation directions: 30 and 150 degrees from the x axis, in deep water. The pressure array 140 Kinematics at the center of a pressure array Figure 6.4: LFI predictions (J = 3) and exact kinematics at the center of the array at z = —3 m for a standing wave. 141 300 Figure 6.5: Water surface of short crested wave at t=0 is located at 10m below the MWL. The LFI method applied to third order (J = 3) again finds the kinematics for this wave almost exactly. In this case, the standard window width was also ls (7) = 0.17.), with four water surface nodes (M = 4) and six samples on the pressure records (J = 6) distributed uniformly in time in each window, resulting in 26 equations in 24 unknowns. Figure 6.7 is the same short crested wave, but computed with the single wave method of the LFI solution. In this case, the LFI solution still matches the measured pressure record very well, as is virtually always the case, as those measurements are part of the optimization. The predicted water surface is also fairly accurate, but with the predicted crest slightly underestimated. The vertical velocity is also fairly accurate. The LFI predictions for horizontal velocities, on the other hand, are very different from the Ohyama solution. The resulting discontinuities in the predictions for the horizontal velocities make it clear that something important is missing from the solution. The single wave method is simply not able to capture a short crested sea made up of two distinct directional components. Kinematics at the center of a pressure array 0 1 2 3 4 5 6 7 8 9 — Ohyama etal. Oo LFI Figure 6.6: LFI predictions (J = 3) and exact kinematics at the center of the array at z = —3 m for a short crested wave. 143 Kinematics at the center of a pressure array la Figure 6.7: LFI predictions (J = 3) and exact kinematics at the center of the array at z = —3 m for a short crested wave computed with the single wave method. 144 Choice of Order The choice of order for the analysis of arrays of pressure measurements is essen- tially the same as for a single pressure measurement. See chapter 3 for a detailed discussion. In general, lower orders are necessary for deep water than shallow water, with higher than third order unlikely to be necessary in deep water. Shallow water may require up to fifth or sixth order in order to capture near limit waves. Choice of Window Width The choice of window width is also very similar for the analysis of pressure arrays as it is for the analysis of a single pressure measurement, or arrays of water surface measurements. Window widths of a minimum of about 1/10 the zero crossing period are required, with the occasional need to widen the windows near zero crossings. The additional data provided by the array of pressure measurements allowed a solution to a very steep shallow water wave at high order with a window width of 0.17,, where a similar solution required the window width to be doubled to 0.27, when there was only a single point pressure measurement available (Section 3.4). Number of Nodes on the Water Surface and Pressure Records The criteria developed in the previous chapters for the number of water surface nodes are applicable here. J+1 points are required to specify uniquely a Fourier series of order J. In order to keep the order of the water surface consistent with the order of the potential function, at least M = J+1 water surface nodes should be used. The same number of points on the measured pressure records is usually adequate. For the standing and short-crested waves, additional equations were needed to specify the solution, so J = J + 3 points on the pressure records were used. The results are not highly sensitive to this parameter, provided sufficient samples are used to define the local curvature and specify the system of equations. More points on the measured pressure records will assist in accommodating noise in the record. 145 6.5 Field Measurements The following results are from data provided by Gary L. Howell of the US Army Corps of Engineers Waterways Experiment Station, Coastal Engineering Research Center (CERC). The data were collected as part of CERC’s Coastal Inlets Research Program. They were collected near the Ponce de Leon inlet, Florida, by a DWG-1 with three absolute pressure gauges, data sampled at 5 Hz. The DWG-1 was resting on the bottom, so that the pressure sensor diaphragms were positioned 0.21m from the bed. 6.5.1 Field Results Figure 6.8 is a segment of a record collected on August 20, 1996. The statistics of the record, based on about one hour of record are: mean wave height = 1.05m, peak period = 5.6s, peak direction = 83° from true north, in a depth of approximately 7.9m. The top plot is the dynamic pressure at the center of the array. The solid line is the approximate measured pressure, computed by linear interpolation from the three measured points in the array. The circles are the values predicted by the LFI method. The other four plots are the predicted water surface, and the three orthogonal velocities, as predicted by the LFI computation. The single wave method was used, with the following parameters: third order (J = 3), four water surface nodes and four samples on the pressure records (J = M = 4), and a window width of 1.42s (0.37-), resulting in 20 equations in 11 unknowns. There were no data available about the Eulerian current, so it is taken as zero. Window Width In the previous section, results on theoretically generated records . indicated the successful solutions were possible with quite narrow windows, gener- ally one tenth of the zero crossing period. When working with this field record, the optimization converged with a window width that narrow, but it resulted in wild fluc- tuations in the predicted propagation direction from window to window. Figure 6.9 is the same segment of the record as in figure 6.8, computed with the same parameters, 146 Predicted kinematics at the center of the array Figure 6.8: LFI predictions (J = 3) of kinematics at the center of the array at the water surface for a field record. 147 except that the window width used was much narrower, 7 = 0.1T7-.). There are large discontinuities in the horizontal velocity in the z direction (third plot), near both crests. Figure 6.10 shows the parameters of the potential function from this solution. 6, where 0 = tan7!(k,/kz), is the local direction of propagation of the wave. There is a sudden shift in direction of 80° at the first crest, and over 100° at the second crest. This would result in unrealistic predictions of huge accelerations. Figure 6.11 shows the parameters of the solution computed with 7) = 0.37. With the wider window, the propagation direction varies smoothly around the peak direction of the record, Q = 83°. Array Size The need for a fairly wide window for the field solution is likely due to the small size of the pressure gauge array. Compact size is one of the major advantages of the DWG-1 array, as it is a single unit that is easily deployed, but it makes it susceptible to difficulties with local analysis. The propagation direction within each window is determined by the spatial gradients of the dynamic pressure, as estimated by the measurements. In segments of the record where the gradients are small, the predicted direction of propagation can fluctuate a great deal with only small change in any of the measured values. These changes could be the result of measurement error, or be small fluctuations in the actual dynamic pressure. In either case, the large resulting change in predicted propagation direction can be a source of error. The linear dispersion relation predicts that the wavelength corresponding to the peak period of this record in this depth of water would be over 42 m. The DWG-1 array is only 1.6 mon a side, which is about 0.04 of the approximate wavelength of the peak period waves. As a result, near the crest of the wave, where the water surface is close to horizontal, all three gauges are near this point, and the local estimate of the gradient of the pressure is very sensitive to small fluctuations. This effect is mitigated by using a larger window in time, so that a part of the wave with a higher gradient is part of the window, allowing the propagation direction to be better defined, and smoothing out any measurement errors. As minimum window size has been determined to be about 1/10 of the zero crossing frequency, an array size of 148 Predicted kinematics at the center of the array tw, Figure 6.9: LFI predictions (J = 3) of kinematics at the center of the array at the water surface for a field record, using a narrow window. 149 Solution Parameters (0) 2 4 6 8 10 12 14 16 ndt Figure 6.10: Parameters of the LFI solution for a field record, using a narrow window (To = UEIIE 150 Solution Parameters nat Figure 6.11: Parameters of the LFI solution for a field record, using a wide window (To = Orsi): 151 about 1/10 of a zero crossing wavelength would probably result in more stable results with the smaller windows. 6.6 Discussion The LFI method can be applied to interpretation of records from arrays of pressure measurements. It is effective for records from primarily unimodal seas, as well as the bimodal seas that are generated near reflecting surfaces. The result of the analysis is a complete description of the kinematics of the waves throughout the water depth, in the vicinity of the array. The method is local in that a separate potential function is used to describe each small segment of the record, using information from only that segment. This approach is rational in that it seeks to satisfy the governing equations of gravity waves, including the full nonlinear free surface boundary conditions. Using a potential function representing a single progressive wave, the method has been shown to be effective on theoretically generated records of steady waves in both deep and shallow water. This formulation is appropriate in locations far from reflecting surfaces, where the sea state is expected to be unimodal. When the method is applied with a potential function representing two intersecting waves, it is effective ' in accurately capturing the kinematics of standing and short crested waves, as result from the reflection of steady waves from a reflecting surface, such as a sea wall. This two wave formulation is appropriate near such reflecting surfaces, where the sea state is likely to be bimodal. A number of parameters must be specified to establish a solution, including the order, window width, number of water surface nodes defined, and number of samples on the measured records. The examples in this chapter indicate that the criteria used to define these parameters are similar to those developed in the analysis of point pressure records and arrays of water surface measurements. Fairly low order solutions (J = 2 or 3) are effective in deep water, while shallow water requires higher order solutions (J = 5 or 6). Window widths of 1/10 the zero crossing period of the record are effective on theoretical records. With field records, however, when the array used is small (less 152 than 1/10 a typical wavelength), the predictions of the propagation direction can be very sensitive to subtle fluctuations in the measure record. In order to find a stable solution, the window width must be increased to about 1/3 of the zero crossing frequency. This includes more of the record in the window, allowing the propagation direction to be more clearly defined. The examples in this Chapter were all computed using MATLAB running under the Linux operating system on an Intel 90MHz Pentium processor based PC. The shallow water steady wave (Fig 6.2) took the longest to compute, about 64min., and 471x10° floating point operations (flops). The steady deep water wave (Fig. 6.3) took the least time, 12min. and 22x10° flops, the standing wave (Fig. 6.4) took 16min. and 131x10® flops, the short crested wave (Fig. 6.6) took 13min. and 133x10° flops, and the field record (Fig. 6.8) took 17min. and 46 x10® flops. 153 Chapter 7 Conclusions This dissertation examines the problem of determining the kinematics of three di- mensional irregular seas from a variety of wave measurements. The presented methods and results lead to the following conclusions: e Currently used global methods are inadequate for determining accurately the kinematics of irregular waves. — Linear superposition does not satisfy the nonlinear free surface boundary conditions, resulting in large errors in the predicted kinematics near the water surface. — Directional spectral methods generally disregard the phase information, making it impossible to determine the detailed kinematics of directional seas. — Global methods of high enough order to accurately solve the free surface boundary conditions in three dimensions are prohibitively computationally intensive. e The Local Fourier Method for Irregular waves (LFI) can accurately describe the kinematics of two and three dimensional irregular seas by satisfying the full free surface boundary conditions. 154 — The local nature of the LFI method allows it to satisfy the nonlinear free surface boundary conditions at low order by finding separate solutions in each of a sequence of small windows in time. — The LFI method provides an accurate two dimensional description of the kinematics of irregular seas measured by a single water surface gauge or subsurface pressure gauge. — The LFI method provides an accurate three dimensional description of the kinematics of irregular seas measured by arrays of water surface gauges or arrays of subsurface pressure gauges, including the bimodal seas resulting from reflection from a vertical surface, such as a sea-wall or breakwater. — When working with subsurface pressure records, the LFI method is limited in its capability to capture the high frequency components of the sea state that decay rapidly with depth. — In order to accurately capture the detail of the measured records, the LFI method requires very accurate data, with high sampling rates. — The LFI method can be adapted to virtually any arrangement of wave measuring instruments. — The LFI method required substantial, but not prohibitive, computational resources. 7.1 Future Work The results presented indicate the promise of the LFI method, and local methods in general. There is still much work to be done before the method could be considered ” code. usable as a “black box’ Much of this work revolves around determining the numerical details of the nonlin- ear optimization. While the presented results provide some guidelines to determining appropriate parameters, considerable judgment and experimentation with any given record is required. In the future, it may be possible to establish criteria for defining the following numerical solution parameters: 155 e Window width e Order of the potential function e Number of water surface nodes Number of samples of the measured record e Number of iterations allowed in the optimization Some of the criteria that might be used to help define these values are: e Non-dimensional depth e Ursell number e Local curvature of record In principle, the LFI method may be extended to additional arrangements of in- struments. While it is straightforward, in theory, to adapt the method to virtually any arrangement of instruments, each different arrangement is likely to present a new set of numerical solution difficulties. One of the sources of this is the use of non- _ linear optimization. Different optimization problems tend to be unique, and require substantial experience in order to determine the appropriate approach. The experience gained in applying the LFI method to a variety of arrangements of instruments helps demonstrate the need for appropriate data. Whether this par- ticular method of analysis is employed, or any other nonlinear technique, the need for accurate data is paramount. The design of a field data collection program inherently includes assumptions as to how the resulting data will be analyzed. Currently used data collection programs are often designed with the assumption that the data will be analyzed with the common methods of linear spectral analysis. As a result the data may not be appropriate for the nonlinear analysis needed to determine the detailed kinematics of a measured wave field. Some of the considerations that should be included in the design of a field mea- surement program are: 156 e Array size appropriate for the analysis method chosen and wave field expected. e Measurement or prediction of the local current field. e Sampling rate that captures the detail of the record. Reduction of sensor noise and increased instrument accuracy. As the methods for measuring waves and interpreting those measurements become more accurate, the understanding of the complex processes in the ocean and coasts that are molded by waves will increase as well. This understanding should lead to safer coastal and offshore structures, and well as reducing the impact that these structures have on the sensitive coastal environment. Bibliography Baldock, T. E. and C. Swan (1994). Numerical calculations of large transient water waves. Applied Ocean Research 16, 101-112. Baldock, T. E. and C. Swan (1996). Extreme waves in shallow and intermediate water. Coastal Engineering 27, 21-46. Barker, C. H. and R. J. Sobey (1996). Irregular wave kinematics from a pressure record. In Proceedings of the International Conference on Coastal Engineering. ASCE. Barker, C. H. and R. J. Sobey (1997). Nonlinear kinematics from directional wave measurements. In E. Mansard (Ed.), JAHR Seminar: Multidirectional waves and their interaction with structures., pp. 343-356. IAHR. Bath, M. (1974). Spectral analysis in geophysics. Amsterdam: Elsevier. Battjes, J. A. (1982). Effects of short-crestedness on wave loads on long structures. Applied Ocean Research 4, 165-172. Bishop, C. T. and M. A. Donelan (1987). Measuring waves with pressure transduc- ers. Coastal Engineering 11, 309-329. Bodge, K. R. and R. G. Dean (1984). Wave measurement with differential pressure gauges. In Proceedings of the International Conference on Coastal Engineering, pp. 755-769. ASCE. Bouws, E., H. Gunther, and C. L. Vincent (1985). Similarity of the wind wave spectrum in finite depth water. Journal of Geophysical Research 90, 975-986. Chaplin, J. and K. Anastasiou (1980). Some implications of recent advances in wave 158 theories. In Proceedings of the International Conference on Coastal Engineering, Volume 1, pp. 31-49. ASCE. Chaplin, J. R. (1980). Developments of a stream-function theory. Coastal Engi- neering 3, 179-205. Chappelear, J. E. (1961). Direct numerical calculation of wave properties. Journal of Geophysical Research 66(2), 501-508. Chen, Y. Y. (1988). The standing gravity waves. Technical report, The Institute of Harbour and Marine Technology in Taiwan. in Chinese. Daemrich, k.-F. and W.-D. Eggert (1980). Investigations on irregular waves in hydraulic models. In Proceedings of the International Conference on Coastal Engineering, Volume 1, pp. 186-203. ASCE. de Boor, C. (1978). A Practical Guide to Splines. New York: Springer-Verlag. Dean, R. G. (1965). Stream function representation of nonlinear ocean waves. Jour- nal of Geophysical Research 70, 4561-4572. Dean, R. G. (1977). Hybrid method of computing wave loading. In 9th Offshore Technology Conference. Dean, R. G. (1990). Stream function wave theory and applications. In J. B. Herbich (Ed.), Handbook of coastal and ocean engineering, Vol. 1, pp. 63-94. Houston, Texas: Gulf Publishing. Dean, R. G. and R. A. Dalrymple (1991). Water wave mechanics for engineers and scientists. World Scientific. Dold, J. and D. H. Peregrine (1984). Steep unsteady water waves: an efficient com- putational scheme. In Proceedings of the International Conference on Coastal Engineering, Volume 1, pp. 955-967. ASCE. Duncan, P. E. and K. R. Drake (1995). A note on the simulation and analysis of irregular non-linear waves. Applied Ocean Research 17, 1-8. Duncan, P. E. and Drake, K. R. 159 Fenton, J. D. (1985a). A fifth-order Stokes theory for steady waves. Journal of Waterway, Port, Coastal and Ocean Engineering 111(2), 216-234. Fenton, J. D. (1985b). Wave forces on vertical walls. In Proceedings of the Interna- tional Conference on Coastal Engineering, Volume 2, pp. 315-324. ICCE. Fenton, J. D. (1985c). Wave forces on vertical walls. Journal of Waterway, Port, Coastal and Ocean Engineering 111(4), 693-718. Fenton, J. D. (1986). Polynomial approximation and water waves. In Proceedings of the International Conference on Coastal Engineering, pp. 384-396. ASCE. Fenton, J. D. (1990). Nonlinear wave theories. In B. LeMéhauté and D. M. Hanes (Eds.), The Sea: Ocean Engineering Science, Vol. 9, pp. 3-25. Wiley. Fenton, J. D. and C. D. Christian (1989). Inferring wave properties from sub-surface pressure data. In Proceedings of the 9th Australasian Conference on Coastal and Ocean Engineering, Adelaide, pp. 380-384. Fenton, J. D. and W. D. McKee (1990). On calculating the lengths of water waves. Coastal Engineering 14, 499-513. Fenton, J. D. and M. M. Rienecker (1980). Accurate numerical solutions for non- linear waves. In Proceedings of the International Conference on Coastal Engi- neering, pp. 50-69. ASCE. Forristall, G. Z. (1985). Irregular wave kinematics from a kinematic boundary con- dition fit. Applied Ocean Research 7, 202-212. Forristall, G. Z. (1986). Kinematics in the crests of storm waves. In Proceedings of the International Conference on Coastal Engineering, pp. 208-222. Forristall, G. Z., E. G. Ward, V. J. Cardone, and L. E. Borgmann (1978). The directional spectra and kinematics of surface gravity waves in tropical storm Delia. Journal of Physical Oceanography 8, 888-909. Goda, Y. (1967). The fourth order approximation to the pressure of standing waves. Coastal Engineering in Japan 10, 1-11. 160 Grace, A. (1992). MATLAB Optimization Toolbox User’s Guide. The Mathworks, Inc. Grace, R. A. (1978). Surface wave heights from pressure records. Coastal Engineer- ing 2, 55-67. Gudmestad, O. T. and J. J. Connor (1986). Engineering approximations to non- linear deepwater waves. Applied Ocean Research &, 76-88. Hamm, L., P. A. Madsen, and D. H. Peregrine (1993). Wave transformation in the nearshore zone: a review. Coastal Engineering 21, 5-39. Hammack, J., D. McCallister, N. Scheffner, and H. Segur (1995). Two-dimensional periodic waves in shallow water part 2. asymmetric waves. Journal of Fluid Mechanics 285, 95-122. Hammack, J., N. Scheffner, and H. Segur (1989). Two-dimensional periodic waves in shallow water. Journal of Fluid Mechanics 209, 567-589. Herbers, T. H. and R. T. Guza (1989). Estimation of wave radiation stresses from slope array data. Journal of Geophysical Research 94, 2099-2104. Herbers, T. H. and R. T. Guza (1990). Estimation of directional wave spectra from multicomponent observations. Journal of Physical Oceanography 11, 1703-1723. Herbers, T. H. and R. T. Guza (1991). Wind-wave nonlinearity observed at the sea floor. part IJ: wave numbers and third-order statistics. Journal of Physical Oceanography 22, 489-503. Howell, G. L. (1992). A nearshore directional wave gage. In Proceedings of the In- ternational Conference on Coastal Engineering, Volume 1, pp. 295-307. ASCE. Hsu, J. R. (1990). Short-crested waves. In J. B. Herbich (Ed.), Handbook of Coastal and Ocean Engineering, Vol. 1, pp. 95-175. Houston, Texas: Gulf Publishing. Hsu, J. R. and Y. Y. Chen (1992). Frequency modulations between two intersecting waves in deep water. In L. Debnath (Ed.), Nonlinear Dispersive Wave Systems, pp. 299-328. World Scientific. 161 Hsu, J. R., Y. Tsuchiya, and R. Silvester (1979). Third-order approximation to short crested waves. Journal of Fluid Mechanics 90, 179-196. Ippen, A. T. and F. Raichlen (1957). Turbulence in civil engineering: measurements in free surface streams. Journal of the Hyraulic Division of ASCE 83(HY5), 1392.1-1392.17. Kishida, N. and R. J. Sobey (1988). Stokes theory for waves on a linear shear current. Journal of Engineering Mechanics 114, 1317-1334. Lambrakos, K. F. (1981). Extended velocity potential wave kinematics. Journal of Waterway, Port, Coastal and Ocean Engineering 107, 159-174. Lee, D.-Y. and H. Wang (1984). Measurement of surface waves from subsurface gage. In Proceedings of the International Conference on Coastal Engineering, pp. 271-286. ASCE. Lo, J.-M. and R. G. Dean (1986). Evaluation of a modified stretched linear wave theory. In Proceedings of the International Conference on Coastal Engineering, pp. 522-536. ASCE. Longuet-Higgins, M. S. (1962). Resonant interactions between two trains of gravity waves. Journal of Fluid Mechanics 12, 321-332. Longuet-Higgins, M. S. (1974). On the mass, momentum, energy, and circulation of a solitary wave. Proceedings of the Royal Society of London, Series A 337, IS}. Longuet-Higgins, M. S. (1975). Integral properties of periodic gravity waves of finite amplitude. Proceedings of the Royal Society of London, Series A 342, 157-174. Longuet-Higgins, M. S. (1987). The propagation of short surface waves on longer gravity waves. Journal of Fluid Mechanics 177, 293-306. Longuet-Higgins, M. S., D. E. Cartwright, and N. D. Smith (1963). Observations of the directional spectrum of sea waves. In Ocean Wave Spectra, pp. 111-136. Prentice-Hall. 162 Longuet-Higgins, M.S. and J. D. Fenton (1974). On the mass, momentum, energy, and circulation of a solitary wave II. Proceedings of the Royal Society of London, Series A 340, 471-493. Longuet-Higgins, M.S. and O. M. Phillips (1962). Phase velocity effects in tertiary wave interactions. Journal of Fluid Mechanics 12, 333-336. Mansard, E. (Ed.) (1997). JAHR Seminar: Multidirectional Waves and their In- teraction with Structures, Canadian Hydraulics Center, Ottawa, On K1A OR6. The National Research Council of Canada. Mei, C. C. (1990). Basic gravity wave theory. In J. B. Herbich (Ed.), Handbook of coastal and ocean engineering, Vol. 1, pp. 1-62. Houston, Texas: Gulf Publish- ing. Meiron, D. I., P. G. Saffman, and H. C. Yuen (1982). Calculation of steady three- dimensional deep-water waves. Journal of Fluid Mechanics 124, 109-121. Newland, D. E. (1993). Random Vibrations, Spectral and Wavelet Analysis. Long- man Scientific and Technical. Nielsen, P. (1986). Local approximations: A new way of dealing with irregular waves. In Proceedings of the International Conference on Coastal Engineering, pp. 633-646. ICCE. Nielsen, P. (1989). Analysis of natural waves by local approximations. Journal of Waterway, Port, Coastal and Ocean Engineering 115, 384-396. Ohyama, T., D.-S. Jeng, and J. R. Hsu (1995a). Fourth-order theory for multiple- wave interaction. Coastal Engineering 25, 43-63. Ohyama, T., D.-S. Jeng, and J. R. Hsu (1995b). Fourth-order theory for multiple- wave interaction: Erratum. Coastal Engineering 26, 297-298. Ohyama, T., W. Kioka, and A. Tada (1995). Applicability of numerical models to nonlinear dispersive waves. Coastal Engineering 24, 297-313. Ohyama, T. and K. Nadaoka (1994). Transformation of a nonlinear wave train 163 passing over a submerged shelf without breaking. Coastal Engineering 24, 1- Wie O’Reilly, W. C., T. H. Herbers, R. J. Seymour, and R. T. Guza (1996). A com- parison of directional buoy and fixed platform measurements of Pacific swell. Journal of Atmospheric and Oceanic Technology 13, 221-238. Phiillips, O. M. (1960). On the dynamics of unsteady gravity waves of finite am- plitude. Journal of Fluid Mechanics 9, 193-217. Prislin, I. and J. Zhang (1997). Prediction of wave kinematics under the wave crest in short-crested irregular gravity waves. In E. Mansard (Ed.), JAHR Seminar: Multidirectional Waves and their Interaction with Structures, pp. 357-370. Prislin, I., J. Zhang, and R. J. Seymore (1997). Deterministic decomposition of deep water short-crested irregular gravity waves. Journal of Geophysical Re- search 102(C6), 12677-12688. Raichlen, F., J. D. Ramsden, and J. R. Walker (1990). Bottom pressures due to long waves: laboratory and field measurements. In Proceedings of the International Conference on Coastal Engineering, pp. 1144-1159. ASCE. Rienecker, M. and J. Fenton (1981). A Fourier approximation for steady water waves. Journal of Fluid Mechanics 135, 301-321. Roach, P. J. (1976). Computational Fluid Dynamics. Albuquerque: Hermosa Pub- lishers. Roberts, A. J. (1983). Highly nonlinear short-crested water waves. Journal of Fluid Mechanics 135, 301-321. Roberts, A. J. and L. W. Schwartz (1983). Calculation of nonlinear short-crested gravity waves. Physics of Fluids 26, 2388-2392. Rodenbush, G. (1986). Random directional wave forces on template offshore plat- forms. In Offshore Technology Conference, pp. 147-156. Rottman, J. W. (1982). Steep standing waves at a fluid interface. Journal of Fluid Mechanics 124, 283-306. 164 Salvadori, M. G. and M. L. Baron (1961). Numerical Methods in Engineering. Prentice-Hall. Sheffner, N. (1986). Biperiodic waves in shallow water. In Proceedings of the Inter- national Conference on Coastal Engineering, pp. 724-736. ASCE. Skjelbreia, L. and J. Hendrickson (1962). Fifth order gravity wave theory. In Pro- ceedings of the 7th conference on coastal engineering, pp. 293-306. ICCE. Sobey, R. J. (1989). Variations on Fourier wave theory. International Journal for Numerical Methods in Fluids 9, 1453-1467. Sobey, R. J. (1990). Wave theory predictions of crest kinematics. In A. Torum and O. T. Gudmestad (Eds.), Water Wave Kinematics, pp. 215-231. Netherlands: Kluwer Academic Publishers. Sobey, R. J. (1992). A local Fourier approximation for irregular wave kinematics. Applied Ocean Research 14, 93-105. Steele, K. E., C.-C. Teng, and D. W. C. Wang (1992). Wave direction measurements using pitch-roll buoys. Ocean Engineering 19, 349-375. Stokes, G. G. (1880). Mathematical and Physical Papers: Vol I, Chapter On the theory of oscillatory waves, pp. 197-229. Cambridge University Press. Su, M.-Y. (1982). Three dimensional deep-water waves. part 1. experimental mea- surement of skew and symmetric wave patterns. Journal of Fluid Mechan- ics 124, 73-108. Sultan, N. J. (1992). Irregular wave-induced velocitites in shallow water. Technical Report CERC-92-9, Coastal Engineering Research Center, US Army Engineer- ing Waterways Experiment Station, Vicksburg, Mississippi USA. Sultan, N. J. and S. A. Hughes (1993). Irregular wave-induced velocities in shallow water. Journal of Waterway, Port, Coastal and Ocean Engineering 119, 429- 447. Swan, C. (1990a). A viscous modification to the oscillatory motion beneath a series of progressive gravity waves. In A. Torum and O. T. Gudmestad (Eds.), Water 165 Wave Kinematics, pp. 313-329. Netherlands: Kluwer Academic Publishers. Swan, C. (1990b). Wave kinematics within the crest to trough region. In Environ- mental Forces on Offshore Structures and their Prediction, Vol. 26, pp. 313-329. Netherlands: Society for Underwater Technology. Tsai, C.-P. and D.-S. Jeng (1994). Numerical Fourier solutions of standing waves in finite water depth. Applied Ocean Research 16, 185-193. Tsuchiya, Y. and T. Yasuda (1984). A dynamical expression of waves in shallow water. In Proceedings of the International Conference on Coastal Engineering, pp. 435-451. ASCE. Vis, F. C. (1980). Orbital velocities in irregular waves. In Proceedings of the Inter- national Conference on Coastal Engineering, pp. 173-185. ASCE. Wheeler, J. D. (1969). Method for calculating forces produced by irregular waves. In Offshore Technology Conference, pp. 71-82. Zelt, J. A. and O. T. Gudmestad (1995). Fluid accelerations under irregular waves. Applied Ocean Research 17, 43-54. 166 Appendix A Intersecting wave theory A.1 Introduction The theory of multiple wave interactions presented by Ohyama, Jeng, and Hsu (1995a) has been used for two purposes in this work. The most important function was to suggest a form for a general three dimensional potential function representing interacting waves. The form chosen is presented in Chapter 4. It also provided a theoretical method for generating a complete, consistent set of records to use for the development and testing of the LFI method. This appendix presents a review and detailed critique of the intersecting wave theory. A.2 Theoretical Background Ohyama et al.’s method is an extension of well accepted and verified methods. It is essentially a Stokes type expansion in the wave steepness. As such, it is ex- pected have greatest applicability in deep water. The flow is taken to be irrotational and incompressible, and thus the governing equations are the same as those used in chapter 4 and for most finite depth irrotational wave theory: e Mass conservation (Laplace equation) (Eq. 4.2) e Bottom boundary condition (Eq. 4.3) 167 e Dynamic free surface boundary condition (Eq. 4.4) e Modified kinematic free surface boundary condition (Eq. 4.6) There is no accommodation for an Eulerian current, although it might easily be included in future work. The non-dimensional potential function, water surface, and frequencies are expanded in a power series in the small parameter, (€ = ka), where k is a typical wave number and a is a typical amplitude of the first order component of the water surface. b=) + 2d + Hd) + AG + O(€) (A.1) p= GB Reg? + eg” + Gg £0) (A.2) w= eu”) + ew?) + eu” + Au) + O(e?) (A.3) Where the dé”, n'” are the potential function and water surface at order n, and ual”) is the frequency of the zth wave at order n. The steepness of all the intersecting waves is taken to be of the same order. Modern Stokes theories have found it necessary to use an expansion parameter based on the physical wave height, rather than the amplitude of the first order component (Fenton 1990). This is not feasible with intersecting waves, as there is no clearly defined wave height. While not ideal, Stokes methods based on this first order expansion parameter can be effective (Skjelbreia and Hendrickson 1962), as long as they are algebraically correct and not pushed to near limit waves or shallow water. Care must be taken when computing a given order solution when multiple parts of the solution are all expanded in a perturbation series. In this case, the potential function, water surface, and frequencies are expanded separately, but the potential function and water surface are dependent on the full order frequency. Ideally, a given order component will depend only upon parameters of the same order. In this case, however, the full order frequency must first be computed, and then used as the frequency at all orders. If the order of the frequency used was matched to the order of 7 or ¢ being computed, the resulting solutions would be out of phase, with, 168 for example, 6") having a different frequency to 6). It could be argued that this approach would be accurate to second order, but the components would be shifted out of phase as time progressed. It is more accurate to compute the results using the full order frequencies at all orders of ¢ and 7. The first order solution for N intersecting waves is the familiar superposition of linear waves: N (@) 2 a; cosh k(h + 2) a (if ree 4 ro) De i a eee sin (k,z + kyy — wt + a) (A.4) 7) = yy a; cos (ky + kyy — wt + a) (A.5) i=1 w'') = ./gk; tanh kh; (A.6) The higher order solutions are computed by applying the Laplace equation, the bottom boundary condition, and the free surface boundary conditions, expanded in Taylor series about the mean water level. The algebra and the resulting solutions are long and extremely tedious, are presented in Ohyama et al. (1995a), and they will not be repeated here. It should be noted that the frequency was only solved to third order, as the fourth order frequency would require a solution to the fifth order potential function, which was not presented. A.3 Analytical Verification Ohyama et al. provided a number of comparisons with other wave theories. When the method is applied for a single wave component, it represents a steady wave, and as such should be expected to provide the same solution as conventional Stokes theory. The authors presented a comparison with the Fenton (1985a) solution. The Fenton solution is based on an expansion parameter based on the wave height, rather than the first order amplitude: € = kH/2. By solving for € in terms of the «, The two solutions were found to coincide exactly. 169 When applied for two wave components of the same wavelength, same amplitude, and opposing directions, the solution is a standing wave. The resulting solution coincides almost exactly with a standing wave solution given by Chen (1988), except for one coefficient. The authors suggest that the slight discrepancy is the result of a typographical error in Chen’s paper. Ohyama et al. also indicated that they compared their solution with those for short-crested seas (Hsu 1990; Hsu and Chen 1992; Fenton 1985c) and once again found essentially perfect agreement. A.4 Numerical Verification The above theoretical verifications indicate that the method is correct when re- duced to two particular simplifications. It remains to verify that the method is accurate when applied to a more complex situation, with a number of intersecting waves of different heights, directions and wavelengths. The actual computation of the general solution is quite complicated and contains a very large number of coefficients that must be computed. Despite this complexity, once the code is written and debugged, it is a trivial matter for modern computers - to calculate the full solution. It is not, however, a trivial matter to write the code without errors or, indeed, to find the likely typographical errors in the published paper. A number of these typographical errors have been presented in the published erratum (Ohyama et al. 1995b). An additional error has been found and confirmed by personal communication with the authors. pp. 55 Eq. 47 in Ohyama, Jeng, and Hsu (1995a) should read: ¢ = {£?Bo20 = Baro} cos 2kx + €4Bs49 cos 4kr + {2 2s "Baus } cos kx cos wt + ror + <“Bixy} cos 2kz cos 2wt + €°Ba13 cos kx cos 3wt SES hee cos 3kz coswt + €? Boga cos 3kz cos 3wt + a" Bana cos 2kz cos 4wt +€4B442 cos 4k cos Qwt + €4B444 cos 4k cos 4wt + Ofe*} As a result of the complexity of the computational code, it is advisable to use numer- 170 ical methods to verify the validity of the solution and resulting computational code. The following results were produced by a Fortran code derived from the code written and generously provided by Takumi Ohyama that was used to compute the water surface plots given in (Ohyama et al. 1995a). The code was expanded to compute the full kinematics. A.4.1 Richardson Extrapolation Fenton (1985a) presented a variation of Richardson extrapolation that can be used for numerical verification of wave theory code. Richardson extrapolation (also known as extrapolation to the limit) is a method traditionally used to increase the accuracy of finite difference computations (Salvadori and Baron 1961; Roach 1976). In Fenton’s adaptation, the method is used to simultaneously check all the terms in the solution, and to provide an estimate of the order of the accuracy of the result. The form of Ohyama et al.’s solution solves the field equation (Laplace) and the bottom boundary condition exactly. It is expected to satisfy the free surface boundary conditions to the order it is applied. The following method calculates the order of accuracy of the free surface boundary conditions. Classical Richardson extrapolation is used to evaluate the error of a solution at a given point in time and space. By taking advantage of the periodicity of wave motion, Fenton’s variation can be used to evaluate the solution throughout a complete wave. Single Wave The solution for a single wave is periodic and symmetric about the crest, and thus the error in each of the boundary conditions can be expanded in a single sided Fourier series: ges ge” (A.7) 171 [ Fourier Component, (| 0] 1)2)3]4]5] 6] [_ Order of KFSBC error [2.0 [2.0/2.0] 2.02.01 2.0] 2.0] [_Order of DFSBC error [3.0] 2.2 [2.0 [2.02.0] 2.0] 200 | Table A.1: Order of errors for a single wave computed to order 1 [Fourier Component, (K)] 0] 1]2]3]4][5] 6 [_ Order of KFSBC error | 3.0/3.0 30 [30/30] 3.0/3.0] [Order of DFSBC enor [2.9/3.0] 3.0/3.0] 3.0/3.0] 3.0] Table A.2: Order of errors for a single wave computed to order 2 The magnitude of this error, and thus each of the coefficients can be expanded in a first order Taylor series. Cy(Ex) = Bub. + O (&) (A.8) These coefficients are a function of the small parameter, € and thus are a function of €” at nth order, &, =e". C.(€) = Bre +O (e*) (A.9) When the discreet transforms of the errors are computed for two different values of €, an approximation for the n, can be computed, giving an approximation to the order of the errors of each harmonic in each boundary condition. iy, = GaGa + O(€1, €2) (A.10) Table A.1 presents the results for a single wave, computed to order | in deep water (kh = 100, & = 0.01, €2 = 0.02). These results show that the error is of at least second order in both of the boundary conditions, indicating that the computation is accurate to first order. Tables A.2 through A.4 present the results for the same wave, computed to sec- ond through fourth order. These results indicate that the error in both boundary 172 [Fourier Component, ()] 0] 1]2]3]4]5] 6 [_ Order of KFSBC error [4.0 [4.0] 4.0 [4.0 [4.0 [4.0/4.0] [Order of DFSBC error | 4.0/3.9] 4.0] 41 [4.0] 40] 40] Table A.3: Order of errors for a single wave computed to order 3 [Fourier Component, (K)] 0] 1]2[3]4]5] 6] [Order of KFSBC eor [5.0 [5.0 [5.0] 5.0] 5.0 |5.0 | 5.0 | [_ Order of DFSBC error [48 [5.0/5.0] 5.0 [5.0 [5.050] Table A.4: Order of errors for a single wave computed to order 4 conditions is of at least one order higher than the order used to compute the results. These results serve to confirm the efficacy of the Richardson extrapolation method, as well as the accuracy of the wave theory up to fourth order. Standing Wave Ohyama et al.’s method can also be used to compute a standing wave by defining two waves of the same wavelength and amplitude, and opposing directions. In this case, the wave is not steady, and the errors in the free surface boundary conditions must be considered for a full period in time and full wavelength in space. The ex- trapolation to the limit method can be extended to this case by expanding the free surface errors in a two dimensional Fourier series (Newland 1993). =) Ss Cp (jee e/E+H/7) (A.11) k=0 I=0 An estimate of the order of the errors can then be computed in a similar manner as Eq. A.10. The results of this computation for a standing wave in deep water (kh = 100, €; = 0.01, €2 = 0.02) are given in Tables A.5 and A.6. There are no values given for the zeroth order in time of the kinematic boundary condition because C(€),,o for both values of € are zero to the precision of the calculation, so the values computed are spurious. These results indicate that Ohyama et al.’s method is also accurate to 173 [Fourier Component, [0] 1]2[3]4)]5] 6] (ea [1 [50 [50 [50 [50 [5050/50 en dd Ne [Es 5050/8050 [5.05.0 | 5 | a [ES 49 [51 [as [50 [51 50) 5.0 [Ee [5050 [51 [5.0 [5.0 [50] 5.3 | Table A.5: Order of KFSBC errors for a standing wave computed to order 4 [Fourier Component. (0 ]1]2[/3]4]5] 6] ee EES a4 [50/4150] 49]50 [50 | [61 [52 [6.5 [5.0] 5.8 [5.5 [5.3 | [5 [495.0 [5.2] 5152] 8.049 | [E853 [5.0] 5.95.0] 47 [8.0 | 53 | Table A.6: Order of DFSBC errors for a standing wave computed to order 4 approximately fifth order for a standing wave. A.4.2 More Complex Seas Fenton’s method relies on the periodicity of the solution to be tested. If the solution is not periodic in space and time, then the errors in the boundary conditions will not be periodic, and the Fourier expansion cannot be strictly applied. Fourier coefficients could still be computed by assuming periodicity, but the discontinuity generated by the assumption of periodicity would generate spurious results. These could be minimized by taking a large number of points over a large range in time and space, but this would become very computationally intensive. The accuracy of any approximate solution can be measured by examining the 174 oder (a) 342 x 107 Table A.7: RMS errors in free surface boundary conditions for a short-crested sea errors in the governing equations resulting from the solution. Fenton’s method is a particularly powerful method for examining these errors, in that it considers an entire wavelength and period at once and gives an estimation of the order of accuracy of the solution. A simpler approach can also be useful and be universally applied. In the case of Ohyama et al.’s method, the field equation and bottom boundary condition are exactly satisfied at all orders of approximation. If the solution is accurate, the errors in the free surface boundary conditions will decrease with each increase in the order of approximation. While this method does not guarantee that the solution is correct, most errors in the algebra or the computational code will result in a decrease in the accuracy of the result at the order in which the error occurs. Table A.7 contains the errors in the free surface boundary conditions at orders one through four for a short crested sea. This sea is created by two waves of the same height and wavelength (« = 0.1, kh = 100) intersecting at an angle of 30 degrees. The values given are the root mean square of the errors computed by averaging over a grid of points one wavelength in both horizontal directions, and over one period. In this case the errors in the free surface boundary conditions decrease with increasing order of solution. It should be noted, however, that the errors do not decrease as much as e”. This might be partially due to the fact that the frequency is expressed to only second order. Including the fourth order frequency would probably reduce the error at fourth order. This behavior warrants closer examination, particularly if the method were to be expanded to higher order. 175 A.5 Depths of Validity As the method is very similar to Stokes theory, it is expected to be valid in a similar range of depths. Stokes theory is valid if both ¢ and the parameter, ¢/(kh)*, similar to the Ursell number, are small (Fenton 1985a). This indicates that the method has optimal applicability in deep water, and should not be applied in shallow water. Vie re ) , —_______— REPORT DOCUMENTATION PAGE ; Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining l | the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions | | for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the | | Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. | 1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED | September 1998 Final report 14. TITLE AND SUBTITLE 5. FUNDING NUMBERS | Directional Irregular Wave Kinematics 6. AUTHOR(S) Christopher H. Barker, Rodney J. Sobey 8. PERFORMING ORGANIZATION REPORT NUMBER Technical Report CHL-98-24 | 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) University of California at Berkeley Berkeley, CA 94720 10. SPONSORING/MONITORING AGENCY REPORT NUMBER |9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) U.S. Army Corps of Engineers, Washington, DC 20314-1000 U.S. Army Engineer Waterways Experiment Station Road, Vicksburg, MS__ 39180-6199 12b. DISTRIBUTION CODE 7 12a. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; distribution is unlimited. }13. ABSTRACT (Maximum 200 words) Coastal and ocean processes are heavily influenced by the kinematics of waves. In order to understand these processes, researchers place a variety of instruments in the sea in an attempt to measure the waves. These instruments all measure a small set of physical quantities at a small number of locations. The balance of the kinematics must be predicted through analysis of the measured records. Most of the currently used methods of analysis rely on the superposition of linear waves to recreate complex seas. These methods are compromised by linearizing approximations to the free surface boundary conditions. Fidelity in the interpretation of wave measurements is enhanced by insisting that the analysis satisfies the full nonlinear free surface boundary conditions. The Local Fourier method for irregular wave kinematics (LFI) is introduced and expanded to include the interpretation of records from arrays of instruments. It is a local method, in that a separate solution is sought that fits the measured record(s) in a small local window in time, rather than attempting to find a single solution for a large segment of the record. Each window solution satisfies the full set of governing equations for gravity waves, including the nonlinear free surface boundary conditions. The solution in each window is a potential function whose form is based upon a Stokes style expansion for intersecting waves. The parameters of the potential function are found by a nonlinear optimization that seeks the solution that matches the measured record and satisfies the full free surface boundary conditions. (Continued) . SUBJECT TERMS Crest kinematics Wave kinematics Local fourier method Wave measurements Nonlinear waves Waves Wave gauges Standard Form 298 (Rev. 2-89) Prescribed by ANSI Std. Z39-18 298-102 13. (Concluded). The LFI method was introduced by Sobey (1992) as a two dimensional method for interpreting a single water surface measurement of irregular waves. In this dissertation, the two dimensional method is expanded to include the analysis of subsurface pressure records. The method is then extended to the prediction of directional irregular wave kinematics from the measurements of arrays of instruments, including wave staffs and subsurface pressure gauges. In this extended analysis, the LFI method includes nonlinear wave-wave interactions, in addition to the self-wave interactions of the two dimensional form. Results of the analysis on theoretical, laboratory and field records are provided. Destroy this report when no longer needed. Do not returm it to the originator. SEP 13'93 MS DEPARTMENT OF THE ARMY WATERWAYS EXPERIMENT STATION, CORPS OF ENGINEERS 3909 HALLS FERRY ROAD SPECIAL VICKSBURG, MISSISSIPPI 39180-6199 FOURTH CLASS BOOKS/F ILM Official Business 256/L25/ 1 DATA/DOCUMENT LIBRARY, WHOI MCLEAN LAB, MS #8 360 WOOD HOLE ROAD WOODS HOLE MA 02543-1539 Be nirh. apr MMETER 576915