- Email: [email protected]
Numerical simulation of a wave channel
Numerical simulation of a wave channel Raymond Cointe Institut Francais du P~trole, Rueil-Malmaison, France Presently: Bassin arEssais des CarOnes, 8, Bd Victor, 75732 Paris Cedex 15, France The application of the Mixed Eulerian-Lagrangian method to the simulation of transient free surface flows in the vicinity of a free surface piercing structure is considered. Particular 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 approximate theories or experiments are shown. They demonstrate 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. Key words: free surface flow, nonlinear water waves, numerical simulation, potential flow.
1. I N T R O D U C T I O N The direct numerical simulation of unsteady two-dimensional potential free surface flows using a Mixed Eulerian-Lagrangian method (MEL) has received considerable attention since the pioneering work of LonguetHiggins and Cokelet 1. 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 2, simulation of breaking in a tank 3, forced heaving of a cylinde?, etc...). Compared to the initial application of the MEL to the study of steep periodic waves (e.g., Ref. 1), several new difficulties 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 accuracy 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 computed and quantitative validations of the numerical procedure are almost non-existent, even for weakly 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 developed for the purpose of simulation of a two-dimensional wave tank using potential flow theory. Particular 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
© 1990Computational Mechanics Publications
comparison easier, an option of the code allows the linear problem 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 pistontype wavemaker. These waves can then interact with submerged or free surface piercing cylinders in forced or free motion. Reflection from the wall opposite tQ the wavemaker is avoided by the use of a damping zone--see Fig. 1. The results of the classical first- and second-order diffraction radiation theories can be recovered using the simulation. 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.
Y
2D O u r n r t ' / ¢ o l
toov~
I g-1
l'onR
x
1
E Plston-WPe wavemaker
Tested body
"Beach" (absorbing zone]
Fig. 1. Sindbad---geometric definitions
Enoineerino Analysis with Boundary Elements, 1990, 1Iol. 7, No. 4
167
Numerical simulation of a wave channel: R. Cointe This work has been reported elsewhere while in progress 5'6'7'a. 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 Ref. 9.
3. N U M E R I C A L 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 difliculties 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 yaxis is upward vertical--see Fig. 1 for geometric definitions. The fluid is assumed to be incompressible and the flow irrotational so that the velocity field ~ is given by: = V~,
(1)
with Aq~ =
O.
(2)
The computation is performed in a bounded domain. Along rigid boundaries (2 ~ Fn(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 2 on the free surface (2 ~ Fn(t)) and the associated value of the potential, t~(2). This yields:*
D~= - y -
- ~ qb~ +
D2 O~ = ~bsg + q~"fi'
02,-p+c(t)
(3) (4)
where D is used to indicate a material derivative, g and h are vectors tangent and normal to the free surface, respectively, and ( is an arbitrary constant. 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 c(t). With an appropriate choice of the velocity potential, this function can be taken equal to zero. At a given instant t, if 4, and ~b, are known along the free surface Fd(t), then the right-hand sides of (3)-(4) can be evaluated. The fact that ~b is harmonic in the fluid
*For the sake of simplicity, we use units such that the acceleration of gravity, g, the specific mass of water, p, and the depth of the tank, h, are equal to 1.
168
domain allows the value of q~, along Fa(t) to be computed from the values of q~ along Fa(t ) and of ~b, along F,(t). The kinematical constraint A~b = 0, associated with the boundary condition on F,, permits therefore to express the free surface boundary conditions (3)-(4) as an evolution equation for (qS, 2). 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 ~b knowing 4~ along Fd(t) and ~b, along F,(t).We use the integral equation:
--O(P)dp(P) + ~
dp(Q)G.(P,Q) dsQ
d r d+I'n
=f
gp.(Q)G(P,Q) dsQ,
(5)
J r a+r'n
where P is a point on the boundary, G is the Green function, 0(P) the angle between two tangents of the boundary at P (equal to rr for a smooth curve) and s a curvilinear abscissa along F. Equation (5) is discretized using a standard collection method. The boundary of the domain is approximated by segments and 4~ and ~b, are assumed to vary linearly along each segment. This allows an analytical 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 F.(t)]--Dirichlet [along Fd(t)] boundary value problem. It is well known that the solution of such a "mixed" problem is singular at the intersection points, i.e., at points belonging both to F.(t) and Fd(t)--e.g. Ref. 10. Let us consider, as an example, the wavemaker problem. For a piston type wavemaker and an horizontal free surface (90 degrees intersection angle), the complex potential W = q~ + i qJ 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:
4o(2, t) = d?o(2) + tchl(x) + . . .
(6)
Not surprisingly, performing such an expansion also leads to a z log z behavior for ~bl (e.g. Refs. 11 and 12). 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 nonlinear reoime, i.e. here, for a small acceleration of the wavemaket (relative to that of gravity) s. 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 condition; it is given by: q~,, + ~ = 0.
Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
(7)
Numerical simulation o f a wave channel: R. Cointe
If regularity in time is assumed, the local behavior of the complex potential for this solution can be shown to be in z 2 log z. Not surprisingly, this singularity is the same as that appearing for the harmonic problem (for which q~ = tp exp(icot)) which was studied by Kravtchenko 13. It is much weaker than that of the problem discretized in time. The singularity can also be studied in a similar fashion for an arbitrary angle of intersection 0. Assuming that the complex potential is bounded near the intersection, the leading behavior of the complex potential can be shown to be in z ~/° or z ~/° log z if n/O is an integer.* A consequence is that, for an intersection angle smaller than re, ~b is continuous along the boundary F while ~bs and tk, are piecewise continuous but experience a finite jump at the intersections. The numerical treatment at the intersection has been devised to accommodate such a behavior, referred to as weakly singular. More details on the numerical implementation can be found in Refs. 5 and 9. These results only apply in the weakly nonlinear regime, i.e., when the classical linear solution yields a first reasonable approximation to the problem. Obviously, it would be interesting to study other regimes. The impulsive regime (very large acceleration of the wavemaker) has been studied by several authors. As expected, the local behavior of the classical linear solution cannot yield any information for the nonlinear problem in this casC 4. It seems that Peregrine's solution 1~ corresponds to an outer solution to which an inner solution should be matched. Of course, only this inner solution is relevant to the study of the local behavior. If work has been carried out to find an inner solution 7'12, it leads to serious difficulties and, to our knowledge, 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 9.
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 diffracted 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 nonlinear 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 piston-type wavemaker. A "beach" is used for the absorption of the waves that are generated in the tank--see Fig. 1. It is in fact a damping zone, similar to that used in Ref. 15. In this zone, the free surface boundary conditions are modified by adding a *The same singularity occurs at corners of rigid boundaries.
damping term. We, therefore, write:
oo ot
D~ Dt
-
Y -
-
)
1
~ c~ + ~ c~ -
v(x)(c~ -
~,)
- (q~fi + c~,h - v(x)(5c - fee),
(8) (9)
where the subscript e corresponds to the reference configuration for the fluid. The function v(x) is homogeneous 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 v(x) is equal to zero except in the damping zone (Xo < x < xl), is chosen continuous and continuously differentiable, and is "tuned" to a characteristic wave frequency o9 and to a characteristic wave number k: v(x) = ctco
X - Xo
,Xo < X < Xl = Xo +
.(10)
If the proper scales have been chosen, values of ~t and fl 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 numerical 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 energy 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 boundaries is computed by two different methods. Both of them use Bernoulli's equation, but they differ by the evaluation of the ~bt term. In one case, this term is evaluated by finite differences in time. In the other case, it is obtained by solving the boundary integral equation for ~br 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 nonlinear 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 generation or absorption mechanism used, etc... A first step has therefore been to make proper comparisons
Engineerino Analysis with Boundary Elements, 1990, Vol. 7, No. 4
169
Numerical simulation of a wave channel." R. Cointe with linear results for the transient problem in a tank of finite length. The main interest of the method being to be able to account for nonlinear phenomena, a proper validation of the nonlinear response should be made. For that purpose, a comparison is performed with approximate nonlinear theories, such as second-order theory or shallow water theory. A final test is provided by comparisons with experiments. In our opinion, however, these comparisons should only be made last if it is the accuracy of the numerical scheme that has to be evaluated. In any event, a comparison 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 important 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 modeling.
'q;
/ F~.2. Sloshing--free surface profiles
4. NUMERICAL RESULTS The purpose of this section is to show some numerical results that can appropriately be compared to other theories or to experiments. Numerical instabilities--encountered by many similar methods--may appear for waves of large steepness. In this case, a 5 points smoothing algorithm similar to that used in Ref. 1 is employed.
4.2. Harmonic motion of a piston-type wavemaker The wavemaker problem is interesting because it provides a simple configuration to study most of the difficulties associated 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: l(t) = 0
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:
al
l(t) = - 2 cos(wt)
t >_ 0.
(14)
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 Ref. 3; an alternative solution that seems to have better convergence properties has also been derived s. Since this problem involves a moving rigid boundary that pierces the free surface, it allows the numerical 4.0_SECOND-ORDER __INITIAL SLOPE: 0.10 .... INITIAL SLOPE: 0.,35 2.0" O.
Y - - Ylinear
(~d) 2
(13)
The length of the tank is d. Note that even though the motion is harmonic for t > 0, transient phenomena are expected.
where d is the length of the tank. We show on Fig. 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 Fig. 3 the perturbation elevation, y~ -
t<0
(12)
0.0. to the second-order wave elevation that has been computed quasi-analyticallya. For ( = 0.1, a good agreement is achieved, indicating that both the linear and secondorder component of the wave elevation are accurately computed. For ( = 0.35, breaking occurs and the agreement deteriorates, suggesting that higher order effects become important.
--2.0
0.0"0:2'0:4"0'.6 X
Fig. 3.
Sloshing--perturbation elevation
170 Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
0:8
1.0
Numerical simulation of a wave channel." R. Cointe Table 1. Convergencetable NT~ NX --,
5
10
• 1.0 20
40
80
160 oo oo oo 1.3 1.4 1.4 1.4
4
I00
oo
oo
oo
oo
8 16 32 64 •28.
61.3 63.7 64.0 64.0 64.0 64.0
26.2 29.3 29.9 29.9 29.9 29.9
10.6 11.5 13.7 14.4 14.4 14.4
oo 4.2 5.3 5.4 5.4 5.4
oo 1.8 2.3 2.5 2.5 2.5
_.STEADY STATE LINEAR SOLUTION _LINEAR SIMULATION WITH ABSORPTION ....LINEAR SIMULATION WITHOUT ABSORPTION
z o
0.5 •
W
0.0
~
":
:-
::.
.... -,,"
-0.5
treatment 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 co = n/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 az) as a function of the number of nodes per wave length, N X , and the number of time steps per wave period, NT*. 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 numerical 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-order.** 4.2.2. Linear solution in a tank o f infinite length In order to test the efficiency of the damping zone, we consider now a tank of infinite length. For a piston type wavemaker, the solution far from the wavemaker and for large times is a progressive wave of amplitude: sinh 2 k a = al k + sinh k cosh k"
(15)
where k is the wavenumber associated to co (092 = k tanh k). Note that in deep water (k ---, oo), a = a~. We perform the simulation for the same wavemaker motion as in 4.2.1. but in a tank of length 10 equipped at its end of an absorbing zone. We compare on Fig. 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 ~t = 1 and ~ = 1. It appears deafly that the damping *A uniform grid is used and, thanks to the symmetry, the bottom is not discretized. No smoothing is applied. However, using the 5 points smoothing procedure introduced in Ref. 1, even at each time-step, does not significantly alter the accuracy of the method. **This indicates that numerical errors are mainly due to the numerical 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.
!
i
i
i
i
i
!
TlgE Wave elevation with and without damping zone
~.4.
zone provides a simple and efficient way to avoid any reflection. In order to evaluate more precisely the efficiency of the damping zone, we show on Fig. 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 = a,) at the same point. If long-lasting low frequency oscillations appear, the relative error is only of a few percents and the absorption mechanism 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 shallow water swell. Mei and Unliiata 16 have explained how such a swell can experience very drastic nonlinear deformations when it propagates. Very recently, Chapalain 17 has performed experiments at the Institut de Mrcanique de Grenoble that neatly confirm this bi-harmonic resonant behavior. It appeared therefore interesting to try and reproduce them with the simulation. The only data for the numerical simulation are the geometry of the tank, the law of motion of the wavemaker and the friction coefficient fw used to model
0.1 iz o ix: tw hi 0.0 LIJ
0(.
-0.1 40.0 Fig. 5.
66.0
86.0 " 100.0 1 2 0 . o "rIME Damping zone---relative error vs. time
140.0
Engineering Analysis with Boundary Elements, 1990, Iiol. 7, No. 4
171
Numerical simulation o f a wave channel; R. Cointe
B ^t~
~A
.,'.^
V
'J
V
^
'J
.,~.e
@. ;st.@
I.
;~¢*cI| Oo
..-, A
'.
A
r, A A A A
+..'\. ,,
,,,,,,~
& - , ~ . @ 1i
.;~.¢ I @.
dr.;
.
:
S.
I@*
$,@°@
IS.
A A AA
@°
""~v',Ov
o~¢°@
IB..
I~..
I~.
I$*
i,.
A A AAAAAI
I
v' v v v '-' V V k./~
*1~°@ . ; Il~c.@ ~ . @ I |d.@ .
I~
II. ll.~.@
lh
II.
liB.
II.
I~.
llo
l~.
/I"A'vAvAvA : , A'AAJvAvA"
lo~.CI .~¢°@ * ~I,P..I
:
.
ll.
I$*
:
.
I~*
Ig.
I ¢ . 0 .i,
l.
~-°
.A./~,_AA A
/~
Im.
I'I.
liB.
A A A A A A A A
*+$1;.@ I'*0 O.
I*
II*
Ill
Io1
'~ v
V V
II,
III1,
II.
'I~*
IS.
ll~.
I@~.@
" " t~ "~" z~ @
~
V
~
V
OI
11
ill°
I1~.@
V ~, v -
I~,
ll;l,
/
I1~.
Itl*° •
I
S*
I
ii11
[ Ile
l.a.
,
,
A4
A A A A /~ A /~ A A A A /~ A V / V ~ 'V' V V v V -V 'V 'V' kJ V V
.
-~.e
~',J
,
* ' * ~. . . . ~". 0 ~ .
v
Io@ •
I
!
~
!
III11
II~I. ?[,'5~S 1~ S*
l~.
lrJ.
-I
lo,
C
5
:0
:5
20 1EK)S
~5
30
~5
LO
Fig. 6. Shallow water swell--measured (left) vs. computed (rioht) wave elevations at several points in the tank dissipationt. The experiment was performed in a 40 cm deep tank with a piston-type wavemaker the motion of which was given by (14) with a I = 15.9 cm and co = 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.
tin order to model dissipative effects in a way similar to that used by Chapalain for Boussinesq equations we substitute to Bernoulli equation:
We compare on Fig. 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 Fig. 7. It appears that the simulation is very efficient to model shallow water waves and their generation by a wavemaker.
4.3. Arbitrary motion of a piston-type wavemaker 1
+
~, + _- v 6 . ¢ ~ + y + vO = 0 2
with 4
~1o)
We took the value of f, used by Chapalain, f~ = 0.1. Note that this modeling of dissipative effects, already used in Ref. 18 for the study of sloshing, is similar to the modeling of the damping zone.
172
In deep water, nonlinear phenomena in the propagation of a regular swell are rather long to develop. Numerical methods with periodic boundary conditions 19, 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 breaking. Experiments based on this principle were performed at MIT and a comparison with a numerical
Engineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
Numerical simulation of a wave channel: R. Cointe 025
70.0
60.0
oooI
_SINDBAD __EXPERIMENTS
(CHAPALAIN)
-0.~
I
I
/
I
I
0
10
Eo
30
40
50
o 30.0
0.00
~ 20.0.
-0.~ t 0
10
m
50.0
oo
40.0
0 Z
10.0. 0.0
0.0
510
16.0
1,.4.0 26.0
DISTANCE FROM WAVEMAKER(M)
25.0
, 30, 40,
,
20
50
60
0.25
I
0.00
I
-
-0.~
~
I
I
I
I
I
0
lO
20
30
40
50
0
~o
20
0
1
20
60
Fi9. 7. Shallow water swell--measured vs. computed Fourier components simulation similar to ours made3--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 absorbing zone and a somewhat 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 CRAY 1 for the simulation performed at MIT). The results of our simulation are in good agreement with both experimental and numerical results obtained at MIT up to breaking--see Fig. 8. If overturning develops at the same time, our computation fails sooner than theirs (before the closing of the tube)**. This confirms that the numerical simulation can reproduce accurately nonlinear phenomena observed experimentaly.
4.4. Wave diffraction on a submerged cylinder The case of the wave diffraction on a submerged cylinder allows a first study of wave-structure interaction. This problem has been studied extensively in the past. In particular, Ursell 2° used linear theory and showed that, for a circular cylinder, there is no reflection. Ogilvie 2. extended Ursell's results and computed the second-order vertical drift force. Very recently, Vada 22 computed the second-order potential and calculated the diffraction loads and the diffracted waves to second-order. Chaplin 23 measured diffractions loads in the laboratory while Grue and Granlund 24 measured the diffracted waves. These experiments have partly confirmed the results of first- and second-order theories. They have also exhibited some important nonlinear phenomena not accounted 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 objectives: , to recover the results of first- and second-order theories in order to assess the numerical accuracy of the method;
t"\..
O.O0 -0.~
I
"I °2'i
I
I
30
I
t
I
I
40
50
i
I
i
30
40
50
---
0
i
10
i
N
I
3O
i
I
40
50
6O
Fi9. 8. Steep deep water waves--measured (left, dashed line) vs. computed (solid line, left: MIT, rioht : Sindbad) wave elevations at several points in the tank
, to compare with experimental results in order to determine 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 25'26. 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 m = 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 direct integration of the pressure. In order to compare the results with those of experiments 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: Fx - F~ ) + ~ F~ I cos(n~ot + 0<"))
**Our 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.
60
r3(o 2
F, _ F~O) + ~ F~,)cos(n~ot + 0<")),
F3(D2
(16)
n_>l
(171
n>l
Enoineerino Analysis with Boundary Elements, 1990, Vol. 7, No. 4
173
Numerical simulation of a wave channel." R. Cointe o.o
where the subscript x and y denote the horizontal and vertical components of the force, respectively. As in Ref. 23 we introduce the Keulegan-Carpenter number, Kc. For a linear deep-water wave, K~ is given by: rca Kc = - - exp(kyc).
(18)
r
Table 2 gives the values of F~ ) and F~r") for n = 0, 12, vs. K~ in the case just described (that corresponds to Chaplin case E23). Following Chaplin, we write: F~)= ~ Cx,,,K ~,
-I.O
fJ'
*EXPERIMENTS (CHAPLIN) xSINDBAD . ~ x _SECOND-ORDER (OGILVIE) • x
-
.--2.0 0 -¢ -3.0
-4.0 -5.0
(19)
....................
-2.00
-1.50
-1.00
m_>l
Cy,,,K~,
Z
F~" ) =
(20)
• ........ -0.50
t.CK¢)
0.00
0.50
1.00
Fig. 10. Diffraction loads--vertical drift force
m_>l
The classical inertia coefficients are equal to C=11 et
Cyxl. According to linear theory, these are the only non-zero coefficients and they are equal. Ogilvie21 calculated them; they go to 2 as the immersion depth goes to infinity. We give on Fig. 9 the horizontal and vertical inertia coefficients vs. K~. For a small value of the KeuleganCarpenter number, both experimental and numerical results go the value predicted by linear theory, 2.25. As K¢ increases, however, the sharp decrease of the inertia coefficient observed experimentally is not predicted by the simulation. As argued by Chaplin, this is very likely due to viscous effects (creation of a circulation around the cylinder). 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, second-order theory and the present simulation. It
Table 2. Diffraction loads-- Sindbad ka
Kc
F(°) =x
F(°) --y
0.05 0.07 0.10 0.12 0.14
0.50 0.75 1.00 1.25 1.47
0.00 0.00 -0.01 0.02 -0.03
0.04 0.10 0.16 0.24 0.30
F(I) =x
Fit) --y
F(2) --x
F(2) =y
1.07 1.58 2.07 2.52 2.88
1.07 1.58 2.06 2.51 2.87
0.07 0.15 0.24 0.33 0.41
0.07 0.15 0.25 0.34 0.41
should be stressed that recovering the results from second-order steady theory with the simulation is by no means obvious: 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 24 performed experiments related to incoming deepwater Stokes waves passing over a restrained submerged circular cylinder. For a small cylinder submergence, a strong local nonlinearity is introduced at the free surface above the cylinder and free higher order harmonic waves are generated. We are reminded of the observations of Grue and Granlund concerning the wavefield far away from the cylinder: , upstream of the cylinder: an incoming Stokes wavetrain. No reflected waves, even to higher order; downstream of the cylinder: shorter free second harmonic waves of considerable amplitude are riding on the transmitted Stokes wave. If these trends are well predicted by second-order theory 22, the quantitative agreement with experiments is rather disappointing. The amplitude of the second harmonic free wave, a2, 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 Fig. 12. o.o
3.00
2.50 g,
-1.0 on B ' ~ o ~1=~, =.x
2.00"
x
x
8
x
*EXPERIMENTS (CHAPLIN) xSINDBAD _SECOND-ORDER (VAOA)
.J
v¢ - -2.0
o
- 1.50'
=g,0
t'N
~. -3.o
c~ 1.00 • EXPERIMENTS (CHAPLIN)
0.50
-4.0
xSINDBAD
_UNEAR THEORY(OGILVlE) -5.0
0 . 0 0
.
0.00
Fig. 9.
174
.
.
.
,
0.50
.
.
.
.
~
1.00 Kc
.
.
.
.
i
1.50
Diffraction loads--inertia coefficients
.
.
.
.
2.00
. . . .
-2.00
,
. . . .
-1.50
,
. . . .
-1.00
,
. . . .
-0.50
,
. . . .
0.00
,
. . . .
0.50
1.00
In(Ke)
Fio. 11. Diffraction loads-- response at double frequency
Engineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
Numerical simulation of a wave channel." R. Cointe 3.0
0.06
.... SINDBAD:
0.05
2.0
_SECOND-ORDER (¥ADA)
UNEAR CALCULATION _SINDBAD: ak = 0.03 _SINDBAD: ok = 0.09
0.04 -
1.0-
e~ 0.03-
o
0.02-
/
...~.~°
o
/,t
0.0-
iq
0.01
0.00
0.00
0.05
0.10 ok
0.15
0.20
Fig. 12. Diffracted waves--amplitude o f second-order free wave These findings suggest that the range of validity of second-order theory is quite narrow in this case. Observing that this theory predicts amplitudes of the secondorder free wave as large as that of the incoming wave, this should not appear as totally unexpected. In order to see if nonlinear free surface effects--and not viscous effect --are, indeed, responsible for this deficiency, it appeared interesting 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: tit = a cos(kx - cot + 0) + a t cos 2(kx - cot + O) + ...,
(21)
where a t is the amplitude of the second-order locked wave; , for the diffracted wave: qa = az cos(kx - cot + 01) + a~ cos 2(kx - cot + 01) + a2 cos(4kx - 2cot + 02) + " ' ,
(22)
where a~ 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 necessary to take a rather fine grid. The wavenumber of this free wave is, indeed, four times larger than the wavenumber of the incoming wave. The wavenumber of the incoming wave is chosen equal to k = 3.42 (so that we are in deep 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 y~ = - 0.1755. The simulation was performed 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 Fig. 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 Fig. 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, indicating that the "saturation" is, indeed, a nonlinear
-1.0-2.0
0.0
I'.o
2'.0
3'.0
4'.0
5.0
X
Fi9. 13.
Diffracted waves--free surface profiles
free surface phenomenon not accounted for by secondorder theory. Grue and Granlund observed breaking for ak ~- 0.085, while we were able to perform the numerical simulation up to ak = 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 numerical computation fails for ak = 0.12 after 3.2 periods. The computation does not blow up because of the overturning of the crest, as would have been expected, but because of a concentration of particles just after the cylinder. Physically, 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 the forces very much. 4.5. Wave radiation by a free surface piercino 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 cylinder have been studied quite extensively using linear and second-order theory 27'2s. Fully nonlinear calculations have also been attempted for forced heaving of a circular cylinder by a few authors. Initially 29'3° only the starting phase was considered. Recently, Hwang et al 4. 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 co = 3.16 (so that k = 10). The cylinder center has for elevation: Yc = 0
t <_ 0
Yc = ac sin(cot)
(23) t >_ 0.
(24)
Because there is no wavemaker in this ease, absorbing beaches (with • = fl = 1) are located at each end of the tank. We used 200 nodes on the free surfaced. In order to avoid too small or too large segment sizes near the
Enoineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
175
Numerical simulation of a ,,ave channel: R. Cointe 0.01
Y
hi 0fY
o
0.00 o oe.
-0.01
).0
2'.0 ' 4'.0
o:o
1().0' 12.0
"I3~E Fig. 14.
intersection, a regridding procedure similar to that introduced in Ref. 26 was used when a node came too close or too far from the intersection. We show on Fig. 14 the transient vertical force on the cylinder as a function of time for a~ = 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 Fig. 15. Once a steady state is reached, we write the force as: Fr + 2ry, - F~°) + --rF (t°) sin(tnt) -- --rF(lV)cos(ot) 0.5Irr2 acO92 + ~ F~") sin(n~t + 0(")).
(25)
n=2
The second term on the left-hand side is the linear hydrostatic contribution. The acceleration-phase and velocity-phase components of the force at the forcing frequency would correspond to the added mass and damping coefficients for the linear problem'ft. We give on Table 3 the amplitude of the different harmonics as a function of ~ = a Jr. The Fourier analysis is performed on the four last periods of the signal. For small values of e, the results from linear and second-order theory are recovered 27'28. As the amplitude of the motion increase, the added mass increases while the damping coefficient decreases. This behavior is in agreement with available experimental observations 2s and other numerical results 4. 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 particular in order to assess the range of validity of approximate theories. If these results for forced motions are very promising, more experimental data would be needed to make proper comparisons. Note that for e = 0.6 the numerical computation breaks down before a steady state is reached, apparently because the cylinder goes out of the water.
Table 3. Radiationloads--Sindbad
0,1 0,2 0.3 0.4 0.5
?tNote that for the nonlinear problem the distinction between addedmass and hydrostatic components is somehowarbitrary.
F~°)
F~'")
V~'b)
F(2) _,
F(3) -,
F(r4)
-0.02 -0.03 -0.05 -0.07 -0.10
0.61 0.62 0.63 0.65 0.66
0.39 0.38 0.36 0.35 0.33
0.06 0.12 0.18 0.23 0.30
0.02 0.04 0.08 0.14
0.02 0.05 0.09
4.5.2. Free motion The free motion of a surface piercing circular cylinder was calculated using linear theory in Ref. 31. A fully nonlinear computation was performed in Ref. 30 using periodic boundary conditions. A good agreement between these two theories was achieved. Here, we again performed the same calculation. The body motion is calculated in a way similar to that used in Ref. 30 (see Ref. 9 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 Fig. 16 the subsequent vertical displacement of the cylinder, Yc, for Yo = 0.9 r. The simulation is performed exactly in the same conditions as in the preceding section. The absence of any nonlinear behavior for the cylinder elevation is rather striking and contrasts with the forced motion case.
2.0 -~
o>-. 0 >.-
b
1.0
_SINDBAD(YO = 0.9 r) __LINEAR THEORY (URS~'LL)
0.0
- 1 . 0
176
Forced heaving--free surface profiles
Fig. 15.
Forced heaving--transient force
'. . . . 6.0 8.0 10.0 12.0 1 LO T Free heaving--vertical displacement
0.0
Fig. 16.
Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
. . . . . . 2.0 4.0
,
,
•
',
•
N u m e r i c a l simulation o f a wave channel: R. Cointe
5. C O N C L U S I O N S T h e M E L has been a p p l i e d to the n u m e r i c a l s i m u l a t i o n of a wave t a n k e q u i p p e d at one end with a wave m a k e r a n d at the o t h e r e n d with a d a m p i n g zone. This configura t i o n allows t w o - d i m e n s i o n a l w a v e - s t r u c t u r e i n t e r a c t i o n p r o b l e m s to be studied. Results have been p r e s e n t e d t h a t include a wide r a n g e of a p p l i c a t i o n s , such as wave g e n e r a t i o n a n d a b s o r p t i o n , wave diffraction a n d wave r a d i a t i o n b y s u b m e r g e d o r surface piercing b o d i e s in forced o r free m o t i o n . Results from a p p r o x i m a t e theories (linear, s e c o n d - o r d e r o r shallow w a t e r t h e o r y ) can be recovered, n o t o n l y for s h o r t transient e v o l u t i o n s b u t also for s t e a d y state responses. Some nonlinear phenomena observed experimentally a n d t h a t have n o t been otherwise e x p l a i n e d are also a c c o u n t e d for. This seems to indicate t h a t the s i m u l a t i o n can be used as a " s t a n d a r d " that allows the validity of a p p r o x i m a t e theories to be assessed. A p p l i c a t i o n s of the m e t h o d in t h a t d i r e c t i o n suggest themselves. I n p a r t i c u lar, the roll m o t i o n of barges is being presently studied. If the M E L p r o v i d e s a n efficient a n d versatile w a y to s t u d y t w o - d i m e n s i o n a l free surface p r o b l e m s , it c a n n o t a c c o u n t for viscous effects (except in a very c r u d e w a y for instance for the m o d e l l i n g of b o t t o m friction). This is p a r t i c u l a r l y a p r o b l e m for viscous effects o c c u r r i n g in the n e a r vicinity of the free surface, i.e., viscous effects associated to breaking. B r e a k i n g is a m a j o r l i m i t a t i o n for the s i m u l a t i o n because the c a l c u l a t i o n has to s t o p w h e n e v e r a local b r e a k i n g event occurs. It w o u l d therefore be m o s t useful to be able to simulate b r e a k i n g a n d a s s o c i a t e d dissipative effects, even in a crude way. S o m e h o p e exists 32'33 for spilling breakers, b u t there is an o b v i o u s need for m o r e theoretical a n d e x p e r i m e n t a l w o r k on the subject.
6. A C K N O W L E D G E M E N T S This w o r k is a result of research s p o n s o r e d in p a r t by the " M i n i s t r r e de la Drfense," D R E T , u n d e r c o n v e n t i o n n u m b e r 88/073. This s u p p o r t is gratefully a c k n o w l e d g e d .
REFERENCES 1 Longuet-Higgins, M. S., and Cokelet, E. D. The deformation of steep surface waves on water. 1. A numerical method of computation, Proc. R. Soc. London, 1976, A 364, 1-26 2 Greenhow, M. Wedge entry into initially calm water, Applied Ocean Research, 1987 9, 214-223 3 Dommermuth, D. G., Yue, D. K. P., Rapp, R. J., Chan, E. S. and Melville, W. K. Deep-water plunging breakers: a comparison between potential theory and experiments, J. Fluid Mech, 1988, 189, 423-442 4 Hwang, J. H., Kim, Y. J., and Kim, S. Y. Nonlinear hydrodynamic forces due to two-dimensional forced oscillation, Proceedings, I U T A M Symposium on Nonlinear Water Waves (Tokyo), SpringerVerlag, 1988, 231-238 5 Cointe, R. Remarks on the numerical treatment of the intersection point between a rigid body and a free surface, Third International Workshop on Water Waves and Floating Bodies, Woods Hole, Mass, 1988
6 Cointe, R. Calcul des efforts hydrodynamiques sur un cylindre horizontal en th6orie potentielle non-lin~aire, Deuxi~mes Journ~es de rHydrodynamique, Nantes, 1989, 89-104 (in French) 7 Cointe, R., Jami, A., and Molin, B. Nonlinear impulsive problems, Second International Workshop on Water Waves and Floating Bodies, Bristol, 1987
8 Cointe, R., Molin, B., and Nays, P. Nonlinear and second-order transient waves in a rectangular tank, BOSS'88, 1988, Trondheim 9 Cointe, R. Quelques aspects de la simulation numrrique d'un canal fi houle, Th~se de Doctorat de l'Ecole Nationale des Ponts et Chauss~es. Paris, 1989 (in French) l0 Grisvard, P. Problrmes aux limites dans les polygones--mode d'emploi, Pr~publication du laboratoire de math~matiques de rUniversit~ de Nice, 1984, (in French) 11 Peregrine, D. H. Unpublished note, 1972 12 Lin, W. M. Nonlinear motion of the free surface near a moving body, Ph.D. Dissertation, MIT, Cambridge, Mass, 1984 13 Kravtchenko, J. Remarques sur le calcul des amplitudes de la houle linraire engendrre par un batteur, 5th Conf. Coastal Eng., Grenoble, France, 1954, 50-61 (in French) 14 Roberts, A. J. Transient free surface flows generated by a moving vertical plate, Q. J. Mech. Appl. Math., 1987, 40(1) 15 Baker, G. R., Meiron, D. I., and Orszag, A. Applications of a generalized vortex method to nonlinear free-surface flows, Third International Conference on Numerical Ship Hydrodynamics, Paris, 1981, 179-191 16 Mei, C. C., and Unltiata, U. Harmonic generation in shallow water waves, Waves on beaches, R. E. Meyer, Ed., Academic Press, 1972, 181-202 17 Chapalain, G. Etude hydrodynamique et srdimentaire des environnements littoraux dominrs par la houle, Th~se de Doctorat de l'Universit~ Joseph Fourier, Grenoble, 1988, (in French). 18 Faltinsen, O. M. A numerical method of sloshing in tanks with two-dimensional flow, Journal of Ship Research, 22, 193-202 19 Dold, J. W., and Peregrine, D. H. Water-wave modulation, Bristol University, Report AM-86-03, 1986 20 Ursell, F. Surface waves in the presence of a submerged circular cylinder, I and II, Proc. Camb. Soc., 1950, 46, 141 158 21 Ogilvie, T. F. First- and second-order forces on a cylinder submerged under a free surface, J. Fluid Mech., 1963, 16, 451-472 22 Vada, T. A numerical solution of the second-order wave-diffraction problem for a submerged cylinder of arbitrary shape, J. Fluid Mech., 1987, 174, 23-37 23 Chaplin, J. R. Nonlinear forces on a horizontal cylinder beneath waves, J. Fluid Mech., 1984, 147, 449-464 24 Grue, J., and Granlund, K. Impact of nonlinearity upon waves traveling over a submerged cylinder, Third International Workshop on Water Waves and Floating Bodies (Woods Hole), 1988 25 Stansby, P. K., and Slaouti, A. On non-linear wave interaction with cylindrical bodies: a vortex sheet approach, Appl. Ocean Res., 1984, 6, 108-115 26 Dommermuth, D. G. Numerical methods for solving nonlinear water-wave problems in the time domain, Ph.D. Dissertation, MIT, 1987 27 Papanikolaou, A., and Nowacki, H. Second-order theory of oscillating cylinders in a regular steep wave, Proceedings, 13th Symposium Naval Hydrodynamics, Tokyo, 1980 28 Kyozuka, Y. Experimental study of second-order forces acting on a cylindrical body in waves, Proceedings, 14th Symposium Naval Hydrodynamics, Ann Arbor, USA, 1982 29 Faltinsen, O. M., Numerical solutions of transient nonlinear free surface motion outside or inside moving bodies, 2nd International Conference on Numerical Ship Hydrodynamics, Berkeley, 1977, 347 357 30 Vinje, T., and Brevig, P. Nonlinear, two-dimensional ship motions, Norwegian Institute of Technology, Report R-112.81, 1981 31 Maskell, S. J., and Ursell, F., The transient motion of a floating body, J. Fluid Mech., 1970, 44, 303-313 32 Tulin, M. P., and Cointe, R. Steady and unsteady spilling breakers: theory, Proceedings, I U T A M Symposium on Nonlinear Water Waves, Tokyo, 1987, 159-165 33 Cointe, R. A theory of breakers and breaking waves, Ph.D. Dissertation, University of California, Santa Barbara, 1987
Engineerino A n a l y s i s with B o u n d a r y E l e m e n t s , 1990, Vol. 7, No. 4
177
1. I N T R O D U C T I O N The direct numerical simulation of unsteady two-dimensional potential free surface flows using a Mixed Eulerian-Lagrangian method (MEL) has received considerable attention since the pioneering work of LonguetHiggins and Cokelet 1. 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 2, simulation of breaking in a tank 3, forced heaving of a cylinde?, etc...). Compared to the initial application of the MEL to the study of steep periodic waves (e.g., Ref. 1), several new difficulties 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 accuracy 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 computed and quantitative validations of the numerical procedure are almost non-existent, even for weakly 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 developed for the purpose of simulation of a two-dimensional wave tank using potential flow theory. Particular 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
© 1990Computational Mechanics Publications
comparison easier, an option of the code allows the linear problem 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 pistontype wavemaker. These waves can then interact with submerged or free surface piercing cylinders in forced or free motion. Reflection from the wall opposite tQ the wavemaker is avoided by the use of a damping zone--see Fig. 1. The results of the classical first- and second-order diffraction radiation theories can be recovered using the simulation. 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.
Y
2D O u r n r t ' / ¢ o l
toov~
I g-1
l'onR
x
1
E Plston-WPe wavemaker
Tested body
"Beach" (absorbing zone]
Fig. 1. Sindbad---geometric definitions
Enoineerino Analysis with Boundary Elements, 1990, 1Iol. 7, No. 4
167
Numerical simulation of a wave channel: R. Cointe This work has been reported elsewhere while in progress 5'6'7'a. 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 Ref. 9.
3. N U M E R I C A L 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 difliculties 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 yaxis is upward vertical--see Fig. 1 for geometric definitions. The fluid is assumed to be incompressible and the flow irrotational so that the velocity field ~ is given by: = V~,
(1)
with Aq~ =
O.
(2)
The computation is performed in a bounded domain. Along rigid boundaries (2 ~ Fn(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 2 on the free surface (2 ~ Fn(t)) and the associated value of the potential, t~(2). This yields:*
D~= - y -
- ~ qb~ +
D2 O~ = ~bsg + q~"fi'
02,-p+c(t)
(3) (4)
where D is used to indicate a material derivative, g and h are vectors tangent and normal to the free surface, respectively, and ( is an arbitrary constant. 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 c(t). With an appropriate choice of the velocity potential, this function can be taken equal to zero. At a given instant t, if 4, and ~b, are known along the free surface Fd(t), then the right-hand sides of (3)-(4) can be evaluated. The fact that ~b is harmonic in the fluid
*For the sake of simplicity, we use units such that the acceleration of gravity, g, the specific mass of water, p, and the depth of the tank, h, are equal to 1.
168
domain allows the value of q~, along Fa(t) to be computed from the values of q~ along Fa(t ) and of ~b, along F,(t). The kinematical constraint A~b = 0, associated with the boundary condition on F,, permits therefore to express the free surface boundary conditions (3)-(4) as an evolution equation for (qS, 2). 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 ~b knowing 4~ along Fd(t) and ~b, along F,(t).We use the integral equation:
--O(P)dp(P) + ~
dp(Q)G.(P,Q) dsQ
d r d+I'n
=f
gp.(Q)G(P,Q) dsQ,
(5)
J r a+r'n
where P is a point on the boundary, G is the Green function, 0(P) the angle between two tangents of the boundary at P (equal to rr for a smooth curve) and s a curvilinear abscissa along F. Equation (5) is discretized using a standard collection method. The boundary of the domain is approximated by segments and 4~ and ~b, are assumed to vary linearly along each segment. This allows an analytical 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 F.(t)]--Dirichlet [along Fd(t)] boundary value problem. It is well known that the solution of such a "mixed" problem is singular at the intersection points, i.e., at points belonging both to F.(t) and Fd(t)--e.g. Ref. 10. Let us consider, as an example, the wavemaker problem. For a piston type wavemaker and an horizontal free surface (90 degrees intersection angle), the complex potential W = q~ + i qJ 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:
4o(2, t) = d?o(2) + tchl(x) + . . .
(6)
Not surprisingly, performing such an expansion also leads to a z log z behavior for ~bl (e.g. Refs. 11 and 12). 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 nonlinear reoime, i.e. here, for a small acceleration of the wavemaket (relative to that of gravity) s. 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 condition; it is given by: q~,, + ~ = 0.
Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
(7)
Numerical simulation o f a wave channel: R. Cointe
If regularity in time is assumed, the local behavior of the complex potential for this solution can be shown to be in z 2 log z. Not surprisingly, this singularity is the same as that appearing for the harmonic problem (for which q~ = tp exp(icot)) which was studied by Kravtchenko 13. It is much weaker than that of the problem discretized in time. The singularity can also be studied in a similar fashion for an arbitrary angle of intersection 0. Assuming that the complex potential is bounded near the intersection, the leading behavior of the complex potential can be shown to be in z ~/° or z ~/° log z if n/O is an integer.* A consequence is that, for an intersection angle smaller than re, ~b is continuous along the boundary F while ~bs and tk, are piecewise continuous but experience a finite jump at the intersections. The numerical treatment at the intersection has been devised to accommodate such a behavior, referred to as weakly singular. More details on the numerical implementation can be found in Refs. 5 and 9. These results only apply in the weakly nonlinear regime, i.e., when the classical linear solution yields a first reasonable approximation to the problem. Obviously, it would be interesting to study other regimes. The impulsive regime (very large acceleration of the wavemaker) has been studied by several authors. As expected, the local behavior of the classical linear solution cannot yield any information for the nonlinear problem in this casC 4. It seems that Peregrine's solution 1~ corresponds to an outer solution to which an inner solution should be matched. Of course, only this inner solution is relevant to the study of the local behavior. If work has been carried out to find an inner solution 7'12, it leads to serious difficulties and, to our knowledge, 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 9.
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 diffracted 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 nonlinear 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 piston-type wavemaker. A "beach" is used for the absorption of the waves that are generated in the tank--see Fig. 1. It is in fact a damping zone, similar to that used in Ref. 15. In this zone, the free surface boundary conditions are modified by adding a *The same singularity occurs at corners of rigid boundaries.
damping term. We, therefore, write:
oo ot
D~ Dt
-
Y -
-
)
1
~ c~ + ~ c~ -
v(x)(c~ -
~,)
- (q~fi + c~,h - v(x)(5c - fee),
(8) (9)
where the subscript e corresponds to the reference configuration for the fluid. The function v(x) is homogeneous 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 v(x) is equal to zero except in the damping zone (Xo < x < xl), is chosen continuous and continuously differentiable, and is "tuned" to a characteristic wave frequency o9 and to a characteristic wave number k: v(x) = ctco
X - Xo
,Xo < X < Xl = Xo +
.(10)
If the proper scales have been chosen, values of ~t and fl 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 numerical 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 energy 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 boundaries is computed by two different methods. Both of them use Bernoulli's equation, but they differ by the evaluation of the ~bt term. In one case, this term is evaluated by finite differences in time. In the other case, it is obtained by solving the boundary integral equation for ~br 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 nonlinear 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 generation or absorption mechanism used, etc... A first step has therefore been to make proper comparisons
Engineerino Analysis with Boundary Elements, 1990, Vol. 7, No. 4
169
Numerical simulation of a wave channel." R. Cointe with linear results for the transient problem in a tank of finite length. The main interest of the method being to be able to account for nonlinear phenomena, a proper validation of the nonlinear response should be made. For that purpose, a comparison is performed with approximate nonlinear theories, such as second-order theory or shallow water theory. A final test is provided by comparisons with experiments. In our opinion, however, these comparisons should only be made last if it is the accuracy of the numerical scheme that has to be evaluated. In any event, a comparison 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 important 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 modeling.
'q;
/ F~.2. Sloshing--free surface profiles
4. NUMERICAL RESULTS The purpose of this section is to show some numerical results that can appropriately be compared to other theories or to experiments. Numerical instabilities--encountered by many similar methods--may appear for waves of large steepness. In this case, a 5 points smoothing algorithm similar to that used in Ref. 1 is employed.
4.2. Harmonic motion of a piston-type wavemaker The wavemaker problem is interesting because it provides a simple configuration to study most of the difficulties associated 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: l(t) = 0
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:
al
l(t) = - 2 cos(wt)
t >_ 0.
(14)
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 Ref. 3; an alternative solution that seems to have better convergence properties has also been derived s. Since this problem involves a moving rigid boundary that pierces the free surface, it allows the numerical 4.0_SECOND-ORDER __INITIAL SLOPE: 0.10 .... INITIAL SLOPE: 0.,35 2.0" O.
Y - - Ylinear
(~d) 2
(13)
The length of the tank is d. Note that even though the motion is harmonic for t > 0, transient phenomena are expected.
where d is the length of the tank. We show on Fig. 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 Fig. 3 the perturbation elevation, y~ -
t<0
(12)
0.0. to the second-order wave elevation that has been computed quasi-analyticallya. For ( = 0.1, a good agreement is achieved, indicating that both the linear and secondorder component of the wave elevation are accurately computed. For ( = 0.35, breaking occurs and the agreement deteriorates, suggesting that higher order effects become important.
--2.0
0.0"0:2'0:4"0'.6 X
Fig. 3.
Sloshing--perturbation elevation
170 Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
0:8
1.0
Numerical simulation of a wave channel." R. Cointe Table 1. Convergencetable NT~ NX --,
5
10
• 1.0 20
40
80
160 oo oo oo 1.3 1.4 1.4 1.4
4
I00
oo
oo
oo
oo
8 16 32 64 •28.
61.3 63.7 64.0 64.0 64.0 64.0
26.2 29.3 29.9 29.9 29.9 29.9
10.6 11.5 13.7 14.4 14.4 14.4
oo 4.2 5.3 5.4 5.4 5.4
oo 1.8 2.3 2.5 2.5 2.5
_.STEADY STATE LINEAR SOLUTION _LINEAR SIMULATION WITH ABSORPTION ....LINEAR SIMULATION WITHOUT ABSORPTION
z o
0.5 •
W
0.0
~
":
:-
::.
.... -,,"
-0.5
treatment 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 co = n/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 az) as a function of the number of nodes per wave length, N X , and the number of time steps per wave period, NT*. 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 numerical 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-order.** 4.2.2. Linear solution in a tank o f infinite length In order to test the efficiency of the damping zone, we consider now a tank of infinite length. For a piston type wavemaker, the solution far from the wavemaker and for large times is a progressive wave of amplitude: sinh 2 k a = al k + sinh k cosh k"
(15)
where k is the wavenumber associated to co (092 = k tanh k). Note that in deep water (k ---, oo), a = a~. We perform the simulation for the same wavemaker motion as in 4.2.1. but in a tank of length 10 equipped at its end of an absorbing zone. We compare on Fig. 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 ~t = 1 and ~ = 1. It appears deafly that the damping *A uniform grid is used and, thanks to the symmetry, the bottom is not discretized. No smoothing is applied. However, using the 5 points smoothing procedure introduced in Ref. 1, even at each time-step, does not significantly alter the accuracy of the method. **This indicates that numerical errors are mainly due to the numerical 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.
!
i
i
i
i
i
!
TlgE Wave elevation with and without damping zone
~.4.
zone provides a simple and efficient way to avoid any reflection. In order to evaluate more precisely the efficiency of the damping zone, we show on Fig. 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 = a,) at the same point. If long-lasting low frequency oscillations appear, the relative error is only of a few percents and the absorption mechanism 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 shallow water swell. Mei and Unliiata 16 have explained how such a swell can experience very drastic nonlinear deformations when it propagates. Very recently, Chapalain 17 has performed experiments at the Institut de Mrcanique de Grenoble that neatly confirm this bi-harmonic resonant behavior. It appeared therefore interesting to try and reproduce them with the simulation. The only data for the numerical simulation are the geometry of the tank, the law of motion of the wavemaker and the friction coefficient fw used to model
0.1 iz o ix: tw hi 0.0 LIJ
0(.
-0.1 40.0 Fig. 5.
66.0
86.0 " 100.0 1 2 0 . o "rIME Damping zone---relative error vs. time
140.0
Engineering Analysis with Boundary Elements, 1990, Iiol. 7, No. 4
171
Numerical simulation o f a wave channel; R. Cointe
B ^t~
~A
.,'.^
V
'J
V
^
'J
.,~.e
@. ;st.@
I.
;~¢*cI| Oo
..-, A
'.
A
r, A A A A
+..'\. ,,
,,,,,,~
& - , ~ . @ 1i
.;~.¢ I @.
dr.;
.
:
S.
I@*
$,@°@
IS.
A A AA
@°
""~v',Ov
o~¢°@
IB..
I~..
I~.
I$*
i,.
A A AAAAAI
I
v' v v v '-' V V k./~
*1~°@ . ; Il~c.@ ~ . @ I |d.@ .
I~
II. ll.~.@
lh
II.
liB.
II.
I~.
llo
l~.
/I"A'vAvAvA : , A'AAJvAvA"
lo~.CI .~¢°@ * ~I,P..I
:
.
ll.
I$*
:
.
I~*
Ig.
I ¢ . 0 .i,
l.
~-°
.A./~,_AA A
/~
Im.
I'I.
liB.
A A A A A A A A
*+$1;.@ I'*0 O.
I*
II*
Ill
Io1
'~ v
V V
II,
III1,
II.
'I~*
IS.
ll~.
I@~.@
" " t~ "~" z~ @
~
V
~
V
OI
11
ill°
I1~.@
V ~, v -
I~,
ll;l,
/
I1~.
Itl*° •
I
S*
I
ii11
[ Ile
l.a.
,
,
A4
A A A A /~ A /~ A A A A /~ A V / V ~ 'V' V V v V -V 'V 'V' kJ V V
.
-~.e
~',J
,
* ' * ~. . . . ~". 0 ~ .
v
Io@ •
I
!
~
!
III11
II~I. ?[,'5~S 1~ S*
l~.
lrJ.
-I
lo,
C
5
:0
:5
20 1EK)S
~5
30
~5
LO
Fig. 6. Shallow water swell--measured (left) vs. computed (rioht) wave elevations at several points in the tank dissipationt. The experiment was performed in a 40 cm deep tank with a piston-type wavemaker the motion of which was given by (14) with a I = 15.9 cm and co = 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.
tin order to model dissipative effects in a way similar to that used by Chapalain for Boussinesq equations we substitute to Bernoulli equation:
We compare on Fig. 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 Fig. 7. It appears that the simulation is very efficient to model shallow water waves and their generation by a wavemaker.
4.3. Arbitrary motion of a piston-type wavemaker 1
+
~, + _- v 6 . ¢ ~ + y + vO = 0 2
with 4
~1o)
We took the value of f, used by Chapalain, f~ = 0.1. Note that this modeling of dissipative effects, already used in Ref. 18 for the study of sloshing, is similar to the modeling of the damping zone.
172
In deep water, nonlinear phenomena in the propagation of a regular swell are rather long to develop. Numerical methods with periodic boundary conditions 19, 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 breaking. Experiments based on this principle were performed at MIT and a comparison with a numerical
Engineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
Numerical simulation of a wave channel: R. Cointe 025
70.0
60.0
oooI
_SINDBAD __EXPERIMENTS
(CHAPALAIN)
-0.~
I
I
/
I
I
0
10
Eo
30
40
50
o 30.0
0.00
~ 20.0.
-0.~ t 0
10
m
50.0
oo
40.0
0 Z
10.0. 0.0
0.0
510
16.0
1,.4.0 26.0
DISTANCE FROM WAVEMAKER(M)
25.0
, 30, 40,
,
20
50
60
0.25
I
0.00
I
-
-0.~
~
I
I
I
I
I
0
lO
20
30
40
50
0
~o
20
0
1
20
60
Fi9. 7. Shallow water swell--measured vs. computed Fourier components simulation similar to ours made3--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 absorbing zone and a somewhat 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 CRAY 1 for the simulation performed at MIT). The results of our simulation are in good agreement with both experimental and numerical results obtained at MIT up to breaking--see Fig. 8. If overturning develops at the same time, our computation fails sooner than theirs (before the closing of the tube)**. This confirms that the numerical simulation can reproduce accurately nonlinear phenomena observed experimentaly.
4.4. Wave diffraction on a submerged cylinder The case of the wave diffraction on a submerged cylinder allows a first study of wave-structure interaction. This problem has been studied extensively in the past. In particular, Ursell 2° used linear theory and showed that, for a circular cylinder, there is no reflection. Ogilvie 2. extended Ursell's results and computed the second-order vertical drift force. Very recently, Vada 22 computed the second-order potential and calculated the diffraction loads and the diffracted waves to second-order. Chaplin 23 measured diffractions loads in the laboratory while Grue and Granlund 24 measured the diffracted waves. These experiments have partly confirmed the results of first- and second-order theories. They have also exhibited some important nonlinear phenomena not accounted 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 objectives: , to recover the results of first- and second-order theories in order to assess the numerical accuracy of the method;
t"\..
O.O0 -0.~
I
"I °2'i
I
I
30
I
t
I
I
40
50
i
I
i
30
40
50
---
0
i
10
i
N
I
3O
i
I
40
50
6O
Fi9. 8. Steep deep water waves--measured (left, dashed line) vs. computed (solid line, left: MIT, rioht : Sindbad) wave elevations at several points in the tank
, to compare with experimental results in order to determine 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 25'26. 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 m = 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 direct integration of the pressure. In order to compare the results with those of experiments 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: Fx - F~ ) + ~ F~ I cos(n~ot + 0<"))
**Our 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.
60
r3(o 2
F, _ F~O) + ~ F~,)cos(n~ot + 0<")),
F3(D2
(16)
n_>l
(171
n>l
Enoineerino Analysis with Boundary Elements, 1990, Vol. 7, No. 4
173
Numerical simulation of a wave channel." R. Cointe o.o
where the subscript x and y denote the horizontal and vertical components of the force, respectively. As in Ref. 23 we introduce the Keulegan-Carpenter number, Kc. For a linear deep-water wave, K~ is given by: rca Kc = - - exp(kyc).
(18)
r
Table 2 gives the values of F~ ) and F~r") for n = 0, 12, vs. K~ in the case just described (that corresponds to Chaplin case E23). Following Chaplin, we write: F~)= ~ Cx,,,K ~,
-I.O
fJ'
*EXPERIMENTS (CHAPLIN) xSINDBAD . ~ x _SECOND-ORDER (OGILVIE) • x
-
.--2.0 0 -¢ -3.0
-4.0 -5.0
(19)
....................
-2.00
-1.50
-1.00
m_>l
Cy,,,K~,
Z
F~" ) =
(20)
• ........ -0.50
t.CK¢)
0.00
0.50
1.00
Fig. 10. Diffraction loads--vertical drift force
m_>l
The classical inertia coefficients are equal to C=11 et
Cyxl. According to linear theory, these are the only non-zero coefficients and they are equal. Ogilvie21 calculated them; they go to 2 as the immersion depth goes to infinity. We give on Fig. 9 the horizontal and vertical inertia coefficients vs. K~. For a small value of the KeuleganCarpenter number, both experimental and numerical results go the value predicted by linear theory, 2.25. As K¢ increases, however, the sharp decrease of the inertia coefficient observed experimentally is not predicted by the simulation. As argued by Chaplin, this is very likely due to viscous effects (creation of a circulation around the cylinder). 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, second-order theory and the present simulation. It
Table 2. Diffraction loads-- Sindbad ka
Kc
F(°) =x
F(°) --y
0.05 0.07 0.10 0.12 0.14
0.50 0.75 1.00 1.25 1.47
0.00 0.00 -0.01 0.02 -0.03
0.04 0.10 0.16 0.24 0.30
F(I) =x
Fit) --y
F(2) --x
F(2) =y
1.07 1.58 2.07 2.52 2.88
1.07 1.58 2.06 2.51 2.87
0.07 0.15 0.24 0.33 0.41
0.07 0.15 0.25 0.34 0.41
should be stressed that recovering the results from second-order steady theory with the simulation is by no means obvious: 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 24 performed experiments related to incoming deepwater Stokes waves passing over a restrained submerged circular cylinder. For a small cylinder submergence, a strong local nonlinearity is introduced at the free surface above the cylinder and free higher order harmonic waves are generated. We are reminded of the observations of Grue and Granlund concerning the wavefield far away from the cylinder: , upstream of the cylinder: an incoming Stokes wavetrain. No reflected waves, even to higher order; downstream of the cylinder: shorter free second harmonic waves of considerable amplitude are riding on the transmitted Stokes wave. If these trends are well predicted by second-order theory 22, the quantitative agreement with experiments is rather disappointing. The amplitude of the second harmonic free wave, a2, 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 Fig. 12. o.o
3.00
2.50 g,
-1.0 on B ' ~ o ~1=~, =.x
2.00"
x
x
8
x
*EXPERIMENTS (CHAPLIN) xSINDBAD _SECOND-ORDER (VAOA)
.J
v¢ - -2.0
o
- 1.50'
=g,0
t'N
~. -3.o
c~ 1.00 • EXPERIMENTS (CHAPLIN)
0.50
-4.0
xSINDBAD
_UNEAR THEORY(OGILVlE) -5.0
0 . 0 0
.
0.00
Fig. 9.
174
.
.
.
,
0.50
.
.
.
.
~
1.00 Kc
.
.
.
.
i
1.50
Diffraction loads--inertia coefficients
.
.
.
.
2.00
. . . .
-2.00
,
. . . .
-1.50
,
. . . .
-1.00
,
. . . .
-0.50
,
. . . .
0.00
,
. . . .
0.50
1.00
In(Ke)
Fio. 11. Diffraction loads-- response at double frequency
Engineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
Numerical simulation of a wave channel." R. Cointe 3.0
0.06
.... SINDBAD:
0.05
2.0
_SECOND-ORDER (¥ADA)
UNEAR CALCULATION _SINDBAD: ak = 0.03 _SINDBAD: ok = 0.09
0.04 -
1.0-
e~ 0.03-
o
0.02-
/
...~.~°
o
/,t
0.0-
iq
0.01
0.00
0.00
0.05
0.10 ok
0.15
0.20
Fig. 12. Diffracted waves--amplitude o f second-order free wave These findings suggest that the range of validity of second-order theory is quite narrow in this case. Observing that this theory predicts amplitudes of the secondorder free wave as large as that of the incoming wave, this should not appear as totally unexpected. In order to see if nonlinear free surface effects--and not viscous effect --are, indeed, responsible for this deficiency, it appeared interesting 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: tit = a cos(kx - cot + 0) + a t cos 2(kx - cot + O) + ...,
(21)
where a t is the amplitude of the second-order locked wave; , for the diffracted wave: qa = az cos(kx - cot + 01) + a~ cos 2(kx - cot + 01) + a2 cos(4kx - 2cot + 02) + " ' ,
(22)
where a~ 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 necessary to take a rather fine grid. The wavenumber of this free wave is, indeed, four times larger than the wavenumber of the incoming wave. The wavenumber of the incoming wave is chosen equal to k = 3.42 (so that we are in deep 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 y~ = - 0.1755. The simulation was performed 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 Fig. 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 Fig. 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, indicating that the "saturation" is, indeed, a nonlinear
-1.0-2.0
0.0
I'.o
2'.0
3'.0
4'.0
5.0
X
Fi9. 13.
Diffracted waves--free surface profiles
free surface phenomenon not accounted for by secondorder theory. Grue and Granlund observed breaking for ak ~- 0.085, while we were able to perform the numerical simulation up to ak = 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 numerical computation fails for ak = 0.12 after 3.2 periods. The computation does not blow up because of the overturning of the crest, as would have been expected, but because of a concentration of particles just after the cylinder. Physically, 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 the forces very much. 4.5. Wave radiation by a free surface piercino 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 cylinder have been studied quite extensively using linear and second-order theory 27'2s. Fully nonlinear calculations have also been attempted for forced heaving of a circular cylinder by a few authors. Initially 29'3° only the starting phase was considered. Recently, Hwang et al 4. 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 co = 3.16 (so that k = 10). The cylinder center has for elevation: Yc = 0
t <_ 0
Yc = ac sin(cot)
(23) t >_ 0.
(24)
Because there is no wavemaker in this ease, absorbing beaches (with • = fl = 1) are located at each end of the tank. We used 200 nodes on the free surfaced. In order to avoid too small or too large segment sizes near the
Enoineerin9 Analysis with Boundary Elements, 1990, Vol. 7, No. 4
175
Numerical simulation of a ,,ave channel: R. Cointe 0.01
Y
hi 0fY
o
0.00 o oe.
-0.01
).0
2'.0 ' 4'.0
o:o
1().0' 12.0
"I3~E Fig. 14.
intersection, a regridding procedure similar to that introduced in Ref. 26 was used when a node came too close or too far from the intersection. We show on Fig. 14 the transient vertical force on the cylinder as a function of time for a~ = 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 Fig. 15. Once a steady state is reached, we write the force as: Fr + 2ry, - F~°) + --rF (t°) sin(tnt) -- --rF(lV)cos(ot) 0.5Irr2 acO92 + ~ F~") sin(n~t + 0(")).
(25)
n=2
The second term on the left-hand side is the linear hydrostatic contribution. The acceleration-phase and velocity-phase components of the force at the forcing frequency would correspond to the added mass and damping coefficients for the linear problem'ft. We give on Table 3 the amplitude of the different harmonics as a function of ~ = a Jr. The Fourier analysis is performed on the four last periods of the signal. For small values of e, the results from linear and second-order theory are recovered 27'28. As the amplitude of the motion increase, the added mass increases while the damping coefficient decreases. This behavior is in agreement with available experimental observations 2s and other numerical results 4. 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 particular in order to assess the range of validity of approximate theories. If these results for forced motions are very promising, more experimental data would be needed to make proper comparisons. Note that for e = 0.6 the numerical computation breaks down before a steady state is reached, apparently because the cylinder goes out of the water.
Table 3. Radiationloads--Sindbad
0,1 0,2 0.3 0.4 0.5
?tNote that for the nonlinear problem the distinction between addedmass and hydrostatic components is somehowarbitrary.
F~°)
F~'")
V~'b)
F(2) _,
F(3) -,
F(r4)
-0.02 -0.03 -0.05 -0.07 -0.10
0.61 0.62 0.63 0.65 0.66
0.39 0.38 0.36 0.35 0.33
0.06 0.12 0.18 0.23 0.30
0.02 0.04 0.08 0.14
0.02 0.05 0.09
4.5.2. Free motion The free motion of a surface piercing circular cylinder was calculated using linear theory in Ref. 31. A fully nonlinear computation was performed in Ref. 30 using periodic boundary conditions. A good agreement between these two theories was achieved. Here, we again performed the same calculation. The body motion is calculated in a way similar to that used in Ref. 30 (see Ref. 9 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 Fig. 16 the subsequent vertical displacement of the cylinder, Yc, for Yo = 0.9 r. The simulation is performed exactly in the same conditions as in the preceding section. The absence of any nonlinear behavior for the cylinder elevation is rather striking and contrasts with the forced motion case.
2.0 -~
o>-. 0 >.-
b
1.0
_SINDBAD(YO = 0.9 r) __LINEAR THEORY (URS~'LL)
0.0
- 1 . 0
176
Forced heaving--free surface profiles
Fig. 15.
Forced heaving--transient force
'. . . . 6.0 8.0 10.0 12.0 1 LO T Free heaving--vertical displacement
0.0
Fig. 16.
Engineering Analysis with Boundary Elements, 1990, Vol. 7, No. 4
. . . . . . 2.0 4.0
,
,
•
',
•
N u m e r i c a l simulation o f a wave channel: R. Cointe
5. C O N C L U S I O N S T h e M E L has been a p p l i e d to the n u m e r i c a l s i m u l a t i o n of a wave t a n k e q u i p p e d at one end with a wave m a k e r a n d at the o t h e r e n d with a d a m p i n g zone. This configura t i o n allows t w o - d i m e n s i o n a l w a v e - s t r u c t u r e i n t e r a c t i o n p r o b l e m s to be studied. Results have been p r e s e n t e d t h a t include a wide r a n g e of a p p l i c a t i o n s , such as wave g e n e r a t i o n a n d a b s o r p t i o n , wave diffraction a n d wave r a d i a t i o n b y s u b m e r g e d o r surface piercing b o d i e s in forced o r free m o t i o n . Results from a p p r o x i m a t e theories (linear, s e c o n d - o r d e r o r shallow w a t e r t h e o r y ) can be recovered, n o t o n l y for s h o r t transient e v o l u t i o n s b u t also for s t e a d y state responses. Some nonlinear phenomena observed experimentally a n d t h a t have n o t been otherwise e x p l a i n e d are also a c c o u n t e d for. This seems to indicate t h a t the s i m u l a t i o n can be used as a " s t a n d a r d " that allows the validity of a p p r o x i m a t e theories to be assessed. A p p l i c a t i o n s of the m e t h o d in t h a t d i r e c t i o n suggest themselves. I n p a r t i c u lar, the roll m o t i o n of barges is being presently studied. If the M E L p r o v i d e s a n efficient a n d versatile w a y to s t u d y t w o - d i m e n s i o n a l free surface p r o b l e m s , it c a n n o t a c c o u n t for viscous effects (except in a very c r u d e w a y for instance for the m o d e l l i n g of b o t t o m friction). This is p a r t i c u l a r l y a p r o b l e m for viscous effects o c c u r r i n g in the n e a r vicinity of the free surface, i.e., viscous effects associated to breaking. B r e a k i n g is a m a j o r l i m i t a t i o n for the s i m u l a t i o n because the c a l c u l a t i o n has to s t o p w h e n e v e r a local b r e a k i n g event occurs. It w o u l d therefore be m o s t useful to be able to simulate b r e a k i n g a n d a s s o c i a t e d dissipative effects, even in a crude way. S o m e h o p e exists 32'33 for spilling breakers, b u t there is an o b v i o u s need for m o r e theoretical a n d e x p e r i m e n t a l w o r k on the subject.
6. A C K N O W L E D G E M E N T S This w o r k is a result of research s p o n s o r e d in p a r t by the " M i n i s t r r e de la Drfense," D R E T , u n d e r c o n v e n t i o n n u m b e r 88/073. This s u p p o r t is gratefully a c k n o w l e d g e d .
REFERENCES 1 Longuet-Higgins, M. S., and Cokelet, E. D. The deformation of steep surface waves on water. 1. A numerical method of computation, Proc. R. Soc. London, 1976, A 364, 1-26 2 Greenhow, M. Wedge entry into initially calm water, Applied Ocean Research, 1987 9, 214-223 3 Dommermuth, D. G., Yue, D. K. P., Rapp, R. J., Chan, E. S. and Melville, W. K. Deep-water plunging breakers: a comparison between potential theory and experiments, J. Fluid Mech, 1988, 189, 423-442 4 Hwang, J. H., Kim, Y. J., and Kim, S. Y. Nonlinear hydrodynamic forces due to two-dimensional forced oscillation, Proceedings, I U T A M Symposium on Nonlinear Water Waves (Tokyo), SpringerVerlag, 1988, 231-238 5 Cointe, R. Remarks on the numerical treatment of the intersection point between a rigid body and a free surface, Third International Workshop on Water Waves and Floating Bodies, Woods Hole, Mass, 1988
6 Cointe, R. Calcul des efforts hydrodynamiques sur un cylindre horizontal en th6orie potentielle non-lin~aire, Deuxi~mes Journ~es de rHydrodynamique, Nantes, 1989, 89-104 (in French) 7 Cointe, R., Jami, A., and Molin, B. Nonlinear impulsive problems, Second International Workshop on Water Waves and Floating Bodies, Bristol, 1987
8 Cointe, R., Molin, B., and Nays, P. Nonlinear and second-order transient waves in a rectangular tank, BOSS'88, 1988, Trondheim 9 Cointe, R. Quelques aspects de la simulation numrrique d'un canal fi houle, Th~se de Doctorat de l'Ecole Nationale des Ponts et Chauss~es. Paris, 1989 (in French) l0 Grisvard, P. Problrmes aux limites dans les polygones--mode d'emploi, Pr~publication du laboratoire de math~matiques de rUniversit~ de Nice, 1984, (in French) 11 Peregrine, D. H. Unpublished note, 1972 12 Lin, W. M. Nonlinear motion of the free surface near a moving body, Ph.D. Dissertation, MIT, Cambridge, Mass, 1984 13 Kravtchenko, J. Remarques sur le calcul des amplitudes de la houle linraire engendrre par un batteur, 5th Conf. Coastal Eng., Grenoble, France, 1954, 50-61 (in French) 14 Roberts, A. J. Transient free surface flows generated by a moving vertical plate, Q. J. Mech. Appl. Math., 1987, 40(1) 15 Baker, G. R., Meiron, D. I., and Orszag, A. Applications of a generalized vortex method to nonlinear free-surface flows, Third International Conference on Numerical Ship Hydrodynamics, Paris, 1981, 179-191 16 Mei, C. C., and Unltiata, U. Harmonic generation in shallow water waves, Waves on beaches, R. E. Meyer, Ed., Academic Press, 1972, 181-202 17 Chapalain, G. Etude hydrodynamique et srdimentaire des environnements littoraux dominrs par la houle, Th~se de Doctorat de l'Universit~ Joseph Fourier, Grenoble, 1988, (in French). 18 Faltinsen, O. M. A numerical method of sloshing in tanks with two-dimensional flow, Journal of Ship Research, 22, 193-202 19 Dold, J. W., and Peregrine, D. H. Water-wave modulation, Bristol University, Report AM-86-03, 1986 20 Ursell, F. Surface waves in the presence of a submerged circular cylinder, I and II, Proc. Camb. Soc., 1950, 46, 141 158 21 Ogilvie, T. F. First- and second-order forces on a cylinder submerged under a free surface, J. Fluid Mech., 1963, 16, 451-472 22 Vada, T. A numerical solution of the second-order wave-diffraction problem for a submerged cylinder of arbitrary shape, J. Fluid Mech., 1987, 174, 23-37 23 Chaplin, J. R. Nonlinear forces on a horizontal cylinder beneath waves, J. Fluid Mech., 1984, 147, 449-464 24 Grue, J., and Granlund, K. Impact of nonlinearity upon waves traveling over a submerged cylinder, Third International Workshop on Water Waves and Floating Bodies (Woods Hole), 1988 25 Stansby, P. K., and Slaouti, A. On non-linear wave interaction with cylindrical bodies: a vortex sheet approach, Appl. Ocean Res., 1984, 6, 108-115 26 Dommermuth, D. G. Numerical methods for solving nonlinear water-wave problems in the time domain, Ph.D. Dissertation, MIT, 1987 27 Papanikolaou, A., and Nowacki, H. Second-order theory of oscillating cylinders in a regular steep wave, Proceedings, 13th Symposium Naval Hydrodynamics, Tokyo, 1980 28 Kyozuka, Y. Experimental study of second-order forces acting on a cylindrical body in waves, Proceedings, 14th Symposium Naval Hydrodynamics, Ann Arbor, USA, 1982 29 Faltinsen, O. M., Numerical solutions of transient nonlinear free surface motion outside or inside moving bodies, 2nd International Conference on Numerical Ship Hydrodynamics, Berkeley, 1977, 347 357 30 Vinje, T., and Brevig, P. Nonlinear, two-dimensional ship motions, Norwegian Institute of Technology, Report R-112.81, 1981 31 Maskell, S. J., and Ursell, F., The transient motion of a floating body, J. Fluid Mech., 1970, 44, 303-313 32 Tulin, M. P., and Cointe, R. Steady and unsteady spilling breakers: theory, Proceedings, I U T A M Symposium on Nonlinear Water Waves, Tokyo, 1987, 159-165 33 Cointe, R. A theory of breakers and breaking waves, Ph.D. Dissertation, University of California, Santa Barbara, 1987
Engineerino A n a l y s i s with B o u n d a r y E l e m e n t s , 1990, Vol. 7, No. 4
177