Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Numerical Computations for a Nonlinear Free Surface Flow Problem K. J. Bai and J. W. Kim Seoul National University Seoul, Korea Y. H. Kim Dae Woo Shipbulding and Heavy Machinery Co. Jangseungpo, Korea ABSTRACT A finite-element method is presented for solving a nonlin- ear free surface flow problem for a ship moving in a towing tank. This problem is formulated as an initial/boundary value problem governed by Laplace's equation. The varia- tional functional used in the present paper is basically sim- ilar to the variational functional given by Luke (1967~. In the numerical procedure, the variational problem is reduced to a set of nonlinear ordinary differential equation where the unknown functions appear as the first deriva- tives with respect to time. In this numerical procedure, the governing equation, i.e., Laplace's equation, is treated as a constraint to the ordinary differential equation. In solving this problem, we employed the lumping and the unwinding schemes to maintain the numerical stability in the integra- tion with respect to time. To illustrate this method, a simple wedge-shape ship stretched to the bottom is treated. Com- putations are made for the range of depth Froude number, Fh = 0.7 ~ 1.1. The computed results show a good agree- ment with the previous results obtained by other methods. Even though the present computations are made mainly in the neighborhood of Fh = 1, this method can be applied to general Froude number. 1 Introiluction A free surface flow of an inviscid, incompressible fluid past a ship moving with a constant velocity in a towing tank is described by an initial/ boundary value problem governed by Laplace's equation with a free surface as a part of solu- tion. III the past, the problems of this type were generally treated after the boundary condition on the free surface was linearized. Recently, however, there are growing interests in treating the nonlinear free surface flow problems by various 403 numerical schemes. For example, two-dimensional problem was treated by Washizu et al. (1977), Hess (1977), Korving & Hermans (1977), Salvesen & van Kerczek (1978), Yen, Lee & Akai (1977) and Mori & Shin (1988~. Recently, a three dimensional nonlinear problem is treated by Dommermuth & Yue (1988~. There is another line of investigations on the nonlinear free surface flow problem based on the shallow water approx- imation which results in the Korteweg-de Vries, Boussinesq and Kadomtsev-Petviashivili equations and Green-Naghdi formulation. Many references along this approach can be found in Choi & Mei (1989~. To name few, Choi & Mei (1989), Ertekin, Webster & Wehausen (1984,1986) treated three-dimensional problems. In the present paper, we initially planned to investigate two different computational methods, i.e., the finite-element method and the fundamental singularity-distribution method for a three dimensional nonlinear problem. However, the progress along the second approach is not quite satisfactory at this time. Thus we describe only the application of a finite element method in this paper. One may find the ap- plication of fundamental source distribution method for a sloshing problem in three dimensional rectangular tank in Kim (1989~. To illustrate the present numerical method, computations are made for a wedge-shape ship extended from the bottom to the free surface. Main computational ef- forts are made for the critical Froude number, i.e., Fh = 1.0, just because of less computation time compared to the small depth Froude numbers. To treat the case of small Froude number, one simply has to take more finite elements. The computed results show two dimensional solitons periodically generated in the upstream and complicated three dimen- sional waves in the downstream. The comparisons show a good agreement with the previous results.
and then with respect to by, ~ J{o ~ I/SF /~ID V47 V§g} dV ~ = | dt ~ J~J. (O ~ And) 8+ dS (19) + fff v2¢ b~ ~ ~ . Here, [J = [Js + [J' Equation (18) shows that the stationary condition on J for the variation of ~ recovers the dynamic free surface boundary condition in each time and that the wave elevation at t = 0, t* should be specified as the constraints. Equation (19) shows that the stationary condition on J for the variation of ~ recovers the kinematic condition on SF and the governing equation. The variational form written above is previously given by Miles (1977) and is slightly different from that given by Luke(1967~. In the present variational formulation the wave elevation ~ is assumed to be known at t = 0, t*, whereas Luke assumed the potential ~ to be known at both initial and final times additionally. Specifically, the present varia- tional form is obtained from Luke's form by subtracting the volume integral of the potential resulted in the process of the integration by parts with respect to time. The present vari- ational functional has more advantage over the original Luke variational functional in treating the nonlinear free surface boundary conditions. 3.2 Finit - Element Discretization. The original initial/boundary value problem is well de- fined for fib of which the admissible solution should be twice continuously differentiable in space. However, in the above variational method, it suffices that the admissible trial func- tions have the square integrable properties of the function ~ and only their first derivatives in space. This enables us to look for an approximate solution in a wider class in the variational method. As in usual finite element analysis we discretize the fluid domain into finite number of finite elements. In the present computations, the finite elements are generated such that the projections on the horizontal plane, i.e. z=0, is un- changed while the other coordinate, i.e. the z-axis, is allowed to change in time. This simplifies the mesh generation and the computations considerably. However, this restriction is not necessary for the present method in applying to a gen- eral problem. Then we approximate ~ in N dimensional function space whose basis is continuous in D and has con- tinuous derivatives in each element. However, it is allowed to have a finite discontinuity in the normal derivative across the common boundaries of the adjacent elements. We de- note the basis of this space by {Ni~i=i,...,N and approximate by the span of the restrictions of {Ni~i=~,...,N on SF which is also continuous and piecewise differentiable on SF; ¢(x, y, z, t) = Ditty Ni(z, y, 3; S) (20) <(x,y,t) = Skits Mk~x,y) (21) where Mi~x,y) = ~i`~,Y,8;~z=: k= 1,--.,NF (22) and NF is the number of nodal points on SF and i' is the nodal number of the basis function Ni of which the node coincides with that of the free surface node k. Summa- tion conventions for the repeated indices are used here. It should be noted that the basis function {Ni~i=~,...,N is de- pendent on the free surface shape z = §(x,y,t) but its re- stricition on SF is the function of (x,y) and independent of (. This special property of {Mk~k=~, ,NF is due to the simplicity of the fluid domain D which has no variation in z- direction such that we can change the position of the nodal points only in the z-direction in each time step. Once the trial function is approximated by using the above basis function, we obtain the Lagrangian L for this trial solution as L = ~i,~Til~ (231) - 2¢iKij~j- 2fkPHn where Tkl = ||- MkM' dS SF Kij = ||| VNi VNj dV (24) Pal = |,/- MkMI dS. SF The tensors Kij, Pal are the kinetic and potential energy tensor and Tk' is a tensors obtained from the free surface integral, which can be interpreted as a tensor related to the transfer rate between these two energies. It is of in- terest to note that in Eq. (24), Tk' = Phi. However Tkl will be defined differently from this in the computation through lumping. 405
The stationary condition on J = f Ldt is equivalent to the following Euler-Lagrange equation To ~ = Kim ~i, (25) T,c, ~i' = - 2¢i phi JOPit 0, (26) fork=l, ,NF Kij Hi = 0, i ~ it. (27) Here Eq. (25) and Eq. (26) are the nonlinear ordinary differential equation for Id, Ike=, ,NF and Eq. (27) is the algebraic equations for {¢iJi~i,, which is the constraints for the above two equations. 3.3 Numerical Dispersion Relation. If one solves the above variational formulation numeri- cally, one encounters the effect of numerical dispersion de- pending on the specific numerical schemes employed. How- ever, it is difficult to analyze the numerical dispersion erect in the above nonlinear three-dimensional problem. There- fore, we restrict ourselves to a linear problem and test the numerical dispersion in the following. The linear version of the Eq.'s (25~~~27) is given as Tang = K° j~j (32) K,°j Hi = 0, i ~ it (33) It can be easily shown that the solution of the above Tkl Ail = -Pal ~ (34) discretized problem satisfies the conservation of mass and total energy, i.e. (~ Pap) = 0 (28) dt k,] dd(~qkKijjj+~SkP.l`~) = 0 (29) s,3 k,! This property of the conservations is independent of the ten- sor Tk' if it satisfies ~ Ti' = ~ Pa (30) k k It should be also pointed out that the direct use of (26) leads to some difficulty in the computations. This difficulty arises from the first term in the right-hand side which is the derivative of the kinetic energy tensor with respect to the wave elevation. We have avoided this difficulty by utilizing the fact that (26) is equivalent to the condition of vanishing of the right-hand side in Eq. (18~. In Eq. (18) b; can be regarded as test functions on SF. Then Eq. (26) can be given as Talkie = ||- Mk¢Z§t dS SF - ~ ||- M~IV¢l dS Peal . If the integrals in Eq. (26) and (31) are evalutated ex- actly, they are equivalent. However, in the present compu- tation these integrals are calculated by integral quadrature rules. Therefore, the conservation of energy may not be satisfied exactly due to the error caused by the numerical integration. This test result is given in the next section. for k = 1, · · ·, NF where K,°j is the linearized kinetic energy tensor evaluated in the undisturbed fluid domain. We use the 8-node brick element with linear interpolation along each coordinate. Then by separation of variables the velocity potential ¢(x, y, a, t) can be represented as a prod- uct of the following three sets of functions: (x, t) = Hi (t) Xi~x), i = 1, · · ·, No, (y,t) = ~~.(t)Yj~y), j = 1,···, Ny, ~z,t) = ~k~t)X~(z), k = 1,··., Nz (35) (36) (37) where Xi~x), Yelp), Z~(z) are the one dimensional basis func- tions and Nz, Ny' Nz are numbers of nodal points in x,y and z coordinates, respectively. Taking Fourier transformation from the time domain to the frequency domain, w, we have the following set of eigen-value problems. .~(Xi,X'kZXiXj~dx~jZ = 0 (38) t(YiYJ - k2YiYj~dyy>1. = 0 (39) to <~2 bail - J Size + k2ZiZj~dz ~' = 0 ' , where ~' = LEO) and k2 = kz + k2. Here bit = 1 if i = 1 and zero otherwise. It is understood that the time depen- dent terms are now functions of frequencies without chang- (31) ing their notations. Moreover, we can treat a two-dimensional wave prop- agating in the x direction without loss of generality since the linear wave is isotropic in OX plane as long as we use the same mesh size in x and y directions. To obtain the dispersion relation, we assume the type of solution for By as ~ (x, w) = Ae (41) where A is a constant amplitude and 406
An = , n^z n = 17 2, 3, · · · (42) is the discrete wave numbers and /`x is the horizontal mesh size. By solving Eq. (38) and (40) after substituting Eq. (41) into Eq. (38), the dispersion relation can be obtained as where ~15.~ en ~ 10. k2 - 6 1cos(8n~x) <43~ 5 n '` =2 2 + COS(8n^X) O) = W (en) pN-lGlF2 _ eN-lG2F §2 F1,CNF2 (44) 2 Em/ Pi = 1 +~U /3 - ,Bi(1 - 6 A, (46) Gi = 1-~2/6-,Bi(l+ 3 ), i= li2 (47) with the number of elements in z-direction N = Nz - 1 and = k"Az with vertical mesh size Liz = 1/N. The well-known exact linear dispersion relation corre- sponding to Eq. (43), (44) is k2zact = (9n (48) (~)c tact = k tanh k. (49) The comparisons between these numerical and exact dis- persion relaions are shown in Fig. 1. It can be shown that for Ax << 1 with en held fixed kn = k2ZaC`(l+o(Ax2~), (50) which means that W(8n) has an error of of, however small values of fly we choose. Go .......................... 0.5 ~2 <1 k 0.0 1 O. ~ - -- -- --- -' // 2 ~ / /N '''''''' ''''' '''''''''' ' ''''' ' go'''''''''',? '' ~/~ '' Exact ...~_............................................................................... k 6. 8. 10. Fig. 1 Comparison of numerical and exact dispersion relations. (a): Eq. (43) and Eq. (48) (b): Eq. (44) and Eq. (49) N: No. of elements in vertical direction. Fig. 1 shows that the numerical dispersion relation pre- dicts higher values for co than the exact values. This inher- ent property of finite element method based on variational formulation restricts the numerical stability in time-domain analysis. A simple remedy to correct this inherent draw- back is the use of the so-called 'lumping' in the tensor To which is originally equal to Pal as in Eq. (24~. The lumped tensor Tk' is given as TO = Ail ~ Pam (51) m where Ike is Kronecker's delta. Then the dispersion relation = W ((7n) (52) 2 + COS(8n ~ 27) W `~9 ~ . . , j ~ ;, F.E.M/ / Exact .' / . ~',' ...................... 0.0 0.2 0.4 0.6 0.8 Onyx Or (a) The order of correction term in Eq. (52) is of) which is same as the error in W(8n). For the shortest wave length 2~:z:, W(8n) of the lumped case predict one third of the values of W(8n) of the unlumped case. This means that with lumping we can take time step three times larger than the standard finite element method without loss of the or- der of accuracy in the dispersion relation; The dispersion correction given in Eq. (51) by the lumping also preserves the conservation of the mass and energy, since Eq. (30) is satisfied when Tk' is replaced by Tkl defined in Eq. (51~. 3.4 Treatment of the Convective Term. We assume that a ship is moving with a Etroude number Fh. Then the pure convection terms are added to free sur- face boundary conditions. Since the presence of convection term due to the introduction of a moving coordinate does 407
not change the original dispersion and nonlinearity, we ana- lyze them separately and then add it to the original problem. Therefore if suffices to consider a pure convection problem separately here. Let a convection problem for ~ be given as ~ =Fang on SF. (53) The numerical scheme based on a direct weak formulation for the above equation is known to have a considerable phase error for short waves. To reduce the contamination of the phase errors in the computational domain caused by the short wave components, an unwinding scheme by using asym- metric test function has been successfully employed as in Hughes & Brooks (1982~. The discretized form of Eq. (53) based on this scheme is obtained as T,dn =Fhii (Mi+ 2 a Ma) a Mi dS (54) where ~ is the unwinding parameter and has a value between O and 1. If c' = 0, it is equivalent to the standard Galerkin method based on symmetric test functions. If c' = 1, the right-hand side of equation (54) is equivalent to the forward- difference formula in finite difference method. ~ .o I <' k 3 0.5 3 ............. ...................................... . ~ 0.0 '/ ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''~''''''''''''''''''''''''''''''''''''~''''''''''''''~ 0.0 0.2 0.4 0.6 And 2: 6. 4. 2. Exact (a) Real part 0.8 1.0 / . ~ ..... .......... .... ............ ....... ....... . . 0.4 0.6 0.8 1.0 6'n~ Z: Or (b) Imaginary part Fig. 2 Numerical dispersion relation due to the con- vective term for Fh = 1. (Eq. (55~) In a similar procedure given in Eq. (41) through Eq. (43), the following numerical dispersion relation is obtained for Fh = 1: +i 3 sink - Ax 2 + cos(8n~x) 3a 1cost) . /`x 2 + cost) (55) The real and imaginary parts of the above equations are plotted in Fig. 2a and 2b, respectively. In Fig. 2a the real part of ~ is compared with the exact relation ~ = 6'n. (56) In Eq. (55) it can be seen that the upwinding procedure has no effect on the real part of w. The additonal term in the imaginary part due to the unwinding behaves as a diffusion term and plays a role of damping mainly to the short wave component in time. The computations based on Eq. (54) are made for several values of ~ to test the effect of the numerical damping. The result is given in the next section. 3.5 Treatment of the Radiation Condition. In the present computations, the computation domain is taken sufficiently large in the upstream direction so that the generated solitons travelling upstream does not hit the up- stream boundary at the final time of computations. How- ever, along the downstream radiation boundary, we used a simple boundary condition in a moving frame of reference as ¢'n = 0 at x = XR (57) where x = XR iS the radiation boundary at downstrea~rn. This can be interpreted as the disturbance convected away with the negative velocity of the moving coordinate. The present computed flow cases are such that the Froude num- ber is around one. Accordingly, the velocity of the moving frame of reference is the critical speed. If there are any reflected waves from the downstream radiation boundary, then the speed of the reflected waves will not be larger than the critical speed. Therefore one may expect the reflected wave cannot contaminate the solution in the computation domain. Strictly speaking this is true only for a linear prom fem. However, even in nonlinear case the speed of waves at downstrem do not exceed the critical speed in most cases. This simple numerical treatment of the downstream radia- tion boundary condition is tested successfully and the result is given in the next section. 3.6 Time Integration on the Free Surface. Up to the previous subsections, we presented the numer- ical treatments of the governing equations, the convective terms and the radiation conditions as boundary value prob- lem. Once we discretize the computation domain by the 408
finite elements and obtain the discretized form of matrix equation after integrating out all the space variables, then the problem is reduced to a set of nonlinear ordinary dif- ferential equations to solve. In this reduced equations, the time derivatives are present in the potential and the free sur- face elevation, both defined on the free surface. The reduced set of ordinary differential equations becomes the following form: Tklp = FhCt ~ (58) +Ki~j~j + fin, T`,` bib = FhC~j~jPad (59) Hi Dsk'~ 2 Hi D2kj 45j for k= 1,-~.,NF, where Kiwi = -f', i ~ it (60) Hi, = ifs (Mi + 2 azMi)azM' dS Ckj = ii (Me + 2 ~ Mi)3 Nj dS Dsk' = ifs MiMia3Ni dS (61) Dskj = /! M6VNi · VNj dS fi = Fhii Nine dS. so It should be pointed out that Eq. (60) can be interpreted as a constraint to Eq. (58) and (59~. Equation (60) is obtained from the boundary value problem for Laplace equa- tion with an essential condition on free surface and a natural condition on the body surface. In the solution procedure for this problem, the constraint, i.e., Eq. (60) is first solved by the Jacobi conjugate gradient method. Then the matrix Tkl is inverted. Here the other matrices, which are dependent on the free surface shape, are treated as known from the previous time step. The fi- nal form in Eq. (58) and (59) is solved by the fourth-order Runge Kutta method. 1 2b ~ ~ ~ ~ 2B -_ L :\\\\\\W y = -B A\\\\\\ , Fig. 3 Sketch of ship model and tank in horizon- tal plane view. Here the z-axis is against the gravity. 4 Numerical Results and Discussions The computed model as shown in Fig. 3, is a vertical- sided wedge-shaped ship extended from the free surface to the bottom of the tank. The length of the wedge ship, the length of the parallel middle body, and the beam are denoted by L, Em, and 2b, respectively. The tank has the width of 2B and the mean water-depth of h. The ship is assumed to move with a constant velocity ofU along the centerline of the tank. In presenting our computed results,all the phys- ical quantities are nondimensionalized by h, ph3 and for the length, mass, and time, respectively as mentioned previously. Before we made the computations for the ship in the tank, we tested the case of a two-dimensional free oscillations in a three dimensional rectangular tank as an intial-value problem. The results are given in Fig. 4. Here the motion started from an initial hump on the free surface and let the time increases. The normalized tank length is 40 and the final time in nondimensional scale is 60 with the time step of 0.2 . Fig. 4a and 4b show the results of 15 wave elevations with time increment of 4. The computed wave elevation and the speed of propagation agreed well with the approximate theoretical results for the amplitude and speed of the wave, which is given in Appendix A. ~ . o ._ 0.25 0.00 o ._ - n an t=0 t=36 / V''.'''''..'.'''''''.''''''.'''''''''''''''''''''''''''''''''.'''''.''''''''''''''''''.'''''.''''''''.''''''''''''.'''.''''.'''''''''.~''.''''''.'..'''''''''''.'..'.'''.''''''''' 0- 10. (a) 20. 30. 40. X 0.75 0.50 0.25 0.00, -0.25 t=60 · t t=40 1 ~ 0. 10. 20. X (b) 30. 40. Fig. 4 Numerical test of free oscillation in a tank started from an initial hump at t=O. Tank length = 40, Ax = 0.4 and At = 0.2 409
a, 0.25 ct 0.2= 0.50 o ,,., 0.25 ct 0.00 . (a) t 0. 10. 20. -0.25 20. - 10. 0. x (b) <t = 0.01 10. 20. Fig. 5 shows the test results of the conservation of the total energy, E, which is sum of the kinetic energy K and the potential energy P. The computed case is same as Fig. 4. Here the energy is normalized by the initial energy. This result shows that the total energy is conserved fairly well throughout computations. Lo On l . ~ O. 10. 20. 30. 40. 50. 60. Fig. 5 Numerical test of energy conservation for the case of Fig. 4. K: kinetic energy P: potential energy E = K ~ P: total energy n an- 0.50 o .,. 0.25 Ct C) -0.25 20. - 10. 0. X (e) ct = 0.1 o ·.-, 0.25 Cal - 0.00 0.25 0.50 1 20. - 10. 0. X (C) Cl! = 0.03 0.50 o a, 0.25 - 0.00 -0.25 20. ~ t 10. 20. 10. 0. X (d) a = 0.05 10. 20. 10. 20. Fig. 6 Numerical test of the effect of unwinding pa- rameter c2 for the convective term across the radiation boundary. The initial disturbance is saline as in Fig. 4. (a): Computed in the inertial coordinates. (b)~(e): Computed in the coordinates moving with velocity -1. Fig. 6 shows the numerical test of the effect of the up- winding parameter, c', for the convective term across the radiation boundary. The initial disturbance is same as in Fig. 4. I?ig.6a shows the wave profiles first computed from 410
the formulation defined in the . , . at, , . ~ . inertial coordinates and the results are shifted to`,/ght to compare those given in Fig. 6b through Fig. Be. Fig. 6b through Fig. Be are obtained from the results computed in the coordinates moving with the velocity of-1 for four values of the unwinding parame- ter, cz = 0.01, 0.03, 0.05, 0.1. These comparisons show the dual roles of upwinding on the numerical damping and sta- bility. This means that too much unwinding would give too much numerical damping effect in the numerical results. In the moderate range of the values of a, the unwinding works favorably to the stability. Fig. 7 through Fig.12 are the results for the following computation conditions: Wedge model: Computation Domain: Finite Element Meshes: Upwinding Parameter: Mesh sizes: L=8, b=0.4, Lm=4 x=-30, 30), y=~-4, 0) (Symmetry is used) 150 x 10 x 1 elements c'= 0.05 Ax=0.4, /`y=0.4, /`t=0.2 Fig. 7 throughout Fig. 10 shows the.results of wave profiles as time increases with the time increment of 20, for the cases of Fh =1.0, 0.9, 0.8, 0.7, respectively. The gen- eration of the upstream solitons are pronounced as the Fh approaches to 1. For Fh = 0.7, the generation of upstream solitons is barely noticeable. Fig. 11 shows the results of supercritical Froude number, i.e., Fh=1.Q5. It is of interest to note that the numerical instability shows up when t=60. Further computations for t larger than 60 showed that the water surface hit the bottom locally. At present, our pro- gram cannot incorporate with a local dry bottom. This can be easily incorporated in the near future. Fig. 12a and 12b show the computed wave resistance for Fh=O.9 and 1.0, respectively. Fig. 13a and 13b show the time record of wave elevations for Fh= 0.9 and 1.0, re- spectively. This shows that the time-dependent hydrostatic pressure is not negligible in the resistance computations. Fig. 14a and 14b show the computed results of wave resistance for two different tank conditions: Fig. 14a is for the case of the tank and ship geometry being reduced to one half in the direction of the y-direction. Fig. 14b is for the case of those geometry reduced to one-tenth in that direc- tion. Thus the wave resistances are shown by multiplied by the factor of two and ten, in Fig.14a and 14b, respectively. Table 1 shows the amplitude, speed and the generation period of solitons for the blockage coefficient, S = 0.1. Here the blockage coefficient, S is defined as in Ertekin et al. (1986~. The amplitude A is taken from the first solitons. The period is measured from the computed results at the F.P. Also shown are the computed results by Ertekin et al. who treated a pressure patch in Table 1. Table 2 shows the amplitude, speed and the generation period of solitons obtained in the tank reduced to one half of the tank width for several values of the blockage coefficients when Fh = 1. The amplitudes and the speed increases whereas the generation period decreases as the blockage co- e~cient increases. This trend agrees well with the previous experimental results of Ertekin. Table 1. Amplitude, Speed and generation period of solitons for the blockage coefficient, S=0.1. The values in the parentheses are obtained from Ertekin et al. (1986~. U/~ A/h C/~b UT,/h 0.7 0.139 1.06 14.6 0.8 0.240 1.10 19.4 0.9 0.375 1.16 24.1 (0.5101) (1.224) (20.0) 1.0 0.553 1.24 30.0 (0.6248) (1.280) (29.6) 1.05* 0.677 1.27 34.7 1.1* 0.806 1.32 39.1 (0.7729) (1.390) 39.3 * The values are obtained from the narrow tank reduced to the 1/10 in y-direction. Table 2. Amplitude, Speed and generation period of solitons obtained in the tank reduced to one half of the tank width for U/~iZ = 1. Block. coeff. A/h C/ UT§/h 0.06 0.08 0.1 0.12 0.404 0.475 0.545 0.604 1.18 1.22 1.24 1.26 44.2 36.2 29.8 27.0 As a concluding remark the accuracy of the present com- putations could be improved by employing finer meshes. The computation time for a typical model of 1500 elements was 260 seconds for each time step by IBM PC/XT with T800 Monoputer. A typical number of the total time steps was 500. ACKNOWLEDGEMENTS. This work has been supported by the Nonlinear Ship Hy- drodynamics Program supported by the Korean Science & Engineering Foundation. 411
REFERENCES Choi, H. S. & Mel, C. C. 1989 Wave resistance and squat of a slender ship moving near the critical speed in re- stricted water. Proc. 5th Int. Conf. on Num. Ship Nydrodynam., Hiroshima. Dommermuth, D. G. & Yue D. 1988 The non-linear three- dimensional waves generated by a moving surface dis- turbances. PTOC. 17th Symp. Naval. Hydrodynam., Hague. Ertekin, R. C. 1984 Soliton generated by moving distur- bavces in shallow water: theory, computation and ex- periment. Ph.D. di~eration, University of California, Berkeley. Ertekin, R. C., Webster, W. C. & Wehausen, J. V. 1984 Ship-generated solitons. Proc. 15th Symp. Naval Hy- drodynam., Hamburg, pp. 347-364. Ertekin, R. C., Webster, W. C. & Wehausen, J. V. 1986 Waves caused by a moving disturbance in a shallow channel of finite width. J. Fluid Mech., 169 pp. 347- 364. Green, A. E. & Naghdi, P. M. 1976 A derivation of equations for wave propagation in water of variable depth. J. Fluid Mech. 78 pp. 237-246. Hess, J. 1977 Progress in the calculation of nonlinear free- surface problems by surface-singularity techniques. Proc. 2nd Int. Conf. on Num. Ship Hydrodynam., Berkeley, pp. 278-284. Hughes, T. J. R. & Brooks, A. 1982 A theoretical frame- work for Petrov-Galerkin Methods with discontinuous weighting Functions: Application to the Streamline- upwind procedure Finite Elements in Fluids, vol. 4, John Wiley & Sons Ltd., pp. 47-65. Kim, Y. H. 1989 A numerical analysis of free surface wave problems by source distribution method. M.S. Thesis, Seoul Nat'1 Univ. (in Korean). Korving, C. & Hermans, A. J. 1977 The wave resistance for flow problems with a free-surface. Proc. 2nd Int. Conf. on Num. Ship Hydrodynam., Berkeley, pp. 285- 291. Luke, J. C. 1967 A variational principle for a fluid with a free surface. J. Fluid Mech. 27 pp. 395-397. 412 Miles, J. W. 1977 On Hamilton's principle for surface waves. J. Fluid Mech. 83 pp. 395-397. Mori, K. & Shin, M. 1988 Su~breaking wave: it's character- istics, appearing condition and numerical simulation. Proc. 17th Symp. Naval. Hydrodynam., Hague. Salvesen, N. & van Kerczek, C. 1978 Nonlinear Aspects of Subcritical Shallow- Water Flow Past Two-Dimensional Obstructions. J. Ship Res. vol. 22, No. 4, pp. 179- 200. Washizu, K., Nakayama, T. & Ikegawa, M. 1977 Applica- tion of the finite element method to some free surface fluid problems. Finite Elements in Water Resources, Pentech Press, London, pp. 4.247-4.226 Wu, T. Y. 1981 Long waves in ocean and coastal waters. J. Engng Mech. Div. ASCE 107 pp. 501-522. Wu, D. M. & Wu, T. Y. 1982 Three-dimensional nonlinear long waves due to moving surface pressure. Proc. 14th Symp. Naval Hydrodynam., Ann Arbor, Mich., pp. 103-129. Yen, S. M., Lee, K.D. & Akai, T.J. 1977 Finite-element and finite-difference solutions of nonlinear free-surface wave problems. Proc. 2nd Int. Conf. on Num. Ship Hydrodynam., Berkeley, pp. 305-318.
Fig. 7 Computed free surface for Fh = 413 ~-
Fig. 8 Computed free surface for Fh = 0 9 414
t60 Fig. 9 Computed free surface for Fh = 0.8 ~60 Fig. 10 Computed free surface for Fh = 0 7 415
0.6 ... 0.4 v t=6G Fig. 11 Computed free surface for Fh = 1.05 0.2 1- 0.0 0.6 ....................... 0.4 0.2 : 0.0 1 / ~ ~ / , ~ / , .: ..... : ..... ........ ..... ............................. ............ ..... ... ....... ~ . ~ 0.6 -/ /'''''' ..~ ~ . 1 _ . .~ ...... / .. .. : ~ ............................................ .................................. .............................. .. .... ....................... ......... ...................... ' ~ ~ . ~ ...................................... . ~ ' . ~ , . . ~ 0.2 ~ ' ~ ~ .... ....................... ..... .................. . .......... .... .. ... .. ............ .... .... ..... ........ .. .. ......... .......... ..... . .... ... 1 ' ~ . 0.0- 0. 20. 40. 60. 80. 100. 0. 20. 40. 60. 80. 100. t t (a) (a) ................ ........ .. ..... .. ... .. .............. ..... ....... ..... ..... . .... . .. .... .................. ..... ...... . ... ... ...... ... . . . ...... o in, 0.4 - o .-, 0.4 C: - ~ 0.2 Ct 0. 20. 40. 60. 80. 1 00. 0 - - (b) t Fig. 12 Wave resistance for different Froude num- bers (B = 4, b = 0.4). (a) Fh = 1 (b) Fh = 0.9 .. ... .... '/ \~/' W`~ it_ _ 0.6 2 f 0. 20. 40. 60. (b) t Fig. 13 Wave elevation at F.P of a wedge. (a) Fh = 1 (b) Fh = 0 9 416 80. 1 00.
~2 * an or 0-4 * ~ 0.2 0.4 0.2 / . ~ '. / .. - .................................. / : ......... ~ ............................. ......... ..... .... .. .. ... ....... . / ., ....... .... .......... . . . . . ....... ........ . ...... . .. 40. 60. o.o 1 0. 20. 0.6 ........................................................... 80. 100. (a) t ......................... ' ........................... ........................ . / ~ / o. GN, Wu : L1 0.01 ~ so. loo. L2 20. 40. 60. (b) t Fig. 14 Wave resistance for different tank and hull geometry (Fh = 13. (a) B = 2, b = 0.2 (b) B = 1, b = 0.1 A pp endix '~3 3- A Comparisons with Other Theories. The comparisons are made on the linear dispersion rela- tion and the speed of the nonlinear wave of permanent form among several different approximate theories. We present the results obtained by the equations derived by Green & Naghdi(1976)(hereafter the GNequation)andWu(1981~. Since both equations are derived with the assumption of shallow-water limit or equivalent assumption on the veloc- ity field, we compare them with the result of finite element method with one element in depth-wise direction. For the finite element method we present two different approxima- tion in depth-wise direction. In the first approximation a linear interpolation is used. The results presented in the section 4 were obtained by using this element. In the sec- ond, we use a quadratic interpolation which satisfies the boundary condition on the bottom. The above two approx- imations will be refered to L1 and L2, respectively, hereafter. The approximate velocity potentials of L1 and L2 schemes in two dimensions are given as L1: ¢' (x, z, t) = fit (z, t) + Z f2 (z, t) (A.1) L2: Adz, z, t) = fit (x, t) + (z + 1~2 f2 (z, t). (A.2) A.1 Linear Dispersion Relation The linear dispersion relations of GN and Wu equations are identical and can be found in Ertekin (1984~. The re- lations for L1 and L2 scheme can be derived from Eq. (40) with appropriate basis functions. The results are 2 3k2 3 + k2 2 k2~12 + k2) · ~ = 12 + 4k2 2 _ k2~15 + k2) 15 + 6k2 (A.3) (A.4) (A.5) The above relations are plotted in Fig. 15 compared with the exact relation. It can be found that the finite element method gives the upper bound for the exact values of A, whereas GN and Wu equation gives the lower bound. The L2 scheme predicts more closely to the exact values of whereas the L1 scheme deviates more than others. 6. 5. 4. 2.1........... 0.- O. ... .......... . ...... . . . . .. . . ..... . ...... ... ... . ....... ..~ ............ . ~,_ . ~ ~ ~ ..... . 1. 2. 3. 4. k Fig. 15 Linear dispersion relations in various schemes. 417
A.2 Speed of Permanent Wave Form The speed of solitary waves for a given amplitude can be found in Ertekin (1984) for GN equation. The speed for Wu's equation is given only in limiting case of small amplitude in Wu & Wu (1982~. But one can derive the relation using the integral invariants of Wu's equation in steady state. For L1 and L2 schemes the speed of permanent wave form can be derived using the dynamic free-surface condition and their two integral invariants, i.e., mass and momentum flux. They are given below: ON: C2 = 1 +A (A.6) Wu : c2 = Am+ 2) '(log(1 + A)l + A) (A.7) L1 : C A + (1 + A)H _ A2 = 0 (A.8) L2 : C A + (1 + A)H _ A2 = 0 (A.9) where CA U = 1 + A H = ~+A+~ (A.10) and where C and A are the normalized speed and amplitude. The above results are plotted in Fig. 16. It is of interest to note that in the speed-amplitude relation for L1 and L2 schemes two different speeds are present at an amplitude (note the dotted lines), or vice versa. Presumably this is due to the presence of waves other than the solitary waves. Fig. 16 also shows that the L1 (or L2) scheme gives the bounded amplitude and speed of permanent wave form with A = 0.600 (or 0.714) for the maximum amplitude, and c2 = 1.473 (or 1.591) for the maximum speed. 0.8 ~ 0.6 C~2 * 0.4 0.2 0.0 ............................ ,' ................. , . ~ . . ... . it, ~ at, ·, , . ....:.. ~ ! ~ .... , // A. ....... :::::::::::1 0.0 0.2 0.4 0.6 0.8 1.0 A Fig. 16 Speed of the permanent wave form in vari- ous schemes. 418
DISCUSSION by R.C. Ertekin I think this paper presents an excellent example of how micro computers can be effectively utilized to solve complicated problems in ~ ~ ~ ~ __ solution to the problem is impressive. I would like to comment on some points in the paper. time domain. Their method of The Green-Naghdi equations that you refer to are the simplest equations in a series of equations (or theories) that can be obtained by the Theory of Directed Fluid Sheets. Shields (Ph. D. Thesis, U. C. Berkely, 1985) have derived the theory II and III equations. Their derivations (with Prof. W. C. Webster), in some sense, resemble the Pohlhausen method in laminar boundary layer theory. I think the higher theory Green-Naghdi equations show double valued phase speed (your Fig.16) like Ll and L2. This double-value of c is not very surprising since the highest wave is not the fatest one as shown by Cokelet and Longuet-Higgins in a J. Fluid Mech. article in 70's. If you used a cubic interpolation function, I believe the maximum amplitude would become closer to A~0.78as shown by Fenton in also a JFM article who used a very high order perturbation expansion. As you commented the Green-Naghdi equations can be obtained by variational methods. Solomon and Miles(1985) derived the G-N equations by using Hamilton's principle (J. Fluid Mechanics). I will be happy to furnish these references in a more complete form if you wish. I am not sure that "dry mode" is a good idea. I am not even sure that it happens in nature. We have experienced such dif ficulties also. By changing x,Ay and At we solved the problem. If you reduce your space and time increments your problem of stability may disappear. In free oscillation tests (Fig.4), the oscillations in the trailing part of the waves are perhaps partially numerical. I remember the thesis of H. Schember (Cal. Inst. of Technology, Ph. D. thesis, 1982) who studied a similar problem on these oscillations. Author's Reply Thank you for your nice comments. We have admired your pioneering research on the Green- Naghdi method applied to the nonlinear free- surface flow problems. We were the comment on the dry mode, we partly agree with Professor Ertekin on the possibility of dry bottom due to the numerical instability. However, we have observed the dry bottom in the experiments when we moved a vertical flat plate in the direction normal to the plate and along the centerline of the tank in a shallow water in a wave tank (22cm X 20cm X 110cm). The dry bottom was possible in nature when the flat plate is moving so fast that the water behind cannot fill the empty space fast enough. Then there should be locally dry bottom in nature. On the last comment about Fig.4 for the free oscillation test, we believe that the trailing part of the smaller waves is due to the dispersion of the initial free-surface elevation (hump) into waves of different wave numbers. 419
