Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

Nonlinear Simulation of Transient Free Surface Flows R. Cointe Institut Frangais du Petrole Rueil-Malma~son, France 1 Abstract The application of the Mixed Eulerian-Lagrangian method to the simulation of transient free surface flows in the vicin- ity of a free surface piercing structure is considered. A par- ticular attention is given to the validation of the numerical procedure. Several applications are studied. Comparisons between the results of the numerical scheme and those of approxi- mate theories and/or experiments are shown. They demon- strate the accuracy and versatility of the simulation that can be used as a "standard" to check the applicability of approximate theories. The main limitation of the method is that it cannot account for viscous effects, in particular in the vicinity of the free surface. Approximate ways to simulate dissipative phenomena associated to breaking would be most useful. 2 Introduction The direct numerical simulation of unsteady two-dimen- sional potential free surface flows using a Mixed Eulerian- Lagrangian method (MEL) has received considerable at- tention since the pioneering work of Longuet-Higgins and Cokelet t224. If many codes exist now that use this method, their suitability for the study of nonlinear fluid-structure interaction problems has only been demonstrated for some particular applications (e.g., water impact t15], simulation of breaking in a tank t12], forced heaving of a cylinder t18], etc...~. Compared to the initial application of the MEL to the study of steep periodic waves (e.g., t22~), several new dif- ficulties appear when a structure is present, especially if it pierces the free surface. The first one is related to the proper description of the flow in the vicinity of intersection points between the body and the free surface. A second difficulty is to be able to control the incident wave train that interacts with the body. The main difficulty, however, is probably that the accu- racy of the simulation is difficult to establish for transient flows in the vicinity of free surface piercing bodies, because of a lack of reference cases. Checking such requirements as conservation of fluid or energy might not be sufficient. Surprisingly enough, the linear solution is not always com- puted and quantitative validations of the numerical proce- dure are almost inexistant, even for weal;ly nonlinear flows. These considerations have partly motivated the present work which associates numerical and analytical studies. A code based on the MEL Sindbad has been de- velopped that has for purpose the simulation of a two- dimensional wave tank using potential flow theory. A par- ticular attention is given to the proper validation of the numerical scheme by comparison of its results to those of experiments and of asymptotic studies. To make this com- parison easier, an option of the code allows the linear prob- lem to be solved. Such an approach appeared necessary in order to gain confidence in the code and eventually extend its range of applicability. A direct simulation of experiments that can be carried out in a wave tank is performed. For this purpose, waves are generated in a rectangular tank by use of a piston- type wavemaker. These waves can then interact with sub- merged or free surface piercing cylinders in forced or free motion. Reflection from the wall opposite to the wave- . . . . . maker is avoided by the use of a damping zone see fig- ure 1 The results of the classical first- and second-order dif- fraction radiation theories can be recovered using the sim- ulation. A good agreement with experimental results is obtained in cases where these classical theories fail. The main limitation of the method is that viscous and turbulent effects cannot be accounted for. This problem is crucial when viscous effects occur in the near vicinity of the free surface, in particular during breaking. This work has been reported elsewhere while in progress (hi, Pi, [8], 9. If we focus here on the results of the simulation, a more detailed description of the numerics and of some of the approximate theories referred to here can be found in t74. 239

Y ~Dnu~nerico/u~or~a Bank l g-1 o EM 3 Numerical Scheme 3.1 Outline of the Method Since there exists now quite a large number of codes based on the MEL, the numerical method used will only be briefly outlined. The attention will focus in the next sections on the main difficulties that have been encountered and on the methods used to overcome them. The main idea of the numerical procedure is to choose markers initially at the free surface and to follow them in their motion. We use a coordinate system (x, y). The x-axis coincides with the reference position of the free surface and the y-axis is upward vertical - see figure 1 for geometric definitions. The fluid is assumed to be incompressible and the flow irrotational so that the velocity field v is given by: v = V¢, (1) with: A¢' = 0. (2) The computation is performed in a bounded domain. Along rigid boundaries (x ~ En~t)~, the normal velocity in the fluid is equal to the normal velocity of the boundary. We have, therefore, a Neumann boundary condition. Along the free surface, we use Bernoulli equation and the fact that the free surface is a material surface. The corresponding equations are written for a marker x on the free surface (x ~ I`~(t)) and the associated value of the potential, ¢(x). This yields): Do _y _ ~ ~ _ <' ¢2 + ~ ¢2 _ p + cats (3) Do _ ~ ~ ~ 1 ~ ~ `4) where D is used to indicate a material derivative, s and n are vectors tangent and normal to the free surface, respec- tively, and ~ is an arbitrary constant. - ~ For the sake of simplicity, we use units such that the acceleration of gravity, 9, the specific mass of water, p, and the depth of the tank, h, are equal to 1. Piston-type Tested body wave maker Figure 1: Sindbad geometric definitions 240 h-l "Beach" '(absorbing, zone) This constant specifies the tangential motion of the markers: ~ = 1 identifies markers and particles while ~ = 0 yields a zero tangential motion of the markers. This last choice allows a current to be simulated in the tank. For the applications discussed here, we will however always take = 1. We will assume that the pressure is constant along the free surface. It can therefore be included in the function of time city. With an appropriate choice of the velocity potential, this function can be taken equal to zero. At a given instant t, if o and A', are known along the free surface Fatty, then the right-hand sides of (3~-~4) can be evaluated. The fact that ¢> is harmonic in the fluid domain allows the value of ¢'n along I`~(t) to be computed from the values of o along I,~(t) and of ¢>72 along I`ntt). The kinematical constraint /~¢ = 0, associated with the boundary condition on rtn' permits therefore to express the free surface boundary conditions (3~-~4) as an evolution equation for (¢>, x). This evolution equation can be solved numerically using standard time-stepping procedures, such as a fourth-order Runge-Kutta algorithm. The main numerical difficulty is to be able, at each time-step, to solve for the harmonic function ~ knowing ¢, along Lotte and ¢)n along En~t). We use the integral equation: ~ ~ ~ ~ ~r`~+rn ~ ~ ~ ~ Q = .,lr +r An (Q ~ G(P' Q ~ dog, (5) where P is a point on the boundary, G is the Green func- tion, 8~) the angle between two tangents of the boundary at P (equal to ~ for a smooth curve) and s a curvilinear abscissa along I`. Equation (5) is discretized using a stan- dard collocation method. The boundary of the domain is approximated by segments and ~ and ¢>r' are assumed to vary linearly along each segment. This allows an analyti- cal integration of the Green function, its normal derivative and their products by the curvilinear abscissa so that the calculation of the matrix elements is rather simple (and vectorizes well).

3.2 Numerical Treatment at the Intersections At each time-step, we have to solve a Neumann [along [nit)] - Dirichlet talon" l`~(t)] boundary value problem. It is well known that the solution of such a "mixed" prob- lem is singular at the intersection points, i.e., at points belonging both to I,n~t) and Pelt) e.g., [16;. Let us consider, as an example, the wavemaker problem. For a piston type wavemaker and an horizontal free surface (90° intersection angle), the complex potential ~ = ~ + i ~ solution of such a problem behaves like z log z near an intersection located at z 0. This gives the expected behavior of the problem discretized in time. Discretizing in time is similar to performing a small time expansion, i.e., to write: ¢(x, t) = jO(X) + t ¢~(x) + · · · (6) Not surprisingly, performing such an expansion also leads to a z log z behavior for ¢~ (e.g., t27i, t21~. Does this, however, necessarily imply that this singular behavior has to be expected for the solution of the transient problem? The answer is no, just because the small time expansion is not regular near the intersection point. It is, therefore, improper for a local analysis. A regular expansion can be found in the weakly nor~lin- ear regime, i.e. here, for a small acceleration of the wave- maker (relative to that of gravity) see A. It appears that in this case the first approximation (in an asymptotic sense) is provided by the classical linear solution. The boundary condition for this solution is not a Dirichlet boundary con- dition; it is given by: ¢~+¢y = 0. (7) If regularity in time is assumed, the local behavior of the complex potential for this solution can be shown to be in z2 logz. Not surprisingly, this singularity is the same as that appearing for the harmonic problem (for which = ~ exptiwt)) which was studied by I<ravtchenko t20] as early as 1954. It is much weaker than that of the prob- lem discretized in time. The singularity can also be studied in a similar fashion for an arbitrary angle of intersection 8. Assuming that the complex potential is bounded near the intersection, the leading behavior of the complex potential can be shown to be in z7r/8 or z7r/8 log z if 7r/8 is an integers. A consequence is that, for an intersection angle smaller than or, ~ is continuous along the boundary ~ while Us and In are piecewise continuous but experience a finite jump at the intersections. The numerical treatment at the intersection has been devised to accomodate such a behavior, referred to as weak- ly singular. More details on the numerical implementation can be found in t5] A. These results only apply in the weakly nonlinear regime, i.e., when the classical linear solution yields a first reason- able approximation to the problem. Obviously, it would be interesting to study other regimes. The impulsive regime (very large acceleration of the wavemal;er) has been studied by several authors. As e:;- 2The same singularity occurs at corners of rigid boundaries. pected, the local behavior of the classical linear solution cannot yield any information for the nonlinear problem in this case (e.g., t284~. It seems that Peregrine's solution t27] corresponds to an outer solution to which an inner solution should be matched. Of course, only this inner solution is relevant to the the study of the local behavior. If work has been carried out to find an inner solution~e.g., t21] and [8~), it leads to serious difficulties and, to our knowl- edge, no definitive answer has been provided. It is only for the related but different water entry problem that there exists some results concerning the local behavior of the nonlinear solution (e.g., try. 3.3 Wave Absorption In order to make proper comparisons with experiments and classical theories, it is often necessary to compute a steady state response. When the linear solution is computed, one prescribes an incident wave and write a "radiation" condition that transmitted and radiated waves have to satisfy. Writing a proper radiation condition for the second-order problem has long been a matter of controversy. For the fully non- linear problem, such an approach does not seem possible. In the absence of any mathematically satisfying answer, we have chosen a pragmatic solution similar to that used for an experiment in a tank. This approach does not involve any hypothesis concerning the steepness of the outgoing waves. Waves are generated by a wavemaker, either piston- type or flap-type. A "beach" is used for the absorption of the waves that are generated in the tank see figure 1. It is in fact a damping zone, similar to that used in A. In this zone, the free surface boundary conditions are modified by adding a damping term. We, therefore, write: D! (2 () As + 2 '¢)nv(~x) (~¢~ ~ (~8) DDX = ~ Us s + On nutter`) (<xate) (~9) where the subscript e corresponds to the reference config- uration for the fluid. The function Waxy is homogenous to a frequency. The principle of this damping zone is to absorb the incident wave energy before it can reach the wall. It may be intuited that if the absorption is too weak, part of the incident wave energy will reach the wall and be reflected. Inversely, if the absorption is too strong, part of the energy will be reflected by the damping zone itself. In practice, the damping coefficient Waxy is equal to zero except in the damping zone (x > no)' is chosen continuous and continuously differentiable, and is "tuned" to a char- acteristic wave frequency ~ and to a characteristic wave number k: waxy = <~w t2 x _ pJ2, x > To = 13. (10) If the proper scales have been chosen, values of or and ~ of order 1 should be appropriate for the absorption of a wave train of wave frequency co and wave number k.

3.4 Validation As indicated previously, the proper validation of the nu- merical simulation has been one of our main objectives. Solving the fully nonlinear problem can only be interesting if the accuracy of the numerical scheme is well established. Several strategies exist for this validation. Consistency checks should of course be performed. For that purpose, the volume and the rate of change of en- ergy are computed. When no damping zone is present, this last quantity is compared to the power input from exterior loads. The pressure exerted by the fluid on rigid bound- aries is computed by two different methods. Both of them use Bernoulli's equation, but they diner by the evaluation of the ¢' term. In one case, this term is evaluated by fi- nite differences in time. In the other case, it is obtained by solving the boundary integral equation for hi. In our opinion, consistency checks are not sufficient and comparisons with results from approximate methods are also needed. In the weakly nonlinear regime, the reference solution is the classical linear solution. However, a direct comparison with this linear solution is difficult, for two main reasons: * comparing the nonlinear simulation to a linear result can be confusing since discrepancies might be due to non- linear phenomena. For that reason, a linear version of the code has been derived that only differs from the nonlinear version by the boundary conditions that are satisfied; * usually linear results are for the steady-state response while the simulation is transient. Discrepancies might be due to long-lasting transient phenomena, to the wave gen- eration or absorption mechanism used, etc... A first step has therefore been to make proper comparisons with linear results for the transient problem in a tank of finite length. The main interest of the method being to be able to ac- count for nonlinear phenomena, a proper validation of the nonlinear response should be made. For that purpose, a comparison is performed with approximate nonlinear the- ories, such as second-order theory or shallow water theory. A final test is provided by comparisons with experi- ments. In our opinion, however, these comparisons should only be made last if it is the accuracy of the numerical Figure 2: Sloshing free surface profiles scheme that has to be evaluated. In any event, a compar- ison with a linear result appears necessary in order to be sure that nonlinear phenomena are important. It should also be kept in mind that viscous effects can play an im- portant role and are not accounted for in the simulation so that discrepancies might not be due to the inaccuracy of the numerical scheme but to an improper physical model- ing. 4 Numerical Results The purpose of this section is to show some numerical re- sults that can appropriately be compared to other theories or to experiments. Numerical instabilities encountered by many similar methodsmay appear for waves of large steepness. In this case, a 5 points smoothing algorithm similar to that used in t22] is employed. 4.1 Sloshing in a Tank The study of sloshing provides a first simple test for the accuracy of the simulation. We consider a rectangular tank and the free motion of a fluid initially out of equilibrium. One of the main advantages of this configuration is that it allows quasi-analytical solutions of the transient problem to be derived at first- and second-order A. The efficiency of the numerical scheme for the solution of the linear and nonlinear problems can therefore be directly evaluated. We consider the case of a constant initial slope of the free surface, i.e., an initial elevation given by: y = ~ d (d0.5), (11) where d is the length of the tank. We show on figure 2 the free surface profiles for d = 1 (i.e., a depth equal to the length) and ~ = 0.35. With such an initial amplitude, breaking occurs in the tank. In order to evaluate the accuracy of the simulation, we compare on figure 3 the perturbation elevation, ¢<d~2 (12) to the second-order wave elevation that has been computed quasi-analytically A. For ~ = 0.1, a good agreement is 4.O r . 2.0 -, 0.0 2.C 242 _ SECONDORDER INmAL SLOPE: 0.1O ....lNmAL SLOPE: 0.35 _ .''\ :"'\ O- , , 0.0 0.2 .--------; ."_ - 0.4 0.6 0.8 1.0 x Figure 3: Sloshing perturbation elevation

achieved, indicating that both the linear and second-order component of the wave elevation are accurately computed. For ~ = 0.35, breaking occurs and the agreement deterio- rates, suggesting that higher order effects become impor- tant. More results on this topic can be found in A. 4.2 Harmonic Motion of a Piston-Type Wave- maker The wavemaker problem is interesting because it provides a simple configuration to study most of the difficulties as- sociated with the simulation. In particular, it involves a moving free surface piercing body (the wavemaker itself) that is absent in the sloshing problem. In this section, the wavemaker motion is given by: Arty = 0 t<0 (13) tots = - 2e corset) t > 0. (14) The length of the tank is d. Note that even though the motion is harmonic for t > 0, transient phenomena are expected. 4.2.1 Linear Solution in a Tank of Finite Length For a tank of finite length, it is possible to derive a quasi- analytic solution for the linear transient problem that can be evaluated quite easily. This solution can be found in t124; an alternative solution that seems to have better conver- gence properties has also been derived A. Since this problem involves a moving rigid boundary that pierces the free surface, it allows the numerical treat- ment of the intersection point to be evaluated. Here, we compare the result of the linear simulation and the quasi- analytic result. This ensures that discrepancies are only due to numerical errors. As an example, we consider a wavemaker motion at the frequency ~ = ~r/2. The corresponding wave period and wave length are 4 and 2.5, respectively. The accuracy of the simulation has been evaluated at t = 5 in a tank of length 5. We show on table 1 the root mean square of the error (the reference 100 is taken as the error for the coarsest grid used, equal to 0.0903 at) as a function of the number of nodes per wave length, NX, and the number of time steps per wave period, NT3. VT ~ / NX 4 16 32 64 128 256 5 100 61.3 63.7 64.0 64.0 64.0 64.0 10 oo 26.2 29.3 29.9 29.9 29.9 29.9 20 oo 0.6 11.5 13.7 14.4 14.4 14.4 Table 1: Convergence table 40 . oo . 4.2 . 5.3 . 5.4 5.4 5.4 oo 1 1 1.8 1 1 2.5 1 160 on oo oo 1.3 1.4 .4 .4 3A uniform grid is used and' thanks to the symmetry, the bottom is not d~scretized. No smoothing is applied. However, using the 5 points smoothing procedure introduced in [22], even a each tine-step, does not significantly alter the accuracy of the method. 243 From this table, it appears that: * a sufficiently small time-step has to be used in order for the numerical scheme not to blow up. The influence of the time-step is otherwise rather small; * for a given spatial discretization, results converge as the time-step goes to zero. Note, however, that the numer- ical results do not converge to the exact solution; * as the spatial discretization increases, results seem to converge to the exact solution. However, the convergence rate is not second-order in the grid size. This is very likely due to the numerical treatment at the intersection that is not second-order4. 4.2.2 Linear Solution in a Tank of Infinite Length In order to test the efficiency of the damping zone, we con- sider now a tank of infinite length. For a piston type wave- maker, the solution far from the wavemaker and for large times is a progressive wave of amplitude: e k + sigh k cosh ~ (15) where k is the wavenumber associated to ~ (~2 = k tanh k). Note that in deep water (k ~ oo), a ~ at. We perform the simulation for the same wavemaker mo- tion as in 4.2.1 but in a tank of length 10 equipped at its end of an absorbing zone. We compare on figure 4 the steady state linear solution and the result of the linear simulation with and without an absorbing zone for the wave elevation at the middle of the tank (distance 5 from the wavemaker). The damping coefficient is given by (10) with cat = 1 and = 1. It appears clearly that the damping zone provides a simple and efficient way to avoid any reflection. 1.0 z 0.5 o 0.0 STEADY STATE LINEAR SOLUTION _LINEAR SIMULATION WITH ABSORPTION ....LINEAR SIMULATION WITHOUT ABSORPTION -0.5 1 1 , , , , 1 , 0- x0- Em- ~~. Kim- ~~ 6~ ]0 TIME Figure 4: Wave elevation with and without damping zone ~- 4 this indicates that numerical errors are mainly due to the numer- ical treatment at the intersection. A convergence test for the problem with periodic boundary conditions is therefore not relevant to assess the accuracy of the simulation applied to the wavemaker problem.

of - o to LLI O 1- , . , , , . , . 1 ·4C t.0 60.0 80.0 1 00.0 1 20.0 140.0 TIME Figure 5: Damping zone relative error vs. time In order to evaluate more precisely the efficiency of the damping zone, we show on figure 5 the difference between the steady state linear solution and the result of the linear simulation for the relative wave elevation (wave elevation divided by the wave amplitude ~ at) at the same point. If long-lasting low frequency oscillations appear, the relative error is only of a few percents and the absorption mecha- nism appears to be quite satisfactory. 4.2.3 Nonlinear Solution in Shallow Water In order to estimate the accuracy of the method for the nonlinear computation, we first consider the case of a shal- low water swell. Mel and Unluata t24] have explained how such a swell can experience very drastic nonlinear deformations when it propagates. Very recently, Chapalain t2] has performed experiments at the Institut de Mecanique de Grenoble that neatly confirm this bi-harmonic resonant behavior. It ap- peared therefore interesting to try and reproduce them with the simulation. The only data for the numerical simulation are the ge- ometry of the tank, the law of motion of the wavemaker and the friction coefficient f,,, used to model dissipations. The experiment was performed in a 40 cm deep tank with a piston-type wavemaker the motion of which was given by (14) with an = 15.9 cm and ~ = 2.5 rad/s. The total length of the tank, 36 m, is simulated using 300 nodes on the free surface. The simulation on 30 s took approximately 7 minutes on a CRAY-XMP. sin order to model dissipative effects in a way similar to that used by Chapalain for Boussinesq equations we substitute to Bernoulli equation: with 1 _ _ go+ 2 vet v¢+y+~=o 4 a' w v = 3 fill 2 We took the value of Jo, used by Chapalain, fw = 0.1. Note that this modeling of dissipative effects, already used in [14] for the study of sloshing, is similar to the modeling of the damping zone. _ l ~ n 1 n n 15 n 20.0 2s.0 30.0 1 1 To 1 o.o 1 5.0 20.0 2s.0 30.0 t ~ A A ~ A A A A A Al i ~ A A A ~ A A A A I Ji , ~ ~ 0.0 5.0 1 0.0 1 5.0 20.0 2s.0 30.0 a,/ ~~; A ~~ ~ ~<~ · ~ I. lo. 3°. be. 'I 0.0 5.0 1 0.0 1 S.0 20.0 2s.0 30.0 ~/~/\\,AJ\J\J\J\J\Vl\J\Jl ~~ ·~°- I I 1 1 1 ~ , , , , I a. 1~. 20. I. "~ 0.0 5.0 1 0.0 1 5.0 20.0 2s.0 30.0 . . A ~ /\ /\ ~ ~ ~ /~\ 1\ ~ f\ ~ A ~ A A ~ A ~ . \/ \J \J \J \J A \J V V \) A i A/ \/ ~ V ~ \J A ~ , Figure 6: Shallow water swell measured (left) vs. com- puted (right) wave elevations at several points in the tank We compare on figure 6 the measured and computed wave elevations at several points in the tank. An excellent agreement is achieved. This agreement is confirmed by a Fourier analysis performed once a steady state is reached. The amplitudes of the first three harmonics is plotted as a function of the distance along the tank on figure 7. It ap- pears that the simulation is very efficient to model shallow water waves and their generation by a wavemaker. 70.0 60.0 - ~ 50.0 _' c, 40.0 c' 30.0 20.0 10.0 _SINDBAD EXPERIMENTS (CHAPALAIN) 0.0 0.0 5.0 10.0 15.0 20.0 25.0 DISTANCE FROM WAVEMAKER (M) Figure 7: Shallow water swell measured vs. computed 244

4.3 Arbitrary Motion of a Piston-Type Wave- maker In deep water, nonlinear phenomena in the propagation of a regular swell are rather long to develop. Numerical methods with periodic boundary conditions, such as t10], are probably more suited to the study of this problem than the direct simulation of a wave tank. An appropriate choice for the motion of the wavemaker can however lead to wave focusing that results in break- ing. Experiments based on this principle were performed at MIT and a comparison with a numerical simulation sim- ilar to ours made see t12] where the law of motion of the wavemaker is given. The main drawback of this test case is that it involves a large amount of computer time. We have run Sindbad on the same case, but using an absorb- ing zone and a somehow coarse grid in order to minimize computational effort. The calculation has been performed on a CRAY-XMP and demanded "only" 30 minutes with 250 nodes on the free surface (compared to 30 hours with 500 nodes on a C RAY 1 for the simulation performed at MIT). ._, 0.25 0.00-~ - 0 . 25 - , 0 10 20 30 40 50 ut 0.25 - 0.00 ~ , , -0.25 0 10 20 30 40 50 60 0.25 - 0.00' ~ -0.25 a 0.25 - _ 0.00 - _ .O,,l, , , , , , 1 ~ ~ ~ 20 30 40 50 60 0.25 - 0.00 _ ~ -0.25 - _ o I I ~ I I 1 1 0 20 30 40 50 60 O.25- ' ! ` ~ O DO - - ~ -0 . 25 - , 0 10 20 30 40 50 60 Figure 8: Steep deep water waves measured (left, dashed line) vs. computed (solid line, left: MIT, right: Sindbad) wave elevations at several points in the tank 60ur own experience tends to suggest that "numerical" overturning is very sensitive to the discretization used and more particularly to the node distribution along the free surface. Here again the validity of the simulation is difficult to establish. Our interest has been mainly to perform the simulation up to the point where breaking occurs. The results of our simulation are in good agreement with both experimental and numerical results obtained at MIT up to breaking see figure 8. If overturning develops at the same time, our computation fails sooner than their (before the closing of the tube)6. This confirms that the numerical simulation can repro- duce accurately nonlinear phenomena observed experimen- taly. 4.4 Wave Diffraction on a Submerged Cyl~n- der The case of the wave diffraction on a submerged cylinder al- lows a first study of wave-structure interaction. This prob- lem has been studied extensively in the past. In partic- ular, Ursell t31] used linear theory and showed that, for a circular cylinder, there is no reflection. Ogilvie t25] ex- tended Ursell's results and computed the second-order ver- tical drift force. Very recently, Vada t32] computed the second-order potential and calculated the diffraction loads and the diffracted waves to second-order. Chaplin t3] measured diffraction loads in the labora- tory while Grue and Granlund [17] measured the diffracted waves. These experiments have partly confirmed the re- sults of first- and second-order theories. They have also exhibited some important nonlinear phenomena not ac- counted for by these theories. Such nonlinear phenomena can either be due to nonlinear free surface effects (of third- order or higher) or to viscous effects. Consequently, the present study has two main objec- tives: * to recover the results of first- and second-order theo- ries in order to assess the numerical accuracy of the method; * to compare with experimental results in order to de- termine the relative importance of viscous and free surface effects for the nonlinear phenomena observed. Fully nonlinear simulations similar to ours have been performed by several authors (e.g., t29], [11~. In general, periodic boundary conditions were used. To our knowledge, however, comparisons with second-order theory were not achieved. 4.4.1 Diffraction Loads We consider a circular cylinder of radius r = 0.06. The coordinates of its center are xc = 3.5 and Yc = -0.12. Waves are generated by a piston-type wavemaker moving at the frequency w = 1.85. The simulation is made in a tank of length 10 with 200 markers at the free surface and 60 time-steps per period. The forces acting on the cylinder are computed by di- rect integration of the pressure. In order to compare the results with those of experi- ments or of first- and second-order theory, we use a Fourier series expansion of the transient signal (once a steady-state is reached). This yields: 3 2 = F(°) + ~ F(n) cos(n`~,t + 8(n)) (16) r ~ n>1 3Y2 = F(°) + ~F(n) cos(ncot + 8(n))' (17) 24s

where the subscript x and y denote the horizontal and ver- tical components of the force, respectively. ka Kc 0.05 0.50 0.07 0.75 0.10 1.00 0.12 1.25 0.14 1.47 F(°) F(°) F(1) F(1) F(2) F(2) 0.00 0.04 1.07 1.07 0.07 0.07 0.00 0.10 1.58 1.58 0.15 0.15 -0.01 0.16 2.07 2.06 0.24 0.25 -0.02 0.24 2.52 2.51 0.33 0.34 -0.03 0.30 2.88 2.87 0.41 0.41 3.00 r 2.50 _ 2.00 1 .50 I) 1.00 Table 2: Diffraction loads Sindbad O 50 As in t3] we introduce the Keulegan-Carpenter number, I{c. For a linear deep-water wave, Kc is given by: Kc =exp~kyc). (18) Table 2 gives the values of F(n) and F(n) for n = 0, 1, 2 vs. Kc in the case just described (that corresponds to Chaplin case E t34~. Following Chaplin, we write: F(n) = ~ CxnmKc, F(n) = m>1 ~ CynmI{c m>1 The classical inertia coefficients are equal to Call et Cysts According to linear theory, these are the only non- zero coefficients and they are equal. Ogilvie [25] calculated them; they go to 2 as the immersion depth goes to infinity. We give on figure 9 the horizontal and vertical inertia coefficients vs. Kc. For a small value of the Keulegan- Carpenter number, both experimental and numerical re- sults go the value predicted by linear theory, 2.25. As I(C increases, however, the sharp decrease of the inertia coeffi- cient observed experimentally is not predicted by the sim- ulation. As argued by Chaplin, this is very likely due to viscous effects (creation of a circulation around the cylin- der). Figures 10 and 11 give the second-order vertical drift force and the force at the double frequency, respectively. Here, a good agreement appears between experiments, sec- ond-order theory and the present simulation. It should be stressed that recovering the results from second-order steady theory with the simulation is by no means obvi- ous: it demands a good accuracy and a proper control of the wave generation and absorption mechanism. It can be argued that the results of the simulation are better than those of second-order theory for large values of Kc. 4.4.2 Diffracted Waves Grue and Granlund t17] performed experiments related to incoming deepwater Stokes waves passing over a restrained submerged circular cylinder. For a small cylinder submer- gence, a strong local nonlinearity is introduced at the free surface above the cylinder and free higher order harmonic waves are generated. We remind of the observations of Grue and Granlund quency 246 0.00 ° `0~ x 00 t~ dO R oEXPERIMENTS (CHAPLIN) xSINDBAD _UNEAR THEORY (OGILVIE) . 0.00 0.50 x x - do 0 to 0 e8 0 1 .00 Kc = 1.50 2.00 Figure 9: Diffraction loads inertia coefficients (19) (20) - 1 .O ~2.0 I_ ' -3.0 4 oEXPERIMENTS (CHAPLIN) xSINDBAD '/ _ SECONDORDER (OGILVIE) ~~ x d. Fox Otto 0/ ~0 ~4 ye mix _ _ . . . . . . . 2.00 - 1.50 - 1.00 - 0.50 0.00 0.50 1.00 In(Kc) Figure 10: Diffraction loads vertical drift force 1 .0 _ - c~ ~ -2.0 cod x3.0 _' -4.0 -5.0 oEXPERIMENTS (CHAPLIN) / xSI NDBAD / _SECONDORDER (VADA) ~ Coax a' <~ don Ox y 2~00 - 1.50 - 1.00 - 0.50 0.00 0.50 1.00 In(Kc) Figure 11: Diffraction loads response at double fre-

0.06 - 0.05 - 0.04 - cat 0.03 c, 0.02 0.01 ~~ ° x x Fax 0.00- .k . . . 0.00 0.05 / oEXPERIMENTS (GRUE) / XSINDBAD / _SECONDORDER (VADA) o o . . . . . . . . 0.10 0.15 0.20 ~ k Figure 12: Diffracted waves amplitude of second-order free wave concerning the wavefield far away from the cylinder: * upstream of the cylinder: an incoming Stokes wave- train. No reflected waves, even to higher order; * downstream of the cylinder: shorter free second har- monic waves of considerable amplitude are riding on the transmitted Stokes wave. If these trends are well predicted by second-order the- ory [32i, the quantitative agreement with experiments is rather disappointing. The amplitude of the second har- monic free wave, as, only increases as the square of the amplitude of the incoming wave, a, for very small values of a. A "saturation" rapidly appears; thereafter a2 remains almost constant see figure 12. These findings suggest that the range of validity of second-order theory is quite narrow in this case. Observ- ing that this theory predicts amplitudes of the second-order free wave as large as that of the incoming wave, this should not appear as totally unexpected. In order to see if non- linear free surface effects and not viscous effect are, indeed, responsible for this deficiency, it appeared interest- ing to run Sindbad on this case. In order to compare our results with those of Grue and Granlund, we write (once a steady state is reached): * for the incident wave: Hi = a costedfit+ 8) + at cos 2(kxcat + 8) + · · ·, (21) where at is the amplitude of the second-order locked wave; * for the diffracted wave: y~ = al costedcot + 8~) + ad cos 2(kx - At + 8~ ~ + a2 cos(4kxAct+ 82~+ ·--, (22) where ad and a2 are the amplitudes of the second-order locked and free wave, respectively. In order to exhibit the second-order free wave, it is nec- essary to take a rather fine grid. The wavenumber of this free wave is, indeed, four times larger than the wavenum- ber of the incoming wave. The wavenumber of the incoming wave is chosen equal to ~ = 3.42 (so that we are in deep 3.0 ~ ....SINDBAD: LINEAR CALCULATION _SINDBAD: ok = 0.03 _SINDBAD: ok - 0.09 2.0 c, 1.0 -i 0.0 - -1 .0 2.0 . 0.0 x 3 C On Figure 13: Diffracted waves free surface profiles 5.0 water). With these choices, the results shown by Grue and Granlund correspond to a cylinder radius r = 0.117 and a depth of immersion of the center of the cylinder Yc = -0.1755. The simulation was perfomed in a tank of length 8, with a damping zone of length 3. The cylinder center was located at a distance xc = 2 from the wavemaker. 340 nodes were distributed on the free surface and 60 time steps used per period of the incoming wave. On figure 13, free surface profiles after 7 periods are shown for several amplitudes of the incident wave. The apparition of a perturbation of wavenumber 4k is obvious. However, its amplitude does not increase as the square of the incident wave amplitude. A Fourier analysis of the diffracted wave confirms this trend. We show on figure 12 the comparison between Grue and Granlund's experiments, Vada's second-order theory and the present calculation. The agreement between the numerical simulation and the experiments is very good, in- dicating that the "saturation" is, indeed, a nonlinear free surface phenomenon not accounted for by second-order the- ory. Grue and Granlund observed breaking for ah ~ 0.085, while we were able to perform the numerical simulation up to ah = 0.12. It is rather interesting to note that this does not seem to affect the amplitude of the second-order free wave. A little more surprising is the reason for which the nu- merical computation fails for ah = 0.12 after 3.2 periods. The computation does not blow up because of the overturn- ing of the crest, as would have been expected, but because of a concentration of particles just aft the cylinder. Phys- ically, it seems that particles flow very rapidly over the cylinder and are then decelerated. Here again, the validity of the simulation is difficult to establish. It is rather interesting to note that for waves passing over a submerged cylinder, nonlinear free surface effects are important for the diffracted waves but do not seem to affect very much the forces. 247

0.01 - CY o Lo .$ 0.00- C~ F > -0.01 - at_ T ' I 0 0 2.0 4.0 6.0 8.0 10.0 12.0 TIME Figure 14: Forced heaving transient force 4.5 Wave Radiation by a Free Surface Pierc- ing Cylinder As a last example, we consider the case of a free surface piercing body in forced or free motion. 4.5.1 Forced Motion Forced motions of a free surface piercing circular cylin- der have been studied quite extensively using linear and second-order theory (e.g., t26i, t19~. Fully nonlinear cal- culations have also been attempted for forced heaving of a circular cylinder by a few authors. Initially (~13i, [333) only the starting phase was considered. Recently, Hwang et al. t18] calculated the steady state response and made comparisons with first- and second-order theories and with experiments. As an example, we consider the case of a half-immerged circular cylinder with kr = 1. The simulation is performed in a tank of length 4 with a forcing frequency ~ = 3.16 (so that k = 10~. The cylinder center has for elevation: Yc = 0 t < 0 Yc = ac sin(wt) t > 0. (23) (24) Because there is no wavemaker in this case, absorbing beaches (with cat = ~ = 1) are located at each end of the tank. We used 200 nodes on the free surface. In order to avoid too small or too large segment sizes near the inter- section, a regridding procedure similar to that introduced in t11] was used when a node came too close or too far from the intersection. We show on figure 14 the transient vertical force on the cylinder as a function of time for ac = 0.5 r. If a steady state is rapidly reached, the signal is obviously not monochromatic and harmonics are present. The free surface profiles corresponding to this case are shown on figure 15. Once a steady state is reached, we write the force as: Fy + 2r Yc 0.51r r2 ac w2 y t . x ~ ' it/ f rr _/ f ~~ an-, Figure 15: Forced heaving free surface profiles The second term on the left-hand side is the linear hydro- static contribution. The acceleration-phase and velocity- phase components of the force at the forcing frequency would correspond to the added mass and damping coef- ficients for the linear problems. LIT F(°) F(l a) F(l b) Fy(2) F(3) F;,(4) 0.l -0.02 0.61 0.39 0.06 0.2 -0.03 0.62 0.38 0.12 0.02 0.3 -0.05 0.63 0.36 0.18 0.04 0.02 0.4 -0.07 0.65 0.35 0.23 0.08 0.05 0.5 -0.10 0.66 0.33 0.30 0.14 0.09 Table 3: Radiation loads Sindbad We give on table 3 the amplitude of the different har- monics as a function of ~ = aC/r. The Fourier analysis is performed on the four last periods of the signal. For small values of c, the results from linear and second-order theory are recovered (e.g., t26i, t19~. As the amplitude of the mo- tion increase, the added mass increases while the damping coefficient decreases. This behavior is in agreement with available experimental observations t19] and other numer- ical results t184. The importance of relatively high-order harmonics is quite remarkable. For ~ = 0.5, the ratio of the amplitude of the fourth harmonic to that of the first is almost equal to 15~o. This behavior, that obviously cannot be accounted for using second-order theory, is quite different from that observed for the diffraction over a submerged cylinder. It shows the interest of a fully nonlinear simulation, in partic- ular in order to assess the range of validity of approximate theories. If these results for forced motions are very promis- ing, more experimental data would be needed to make pro- prer comparisons. Note that for ~ = 0.6 the numerical computation breaks down before a steady state is reached, apparently because the cylinder goes out of the water. = F(°) + F(1a) sin(wt)F(1v) cos(wt) oo ~ F(n) sin(nc`;t + 8(n)). rl=2 (25) 7 Note that for the nonlinear problem the distinction between added-mass and hydrostatic components is somehow arbitrary. 248

2.0 1 1.0 0.0 _SINDBAD (YO = 0.9 r) LINEAR THEORY (URSELL) _ 1 1 .0 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 T Figure 16: Free heaving vertical displacement 4.5.2 Free Motion The free motion of a surface piercing circular cylinder was calculated using linear theory in [23~. A fully nonlinear computation was performed in [33] using periodic boundary conditions. A good agreement between these two theories was achieved. Here, we performed again the same calculation. The body motion is calculated in a way similar to that used in t33] (see t7] for details on the numerical implementation). At t = 0, the fluid is at rest and the cylinder, which is half-immerged at equilibrium, has an elevation: Yc(0) = Yo. (26) We show on figure 16 the subsequent vertical displace- ment of the cylindre, tic, for yo = 0.9 r. The simulation is performed exactly in the same conditions as in the preced- ing section. The absence of any nonlinear behavior for the cylinder elevation is rather striking and contrasts with the forced motion case. 5 Conclusions The MEL has been applied to the numerical simulation of a wave tank equipped at one end with a wave maker and at the other end with a damping zone. This configuration al- lows two-dimensional wave-structure interaction problems to be studied. Results have been presented that include a wide range of applications, such as wave generation and absorption, wave diffraction and wave radiation by submerged or sur- face piercing bodies in forced or free motion. Results from approximate theories (linear, second-order or shallow wa- ter theory) can be recovered, not only for short transient evolutions but also for steady state responses. Some non- linear phenomena observed experimentally and that have not been otherwise explained are also accounted for. This seems to indicate that the simulation can be used as a "standard" that allows the validity of approximate theo- ries to be assessed. Applications of the method in that direction suggest themselves. In particular, the roll motion of barges is being presently studied. 249 If the MEL provides an efficient and versatile way to study two-dimensional free surface problems, it cannot ac- count for viscous effects (except in a very crude way for instance for the modeling of bottom friction). This is par- ticularly a problem for viscous effects occurring in the near vicinity of the free surface, i.e., viscous effects associated to breaking. Breaking is a major limitation for the simu- lation because the calculation has to stop whenever a local breaking event occurs. It would therefore be most useful to be able to simulate breaking and associated dissipative effects, even in a crude way. Some hope exists (e.g., [30i, t44) for spilling break- ers, but there is an obvious need for more theoretical and experimental work on the subject. 6 Acknowleclgements This work is a result of research sponsored in pant by the "Ministere de la Defense", DRET, under convention num- ber 88/073. This support is gratefully acknowledged. References t1] BAKER, G.R., MEIRON, D.I., and ORSZAG, A., 1981, "Applications of a generalized vortex method to nonlinear free-surface flows," Third International Conference on Numerical Ship Hydrodynamics, Paris, pp. 179-191. t2] CHAPALAIN, G., 1988, "Etude hydrodynamique et sedimentaire des environnements littoraux domines par la houle," These de Doctorat de l'Universite Joseph Fourier, Grenoble (in French). t3] CHAPLIN, J.R., 1984, "Nonlinear forces on a horizon- tal cylinder beneath waves," J. Fluid Mech., Vol. 147, pp. 449-464. t4] COINTE, R., 1987, "A theory of breakers and break- ing waves," Ph. D. Dissertation, University of Califor- nia, Santa Barbara. t5] COINTE, R., 1988, "Remarks on the numerical treat- ment of the intersection point between a rigid body and a free surface," Third International Workshop on Water Waves and Floating Bodies, Woods Hole, Mass. t6] COINTE, R., 1989, "Calcul des efforts hydrody- namiques sur un cylindre horizontal en theorie potentielle non-lineaire", Deuxiemes Journees de l'Hydrodynamique, Nantes, pp. 89-104 (in French). t7] COINTE, R., 1989, "Quelques aspects de la simulation numerique d'un canal a houle", These de Doctorat de l'Ecole Nationale des Ponts et Chaussees, Paris (in French). t8] COINTE, R., JAMI, A., and MOLIN, B., 1987, "Nonlinear impulsive problems," Second International Workshop on Water Waves and Floating Bodies, Bris- tol. [9] COINTE, R., MOLIN, B., and NAYS, P., 1988, "Non- linear and second-order transient waves in a rectangu- lar tank," BOSS'88, Trondheim.

t10] DOLD, J.W., and PEREGRINE, D.H., 1986, "Water- wave modulation," Bristol University, Report AM-86- 03. t11] DOMMERMUTH, D.G., 1987, "Numerical methods for solving nonlinear water-wave problems in the time domain," Ph.D. Dissertation, MIT. t12] DOMMERMUTH, D.G., YUE, D.K.P., RAPP, R.J., CHAN, E.S., and MELVILLE, W.K., 1988, "Deep- water plunging breakers: a comparison between poten- tial theory and experiments," J. Fluid Mech., Vol. 189, pp. 423-442. t13] FALTINSEN, O.M., 1977, "Numerical solutions of transient nonlinear free surface motion outside or in- side moving bodies," 2nd International Conference on Numerical Ship Hydrodynamics, Berkeley, pp. 347- 357. t14] FALTINSEN, O.M., 1978, "A numerical method of sloshing in tanks with two-dimensional flow," Journal of Ship Research, Vol. 22, pp. 193-202. t15] GREENHOW, M., 1987, "Wedge entry into initially calm water," Applied Ocean Research, Vol. 9, pp. 214- 223. [16] GRISVARD, P., 1984, "Problemes aux limites darts les polygones mode d'emploi", Prepublication du laboratoire de mathematiques de l'Universite de Nice (in French). [17] GRUE, J., and GRANLUND, K., 1988, "Impact of nonlinearity upon waves traveling over a submerged cylinder," Third International Workshop on Water Waves and Floating Bodies (Woods Hole). t18] HWANG, J.H., KIM, Y.J., and KIM, S.Y., 1988, "Nonlinear Hydrodynamic Forces Due to Two- dimensional Forced Oscillation," Proceedings, IU- TAM Symposium on Nonlinear Water Waves (Tokyo), Springer-Verlag, pp. 231-238. [19] KYOZUKA, Y., 1982, "Experimental study of second- order forces acting on a cylindrical body in waves," Proceedings, 14th Symposium Naval Hydrodynamics, Ann Arbor. t20] KRAVTCHENKO, J., 1954, "Remarques sur le Calcul des Amplitudes de la Houle Lineaire Engendree par un Batteur", 5th Conf. Coastal Eng., Grenoble, Fiance, pp. 50-61 (in French). t21] LIN, W.M., 1984, "Nonlinear motion of the free sur- face near a moving body," Ph. D. Dissertation, MIT, Cambridge, Mass. [22] LONGUET-HIGGINS, M.S., and COKELET, E.D., 1976, "The deformation of steep surface waves on water. 1. A numerical method of comptation," Proc. R. Soc. London, Vol. A 364, pp.1-26. t23] MASKELL, S.J., and URSELL, F., 1970, "The tran- sient motion of a floating body," J. Fluid Mech., Vol. 44, pp. 303-313. t24] MEI, C.C., and UNLUATA, U., 1972, Harmonic gen- eration in shallow water waves," Waves on beaches, R.E. Meyer, Ed., Academic Press, pp. 181-202. t25] OGILVIE, T.F., 1963, "First- and second-order forces on a cylinder submerged under a free surface," J. Fluid Mech., Vol. 16, pp. 451-472. t26] PAPANIKOLAOU, A., and NOWACKI, H., 1980, "Second-order theory of oscillating cylinders in a reg- ular steep wave," Proceedings, 13th Symposium Naval Hydrodynamics, Tokyo. t27] PEREGRINE, D.H., 1972, Unpublished note. t28] ROBERTS, A.J., 1987, "Transient Free Surface Flows Generated by a Moving Vertical Plate," Q. J. Mech. Appl. Math., Vol. 40, Pt. 1. [29] STANSBY, P.K., and SLAOUTI, A., 1984, "On non- linear wave interaction with cylindrical bodies: a vor- tex sheet approach," Appl. Ocean Res., Vol. 6, pp. 108- 1 15. [30] TULIN, M.P., and COINTE, R., 1987, "Steady and unsteady spilling breakers: theory," Proceedings, IU- TAM Symposium on Nonlinear Water Waves, Tokyo, pp. 159-165. [31] URSELL, F., 1950, "Surface waves in the presence of a submerged circular cylinder, I and II," Proc. Camb. Phil. Soc., Vol. 46, pp. 141-158. t32] VADA, T., 1987, "A numerical solution of the second- order wave-diffraction problem for a submerged cylin- der of arbitrary shape," J. Fluid Mech., Vol. 174, pp. 23-37. t33] VINJE, T., and BREVIG, P., 1981, "Nonlinear, two- dimensional ship motions," Norwegian Institute of Technology, Rapport R-112.81. 250