Below are the first 10 and last 10 pages of uncorrected machine-read text (when available) of this chapter, followed by the top 30 algorithmically extracted key phrases from the chapter as a whole.

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

Do not use for reproduction, copying, pasting, or reading; exclusively for search engines.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
SESSION 4
WAVY/FREE-SURFACE FLOW: FIELD-EQUATION METHODS

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
This page in the original is blank.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
A Fast Multigrid Method for Solving the Nonlinear Ship Wave Problem with a Free Surface
J.Farmer, L.Martinelli, and A.Jameson
(Princeton University, USA)
Abstract
This paper presents a finite volume method for the solution of the three dimensional, nonlinear ship wave problem. The method can be used to obtain both Euler and Navier-Stokes solutions of the flow field and the a priori unknown free surface location by coupling the free surface kinematic and dynamic equations with the equations of motion for the bulk flow. The evolution of the free surface boundary condition is linked to the evolution of the bulk flow via a novel iteration strategy that allows temporary leakage of mass through the surface before the solution is converged. The method of artificial compressibility is used to enforce the incompressibility constraint for the bulk flow. A multigrid algorithm is used to accelerate convergence to a steady state. The two-layer eddy viscosity formulation of Baldwin and Lomax is used to model turbulence. The scheme is validated by comparing the numerical results with experimental results for the Wigley parabolic hull and the Series 60, Cb=0.6 hull. Waterline profiles from bow to stern are in excellent agreement with experiment. The computed wave drag compares favorably with experiment. Overall, the present method proves to be accurate and efficient.
1
Introduction
It is well established that a complex interaction exists between the viscous boundary layer and wake of a ship hull and the resulting wave pattern [1, 2]. The existence of two similarity parameters, the Froude number (Fr) and the Reynolds number (Re), which do not scale identically between model and full scale hulls, make it difficult to predict the viscous effect on the wave and total drag through model testing. The ship designer may thus resort to numerical simulation, and a great deal of effort has been devoted toward developing numerical tools capable of simulating the flow field about a translating ship. Some of these tools have met with success in capturing the salient features of the flow field, including the difficult-to-model stern region of ship hulls. However, many of the computational methods developed to date, especially those that include viscous effects and a moving free surface, tend to be very complicated and expensive. Thus, the focus of this work is the development of a fast and robust means to compute either viscous or inviscid flow fields about surface piercing ship hulls, and to make comparisons with experimental data.
The method of Hino [3] is a widely used approach for solving incompressible flow problems. This method takes the divergence of the momentum equation and solves implicit equations at each time step for the pressure and velocity fields such that continuity is satisfied. The method is expensive both because of the need to solve implicit equations by an iterative method and because of the cost of calculating the divergence of the momentum equations in a curvilinear coordinate system. Hino uses a finite difference scheme expressed in body-fitted curvilinear coordinates to discretize the solution domain on and below the free surface. The computational grid is not allowed to move with the free surface so an approximation must be employed to model the free surface boundary conditions. A Baldwin-Lomax turbulence model is used in conjunction with the wall function to model the viscous boundary layer. The scheme is first-order accurate in time and requires 104-plus global iterations to reach steady state for simple hull shapes.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
The method of Miyata et al. (1987) [4] uses a similar velocity and pressure coupling procedure but now the grid is allowed to move with the free surface, providing a more exact treatment of the free surface boundary conditions. A sub-grid-scale turbulence model is employed and computations performed for Reynolds numbers up to 105. As with Hino's method, the time-accurate formulation necessitates several thousand time steps to reach steady state solutions. In a later paper, Miyata et al. (1992) [5] present a finite volume approach that substantially improves the computed results over those obtained using the finite difference approach. Simulations using Reynolds numbers up to 106 were made and more complicated hull shapes examined. The method still requires many thousands of time steps to achieve steady state solutions.
The interactive approach of Tahara et al. [6] uses a field method based on the finite-analytic method used by Chen et al. [7] for the viscous region, and a surface singularity method based on the “SPLASH” panel method of Rosen [8] for the inviscid outer domain. The method iterates between the inviscid and viscous regions by adjusting the small-domain panel distribution to allow for the boundary layer displacement thickness determined from the large-domain solution. The free surface boundary conditions are linearized and applied at the mean water elevation surface. Results of this approach appear to be quite promising for the Wigley hull, and a substantial savings in required computational cost is realized over the large-domain approaches of Hino and Miyata.
In this work, a field method is adopted for the entire flow domain like Hino and Miyata. However, the incompressibility constraint is enforced through the method of artificial compressibility, rather than the velocity-pressure coupling method. The method of artificial compressibility was originally proposed by Chorin [9] in 1967 to solve viscous flows. Since then, Rizzi and Eriksson [10] have applied it to rotational inviscid flow, Dreyer [11] has applied it to low speed two dimensional airfoils and Kodama [12] has applied it to ship hull forms with a symmetric free surface. In addition, Turkel [13] has investigated more sophisticated preconditioners than those originally proposed by Chorin. The basic idea behind artificial compressibility is to introduce a pseudotemporal equation for the pressure through the continuity equation. Use of this pressure equation, rather than the velocity-pressure coupling procedure described in references [3]—[7], renders the new set of equations well conditioned for numerical computation along the same lines as those used to calculate compressible flow about complete aircraft [14,15]. When combined with multigrid acceleration procedures [16,17,18] it proves to be particularly effective. Converged solutions of incompressible flows over three dimensional isolated wings are obtained in 25–50 cycles.
The general objective of this work is to build on these ideas to develop a more efficient method to predict free surface wave phenomena, for both inviscid and viscous flows. The viscous solution method introduced in this work is an extension of the inviscid method presented in reference [19, 20]. The nonlinear free surface boundary condition is satisfied by an iterative procedure in which the grid is moved with the free surface. Comparisons of numerical predictions with experimental data, for the Wigley hull and Series 60, Cb=0.6 ship hull, show encouraging results for both waterline profiles and wave drag. Furthermore, it appears that this approach yields a substantial savings in the computational resources required for the simulations.
2
Mathematical Model
Figure 1 shows the reference frame and ship location used in this work. A right-handed coordinate system Oxyz, with the origin fixed at midship on the mean free surface is established. The z direction is positive upwards, y is positive towards the starboard side and x is positive in the aft direction. The free stream velocity vector is parallel to the x axis and points in the same direction. The ship hull pierces the uniform flow and is held fixed in place, ie. the ship is not allowed to sink (translate in z direction) or trim (rotate in x–z plane).
2.1
Bulk Flow
For a viscous incompressible fluid moving under the influence of gravity, the continuity equation and the Reynolds averaged Navier-Stokes equations may be put in the form [3],
ux+νy+wz=0 (1)

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 1: Reference Frame and Ship Location
νt+uνx+ννy+wνz= −ψy+(Re−1+νt) (▽2ν) (2)
wt+uwx+vwy+wwz= −ψz+(Re−1+νt) (▽2w).
Here, u=u(x,y,z,t), ν=ν(x,y,z,t) and w=w(x,y,z,t) are the mean total velocity components in the x,y and z directions. All lengths and velocities are nondimensionalized by the ship length L and the free stream velocity U, respectively. The pressure ψ is the static pressure p minus the hydrostatic component −zFr−2 and may be expressed as ψ=p+zFr−2, where is the Froude number. The pressure variable ψ is nondimensionalized by ρU2. The Reynolds number Re is defined by where ν is the kinematic viscosity of water and is constant. νt is the dimensionless turbulent eddy viscosity, computed locally using the Baldwin-Lomax turbulence model. This set of equations shall be solved subject to the following boundary conditions.
2.2
Boundary Conditions
2.2.1
Free Surface
When the effects of surface tension and viscosity are neglected, the boundary condition on the free surface consists of two equations. The first, the dynamic condition, states that the pressure acting on the free surface is constant. The second, the kinematic condition, states that the free surface is a material surface: once a fluid particle is on the free surface, it forever remains on the surface. The dynamic and kinematic boundary conditions may be expressed as
p=constant
(3)
where z=β(x,y,t) is the free surface location. Equation 3 only permits solutions where β is single valued. Consequently, it does not allow for the breaking of bow waves which can often be observed with cruiser type hulls. Breaking waves are difficult to treat numerically and are not considered in this work.
2.2.2
Hull and Farfield
The remaining boundaries consist of the ship hull, the boundaries which comprise the symmetry portions of the meridian plane and the far field of the computational domain. On the ship hull, the condition is that of no-slip and is stated simply by
u=ν=w=0.
On the symmetry plane (that portion of the (x,z) plane excluding the ship hull) derivatives in the y direction as well as the ν component of velocity are set to zero. The upstream plane has u=U and ψ=0 (p=−zFr−2) with the ν and w velocity components set to zero. Similar conditions hold on the bottom plane which is assumed to represent infinitely deep water where no disturbances are felt. One-sided differences are used to update the flow variables on the starboard plane. A radiation condition should be imposed on the outflow domain to allow the wave disturbance to pass out of the computational domain. Although fairly sophisticated formulations may be devised to represent the radiation condition, simple extrapolations proved to be sufficient in this work.
2.3
Turbulence Model
To model turbulence in the flow field the laminar viscosity is replaced by
μ=μl+μt
where the turbulent viscosity μt is computed using the algebraic model of Baldwin and Lomax [22]. The Baldwin-Lomax model is an algebraic scheme that makes use of a two-layer, isotropic eddy viscosity formulation. In this model the turbulent viscosity is evaluated using

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
where y is the distance measured normal to the body surface and ycrossover is the minimum value of y where both the inner and outer viscosities match. The inner viscosity follows the Prandtl-Van Driest formula,
where
l=ky[1−exp(−y+/A+)]
is the turbulent length scale for the inner region, k and A+ are model constants, |ω| is the vorticity magnitude and is the dimensionless distance to the wall in wall units.
In the outer region of the boundary layer, the turbulent viscosity is given by
(μt)outer=KCcpFwakeFKleb
where K and Ccp are model constants, the function Fwake is
and the function FKleb is
.
The quantities Fmax and ymax are determined by the value and corresponding location, respectively, of the maximum of the function
F=y|ω|[1−exp(−y+/A+)].
The quantity Udif is the difference between maximum and minimum velocity magnitudes in the profile and is expressed as
CKleb and Cwk are additional model constants. Numerical values for the model constants used in the computations are listed here:
A+=26, k=0.4, K=0.0168,
and
Ccp=1.6, Cwk=1.0, CKleb=0.3.
3
Numerical Solution
The formulation of the numerical solution procedure is based on a finite volume method (FVM) for the bulk flow variables (u,ν,w and ψ), coupled to a finite difference method for the free surface evolution variables (β and ψ). Alternative cell-centered and cell-vertex formulations may be used in finite volume schemes [16]. A cell-vertex formulation was preferred in this work because values of the flow variables are needed on the boundary to implement the free surface boundary condition. The bulk flow is solved subject to Dirichlet conditions for the free surface pressure, followed by a free surface update via the bulk flow solution (ie. constant values for the velocities in equation 3). Each formulation is explicit and uses local time stepping. Both multigrid and residual averaging techniques are used in the bulk flow to accelerate convergence.
3.1
Bulk Flow Solution
Following Chorin [9] and more recently Yang et al. [23], the gover ning set of incompressible flow equations may be written in vector form as
wt*+(f−fv)x+(g−gv)y+(h−hv)z=0 (4)
where the vector of dependent variables w and inviscid flux vectors f, g and h are given by
The viscous flux vectors fv, gv and hv are given by
.
where the viscous stress components are defined as
.
Г is called the “artificial compressibility” parameter due to the analogy that may be drawn between the above equations and the equations of

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
motion for a compressible fluid whose equation of state is given by [9]
ψ=Г2ρ.
Thus, ρ is an artificial density and Г may be referred to as an artificial sound speed. When the temporal derivatives tend to zero, the set of equations satisfy precisely the incompressible equations 2, with the consequence that the correct pressure may be established using the artificial compressibility formulation. The artificial compressibility parameter may be viewed as a device to create a well posed system of hyperbolic equations that are to be integrated to steady state along lines similar to the well established compressible flow FVM formulation [18]. In addition, the artificial compressibility parameter may be viewed as a relaxation parameter for the pressure iteration. Note that temporal derivatives are now denoted by t* to indicate pseudo time; the artificial compressibility, as formulated in the present work, destroys time accuracy.
To demonstrate the effect of Г on the above set of equations and to establish the hyperbolicity of the set, the convective part of equation 4 may be written in quasi-linear form to determine the eigenvalues [10]. The eigenvalues are found to be
λ1=U, λ2=U, λ3=U+a, λ4=U−a,
where
U=uωx +νωy+wωz
and
The wave number components ωx, ωy and ωz are defined on −∞≤ωx, ωy, ωz≤+∞. Since the eigenvalues are clearly real for any value of ωx, ωy and ωz, the system of equations 4 is hyperbolic.
The choice of Г is crucial in determining convergence and stability properties of the numerical scheme. Typically, the convergence rate of the scheme is dictated by the slowest system waves and the stability of the scheme by the fastest. In the limit of large Г the difference in wave speeds can be large. Although this situation would presumably lead to a more accurate solution through the “penalty effect” in the pressure equation, very small time steps would be required to ensure stability. Conversely, for small Г, the difference in the maximum and minimum wave speeds may be significantly reduced, but at the expense of accuracy. Thus a compromise between the two extremes is required. Following the work of Dreyer [11], the choice for Г is taken to be
Г2=γ(u2+v2+w2),
where γ is a constant of order unity. In regions of high velocity and low pressure where suction occurs, Г is large to improve accuracy, and in regions of lower velocity, Г is correspondingly reduced.
The choice of Г also influences the outflow boundary condition, or radiation condition. If it can be demonstrated that all system eigenvalues are both real and positive, then downstream or outflow boundary points may be extrapolated from the interior upstream flow. Even though an examination of the eigenvalues reveals that this can never be the case, the condition can be approached by a judicious choice of Г. If Г is large, extrapolation fails because the flow has both downstream and upstream dependence. As Г is reduced, the upstream dependence becomes more pronounced and the downstream is reduced. Eventually, the upstream dependence is sufficiently dominant to allow extrapolation. Hence, all outflow variables are updated using zero gradient extrapolation.
Following the general procedures for FVM, the governing equations may be integrated over an arbitrary volume Λ. Application of the divergence theorem on the convective and viscous flux term integrals yields
(5)
where Sx, Sy and Sz are the projected areas in the x, y and z directions, respectively. The computational domain is divided into hexahedral cells. Application of FVM to each of the computational cells results in the following system of ordinary differential equations,
The volume Λijk is given by the summation of the eight cells surrounding node i,j,k. The convective flux Cijk(w) is defined as
(6)

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
and the viscous flux Vijk(w) is defined as
(7)
where the summation is over the n faces surrounding Λijk.
The projected areas may be computed by taking the cross product of the two vectors joining opposite corners of each cell face in the physical coordinate system. They correspond to the grid metrics Jξx, Jξy, Jξz, etc. appearing in a transformation to a curvilinear coordinate system ξ= ξ(x,y,z), η=η(x,y,z) and ζ=ζ(x,y,z) where J is the Jacobian of the transformation. The flow variables required in the flux evaluation may be averaged on each cell face through the four nodal values associated with each face. Evaluation of the flux terms in equations 6 and 7 may be performed directly, without direct differentiation and without the need to handle grid singularities in a special fashion.
3.1.1
Artificial Dissipation
This scheme reduces to a second order accurate, nondissipative central difference approximation to the bulk flow equations on sufficiently smooth grids. A central difference scheme permits odd-even decoupling at adjacent nodes which may lead to oscillatory solutions. To prevent this “unphysical” phenomena from occurring, a dissipation term is added to the system of equations such that the system now becomes
. (8)
For the present problem a third order background dissipation term is added. The dissipative term is constructed in such a manner that the conservation form of the system of equations is preserved. The dissipation has the form
Dijk(w)=Dξ+Dη+Dζ (9)
where
and
(10)
Similar expressions may be written for the η and ζ directions with , and representing second difference central operators.
In equation 10, the dissipation coefficient α is a scaling factor proportional to the local wave speed, and renders equation 9 third order in truncation terms so as not to detract from the second order accuracy of the flux discretization. The actual form for the coefficient is based on the spectral radius of the system and is given in the ξ direction as
where ũ is the contravariant velocity component
ũ=uξx+νξy+wξz.
Similar dissipation coefficients are used for the η and ζ components in equation 9. The ∊ term is used to manually adjust the amount of dissipation.
3.1.2
Viscous Discretization
The discretization for the viscous fluxes follows the guidelines originally proposed in [24,25] for the simulation of two dimensional viscous flows. The components of the stress tensor are computed at the cell centers with the aid of Gauss' formula. The viscous fluxes are then computed by making use of an auxiliary cell bounded by the faces lying on the planes containing the centers of the cells surrounding a given vertex and the mid-lines of the cell faces. For example, the ux term in may be computed from
(10)
where k=1,6 are the six faces surrounding a particular cell, uk is an average of the velocities from the nodes that define the kth face and Szk are the projected areas in the x direction corresponding to each face. Once the components of the complete stress tensor are computed at the centroids of the cells then the same method of evaluation may be used to compute the viscous fluxes at the vertex through use of equation 7. For this purpose the control volume is now constructed by assembling fractions of each of the eight cells surrounding a particular vertex. The equivalent two dimensional control volume is sketched in the figure below. This discretization procedure is designed to minimize the error induced by a kink in the grid. It has proved to be accurate and efficient in applications to the solution of three dimensional compressible viscous flows [26,27].

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
3.1.3
Time Integration
Equation 8 is integrated in time by an explicit multistage scheme. For each bulk flow time step, the grid, and thus Λijk, is independent of time. Hence equation 8 can be written as
(11)
where the residual is defined as
Rijk(w)=Cijk(w)−Vijk(w)−Dijk(w),
and the cell volume Λijk absorbed into the residual for clarity. If one analyzes a linear model problem corresponding to (11) by substituting a Fourier mode , the resulting Fourier symbol has an imaginary part proportional to the wave speed, and a negative real part proportional to the diffusion. Thus the time stepping scheme should have a stability region which contains a substantial interval of the negative real axis, as well as an interval along the imaginary axis. To achieve this it pays to treat the convective and dissipative terms in a distinct fashion. Thus the residual is split as
Rijk(w)=Cijk(w)+Dijk(w)
where Cijk(w) is the convective part and Dijk(w)=−(Vijk+Dijk) the dissipative part. Denote the time level nΔt by a superscript n, and drop the subscript for clarity. Then the multistage time stepping scheme is formulated as
w(n+1,0)=wn
…
w(n+1,k)=wn−αkΔt(C(k−1)+D(k−1))
…
wn+1=w(n+1,m)
where the superscript k denotes the k-th stage, αm=1, and
C(0)=C(wn), D(0)=D(wn)
The coefficients αk are chosen to maximize the stability interval along the imaginary axis, and the coefficients βk are chosen to increase the stability interval along the negative real axis.
A five-stage scheme with three evaluations of dissipation has been found to be particularly effective. Its coefficients are
α1=1/4 β1=1
α2=1/6 β2=0
α3=3/8 β3=0.56
α4=1/2 β4=0
α5=1 β5=0.44
The actual time step Δt is limited by the Courant number (CFL), which states that the fastest waves in the system may not be allowed to propagate farther than the smallest grid spacing over the course of a time step. In this work, local time stepping is used such that regions of large grid spacing are permitted to have relatively larger time steps than regions of small grid spacing. Of course the system wave speeds vary locally and must be taken into account as well. The final local time step is thus computed as,
where λijk is the sum of the spectral radii of both the convective and viscous flux Jacobian matrices in the x, y and z directions. In regions of small grid spacing and/or regions of high characteristic wave speeds, the time step will be smaller than elsewhere.
3.1.4
Residual Averaging
The allowable Courant number may be increased by smoothing the residuals at each stage using the following product form in three dimensions [18]
where and are smoothing coefficients and the are central difference operators in computational coordinates. Each residual Rijk is thus replaced by an average of itself and the neighboring residuals.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
3.1.5
Multigrid Scheme
Very rapid convergence to a steady state is achieved with the aid of a multigrid procedure. The idea behind the multigrid strategy is to accelerate evolution of the system of equations on the fine grid by introducing auxiliary calculations on a series of coarser grids. The coarser grid calculations introduce larger scales and larger time steps with the result that low-frequency error components may be efficiently and rapidly damped out. Auxiliary grids are introduced by doubling the grid spacing, and values of the flow variables are transferred to a coarser grid by the rule
where the subscripts denote values of the grid spacing parameter (ie. h is the finest grid, 2h, 4h, …are successively coarser grids) and T2h,h is a transfer operator from a fine grid to a coarse grid. The transfer operator picks flow variable data at alternate points to define coarser grid data as well as the coarser grid itself. A forcing term is then defined as
where R is the residual of the difference scheme. To update the solution on the coarse grid, the multistage scheme is reformulated as
…
…
where R(q) is the residual of the qth stage. In the first stage, the addition of P2h cancels R2h(w(0)) and replaces it by ΣRh(wh), with the result that the evolution on the coarse grid is driven by the residual on the fine grid. The result now provides the initial data for the next grid and so on. Once the last grid has has been reached, the accumulated correction must be passed back through successively finer grids. Assuming a three grid scheme, let represent the final value of w4h. Then the correction for the next finer grid will be
where Ia, b is an interpolation operator from the coarse grid to the next finer grid. The final result on the fine grid is obtained in the same manner:
The process may be performed on any number of successively coarser grids. The only restriction in the present work being use of a structured grid whereby elements of the coarsest grid do not overlap the ship hull. A 4-level “W-cycle” is used in the present work for each time step on the fine grid [18].
3.1.6
Grid Refinement
The multigrid acceleration procedure is embedded in a grid refinement procedure to further reduce the computer time required to achieve steady state solutions on finely resolved grids. In the grid refinement procedure the flow equations are solved on coarse grids in the early stages of the simulation. The coarse grids permit large time steps, and the flow field and the wave pattern evolve quite rapidly. When the wave pattern approaches a steady state, the grid is refined by doubling the number of grid points in all directions and the flow variables and free surface location are interpolated onto the new grid. Computations then continue using the finer grid with smaller time steps. The multigrid procedure is applied at all stages of the grid refinement to accelerate the calculations on each grid in the sequence, producing a composite “full multigrid” scheme which is extremely efficient.
3.2
Free Surface Solution
Both a kinematic and dynamic boundary condition must be imposed at the free surface. For the fully nonlinear condition, the free surface must move with the flow (ie. up or down corresponding to the wave height and location) and the boundary conditions applied on the distorted free surface. Equation 3 can be cast in a form more amenable to numerical computations by introducing a curvilinear coordinate system that transforms the curved free surface β(x,y) into computational coordinates β(ξ,η). This results in the following transformed kinematic condition
(12)
where ũ and are contravariant velocity components given by

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
The free surface kinematic equation may now be expressed as
where Qij(β) consists of the collection of velocity and spacial gradient terms which result from the discretization of equation 12. Note that this is not the result of a volume integration and thus the volume (or actually area) term does not appear in the residual as in the FVM formulation. Throughout the interior of the (x,y) plane, all derivatives are computed using the second order centered difference stencil in computational coordinates ξ and η. On the boundaries a second order centered stencil is used along the boundary tangent and a first order one sided difference stencil is used in the boundary normal direction.
As was necessary in the FVM formulation for the bulk flow, background dissipation must be added to prevent decoupling of the solution. The method used to compute the dissipative terms borrows from a two dimensional FVM formulation and appears as follows:
Dij=Dξ+Dη
where
Dξij=dξi+1,j−dξi,j
and
The expression for α may be written as
where J is the sum of the cell Jacobians and is used to manually adjust the amount of dissipation. Hence the system of equations for the free surface is expressed as
where
Rij=Qij−Dij.
The same scheme used to integrate equation 11 is also used here. Once the free surface update is accomplished the pressure is adjusted on the free surface such that
The free surface and the bulk flow solutions are coupled by first computing the bulk flow at each time step, and then using the bulk flow velocities to calculate the movement of the free surface. After the free surface is updated, its new values are used as a boundary condition for the pressure on the bulk flow for the next time step. The entire iterative process, in which both the bulk flow and the free surface are updated at each time step, is repeated until some measure of convergence is attained; usually steady state wave profile and wave resistance coefficient.
Since the free surface is a material surface, the flow must be tangent to it in the final steady state. During the iterations, however, the flow is allowed to leak through the surface as the solution evolves towards the steady state. This leakage, in effect, drives the evolution equation. Suppose that at some stage, the vertical velocity component w is positive (cf. equation 3 or 12). Provided that the other terms are small, this will force βn+1 to be greater than βn. When the time step is complete, ψ is adjusted such that ψn+1>ψn. Since the free surface has moved farther away from the original undisturbed upstream elevation and the pressure correspondingly increased, the velocity component w (or better still q · n where and F=z−β(x,y)) will then be reduced. This results in a smaller Δβ for the next time step. The same is true for negative vertical velocity, in which case there is mass leakage into the system rather than out. Only when steady state has been reached is the mass flux through the surface zero and tangency enforced. In fact, the residual flux leakage could be used in addition to drag components and pressure residuals as a measure of convergence to the steady state.
This method of updating the free surface works well for the Euler equations since tangency along the hull can be easily enforced. However, for the Navier-Stokes equations the no-slip boundary condition is inconsistent with the free surface boundary condition at the hull/waterline intersection. To circumvent this difficulty the computed elevation for the second row of grid points away from the hull is extrapolated to the hull. Since the minimum spacing normal to the hull is small, the error due to this should be correspondingly small, comparable with other discretization errors. The treatment of this intersection for the Navier-Stokes calculations, should be the subject of future research to find the most accurate possible procedure.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 6: Mesh deformation at 8.25, 8.5, 8.75 and 9th cycles.
is two-folded. First, a significant reduction of the amplitude of the hydrodynamic force is observed. Also, a decrease in the difference between positive and negative peaks of this force occurs.
6.2
Heave of submerged cylinder in a tank
Figure 10 shows a typical configuration for the mesh corresponding to the submerged cylinder. For this case, the cylinder radius is unity. The variable NI appearing in this section corresponds to the number of divisions used to discretize the free surface and a quarter of the body surface. Based on the experience gained in the case of the surface-piercing cylinder, in this section we have opted to use quadratic elements only. No smoothing was required for this case.
Figures 11 and 12 respectively show linear and nonlinear free surface time histories for A=0.5. Vertical dimensions have also been magnified by 10 times. Steeper peaks were observed
Figure 7: Velocity vector plots at 8.25, 8.5, 8.75 and pth cycles.
Figure 8: Hydrodynamic forces for surface-piercing cylinder.
to develop in the regions next to the tank walls for nonlinear results. Linear results predicted an almost sinusoidal behavior of the free surface across the width of the tank. As expected, vertical elevations remained lower than that obtained for the case of the surface-piercing cylinder.
Figure 13 shows mesh deformation corre-

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 10: Hydrodynamic forces for surface-peiercing cylinder versus tank dimensions.
Figure 10: Typical mesh for submerged cylinder.
Figure 11: Linear free surface for submerged cylinder (NI=20, Δt=0.1).
Figure 12: Nonlinear free surface time record for submerged cylinder (NI=10, t=0.05).
Figure 13: Mesh deformation at 9.25, 9.5, 9.75 and 10th cycles.
sponding to the 9.25, 9.5, 9.75 and 10th cycles of oscillation. Figure 14 shows corresponding velocity vector plots. Again, velocity fields behaved in a physically acceptable manner with very small differences between specified and calculated values.
For A=0.5, linear and nonlinear prediction of the hydrodynamic forces resulting from vertical motion of the circular cylinder displayed similar trends. Figure 15 shows linear and nonlinear hydrodynamic forces corresponding to an amplitude of motion A=1.5. Even for this extremely large amplitude of motion, both formulations gave similar results. Some differences

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 14: Velocity vector plots at 9.25, 9.5, 9.75 and 10th cycles.
Figure 15: Hydrodynamic force due to sinusoidal excitation, A=1.5.
Figure 16: Free surface frequency spectra due to impulsive motion for circular basin.
can again be seen in the force peaks, but these are significantly less noticeable than for the surface piercing case.
6.3
Impulsive motion in a tank
In this section, we present results corresponding to free surface response due to impulsive excitation of tanks of different shapes. Excitation is prescribed as a step velocity at t=0. We have opted to match the types of excitations given by Yeung and Wu[47]. Application of Fast Fourier Transform to the time record of the free surface elevation measured at a fixed point yields the frequency spectrum. Nondimensional frequencies for linear and nonlinear free surface boundary conditions are presented. Smoothing in the worst of these cases was applied using α=1 every 10th time step.
Figure 16 shows linear and nonlinear frequency spectra results for a circular tank. In this case, frequencies predicted by linear theory tend to be slightly higher than those given by nonlinear theory. Table 1 shows comparison with results given by Chang and Wu[70] (CW) using an integral equation formulation with linear boundary conditions. Their results also exhibit higher values than those predicted by nonlinear theory.
Corresponding results for a parabolic basin are given in Figure 17. For this particular configuration, results tend to show no difference be-

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
tween linear and nonlinear predictions. Table 2 compares our results with those given by Yeung and Wu[47]. In general, good agreement is achieved.
Table 1 Nondimensional frequencies for circular tank
CIRCULAR TANK
Freq.
Lin.
Nonl.
CW
1st
1.7330
1.7180
1.7492
2nd
2.4697
2.4545
2.5632
3rd
3.0066
2.9759
—
Table 2 Nondimensional frequencies for parabolic tank
PARABOLIC TANK
Freq.
Lin.
Nonl.
YW. Lin.
YW. Nonl
1st
1.6720
1.6720
1.6846
1.6846
2nd
2.4390
2.4390
2.4170
2.4170
3rd
3.0066
3.0066
3.0029
3.0029
Figure 18 shows results for a triangular basin. For this case, no definite trend is shown between linear and nonlinear results. As noted in Table 3, this effect is also shown in the results of Yeung and Wu[47]. Results by Lamb[71] are also presented.
Table 3 Nondimensional frequencies for triangular tank
TRIANGULAR TANK
Freq.
Lin.
Nonl.
YW. Lin.
YW. Nonl
Lamb
1st
1.5339
1.5646
1.5381
1.5381
1.5244
2nd
2.3930
2.3930
2.3439
2.3438
2.3447
3rd
3.0066
2.9912
2.9297
3.0029
2.9393
Figure 17: Free surface frequency spectra due to impulsive motion for parabolic basin.
6.4
Slender body calculations
In this section we present results corresponding to the forward motion of a mathematically defined ship hull advancing at zero angle of attack in otherwise calm water. The Wigley parabolic hull, whose main parameters are given in Table 4, is defined as
(19)
Table 4 Wigley hull dimensions
WIGLEY PARABOLIC HULL[m]
Length L
Beam B
Draft D
CB
CWP
CP
2
0.2
0.125
0.444
0.667
0.667
Coordinates z and y remained as defined in pre-

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 18: Free surface frequency spectra due to impulsive motion for triangular basin.
vious sections for surface-piercing cases. The coordinate x is directed along the length of the hull and considered positive towards the stern.
Results presented for trim, sinkage and wave-resistance coefficients have been obtained for the conditions shown in Table 5. The tank width is in metres. The variable Δx refers to the spatial increment used along the x-coordinate in metres. For the length of the Wigley hull used in this study, this increment resulted in 2000 stations. CPU time estimates are presented in Table 6. Smoothing in the worst of these cases was also applied using α=1 every 10th time step.
Table 5 Computational conditions for Wigley hull calculations
COMPUTATIONAL CONDITIONS
Δx
Tank Width
No of elements
No of Nodes
0.001
3.0
1200
3741
Table 6 Computational demand for Wigley hull calculations
COMPUTATIONAL DEMANDS
Machine Type
Mflops
CPU [min/Δx]
SUN MP/670
4.0
5.0
IBM RS/6000
31.0
0.75
Figures 19 and 20 show linear and nonlinear three dimensional views of the free surface elevation for the Wigley hull advancing with a forward speed of . This corresponds with a length Froude number of 0.267. It is clear that nonlinear theory predicts a more realistic shape of the free surface. This is evident in the aft-half of the hull where linear theory is not able to produce the recovery from the second crest of the wave next to the hull. This is more easily seen in Figure 21 where comparison has been made with results presented by Maruo and Song[32]. Nonlinear results show good agreement with experimental values for the wave past the midship section. However, there is a shift of the first crest and trough of the fore-half of the wave. It is conjectured that one possible reason for this shift is due to the at rest (zero-valued) initial conditions used in this work. Figure 21 shows numerical results given by Maruo and Song[32] which clearly display non-homogeneous initial conditions. Free surface elevation contours are also shown in Figure 22 and Figure 23 for linear and nonlinear boundary conditions respectively. Nonlinear results show a ‘jagged' nature due to small wiggles appearing on the free surface. These wiggles were typical of nonlinear free surface profiles and occasionally were the origin of unstable behavior. In order to maintain the size of these wiggles bounded, small time steps (Δx) were used.
Figure 24 shows the wave resistance coefficient as a function of the Froude number Fn for linear and nonlinear formulations. Data given by Aanesland[72] and mean values of experimental data from the 1st and 2nd Workshops on Wave Resistance Computations are presented [73,74]. At low values of Fn nonlinear prediction shows significant improvement over results

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 19: Linear free surface elevation for Wigley hull with forward speed, Fn=0.267.
obtained using linear theory. In fact, linear theory gives very poor agreement with experimental results for this Fn-range. For this set of experimental data, nonlinear results show superior performance throughout the Fn interval. Figure 25 shows the wave resistance coefficient compared with Maruo and Song[ 32]. Again, at low values of the Froude number nonlinear results give a better prediction than linear ones. However, as the Froude number increases, linear results show improvement over their nonlinear counterpart. This is a clear example of the scatter encountered in experimental data for these type of physical problems.
Figures 26 and 27 show sinkage and trim of the Wigley hull as a function of Fn. Improved predictions at the lower end of the Froude number interval used in this work are again realized for nonlinear results. This improvement is more clearly seen in the results for sinkage. These results suggest that nonlinear theory predicts a more accurate magnitude of the forces. This is consistent with results obtained for wave resistance as shown in Figure 24. Trim calculations show little discrepancy between linear and nonlinear results. In turn, this suggests that the distribution of forces is predicted with similar accuracy in both formulations.
CONCLUSIONS
The boundary value problem defined by the Laplace equation for the solution of potential
Figure 20: Non Linear free surface elevation for Wigley Hull with forward speed Fn=0.267
Figure 21: Wave profile at side of Wigley Hull with forward speed, Fn=0.267
Figure 22: Linear free surface elevation contours of Wigley hull with forward speed, Fn=0.267.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
Figure 23: Nonlinear free surface elevation contours of Wigley Hull with forward speed, Fn=0.267
Figure 24: Wave resistace coefficient of Wigley hull versus Fn.
Figure 25: Wave resistance coefficient of Wigley hull versus Fn.
Figure 26: Sinkage of Wigley hull versus Fn.
Figure 27: Trim of Wigley hull versus Fn.
flow in two dimensions has been formulated using a Bubnov-Galerkin weighted residual method to obtain the spatial variation of the velocity potential. The moving boundary problem corresponding to the a-priori unknown dependent and independent spatial variables in the nonlinear boundary conditions has been treated using a semi-Lagrangian semi-implicit scheme. This method proved to be more accurate and stable than predictor-corrector methods. When necessary, dissipation was applied only to the velocity potential.
ACKNOWLEDGEMENT
I (AA) would like to extend my sincere gratitude to Dr. Rodolfo Bermejo for his generous assistance during the development of some of the concepts expressed in this paper. I am also thankful to Dr. Chun-Fa Wu for his help during our phone conversations.

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
REFERENCES
[1] Van Dyke M.D., ‘Perturbation Methods in Fluid Mechanics', Academic Press, 1964.
[2] Cummins W.E., ‘The Wave Resistance of a Floating Slender Body', Unpublished Thesis, American University, 1956.
[3] Kaplan P., ‘Application of Slender Body Theory to Forces Acting on Submerged Bodies and Surface Ships in Regular Waves', Journal of Ship Research, vol. 1, 1957.
[4] Lighthill M.J., ‘Note on the Swimming of Slender Fish', Journal of Fluid Mechanics, vol. 9, 1960.
[5] Vossers G., ‘Wave Resistance of a Slender Ship', Schiffstechnik, vol. 9, 1962.
[6] Maruo H., ‘Calculation of the Wave Resistance of Ships, the draft of which is as small as the beam', Journal of the Society of Naval Architects of Japan, vol. 112 , 1962.
[7] Ursell F., ‘Slender Oscillating Ships At Zero Forward Speed', Journal of Fluid Mechanics, vol. 14, 1962.
[8] Tuck E.O., ‘A Systematic Asymptotic Expansion Procedure for Slender Ships', Journal of Ship Research, vol. 8, 1964.
[9] Newman J.N. and Tuck E.O., ‘Current Progress in the Slender Body Theory for Ship Motions', Procceedings of the 5th Symposium on Naval Hydrodyanmics, 1964.
[10] Newman J.N., ‘The exciting Forces on a Moving Body in Waves', Journal of Ship Research, vol. 9, 1965.
[11] Haskind M.D., ‘The exciting Forces and Wetting of a Ship in Waves', Izvestia Akademii Nauk SSSR, 1957 (in russian). DTMB translation No 307, 1962.
[12] Tuck E.O., ‘Shallow-Water Flows Past Slender Bodies', Journal of Fluid Mechanics, vol. 26, 1966.
[13] Newman J.N., ‘Applications of Slender Body Theory in Ship Hydrodynamics', Annual Review of Fluid Mechanics, vol. 2, 1970.
[14] Chapman R.B., ‘Free Surface Effects for Yawed Surface Piercing Plates', Journal of Ship Research, vol. 20 , 1976.
[15] Jensen J.J. and Pedersen P.T., ‘Wave-Induced Bending Moments in Ships', Transactions of the Royal Institute of Naval Architects, vol. 121 , 1979.
[16] Troesch A.W., ‘Sway, Roll and Yaw Motion Coefficients Based on a Forward Speed Slender Body Theory: Part 1', Journal of Ship Research, vol. 25, 1981.
[17] Troesch A.W., ‘Sway, Roll and Yaw Motion Coefficients Based on a Forward Speed Slender Body Theory: Part 2', Journal of Ship Research, vol. 25, 1981.
[18] Maruo H., ‘New Approach to the Theory of Slender Ships with Forward Speed', Bulletin of the Faculty of Engineering, Yokohama University, vol. 31, 1982.
[19] Noblesse F., ‘A Slender Ship Theory of Wave Resistance', Journal of Ship Research, vol. 27, 1983.
[20] Faltinsen O.M., ‘Bow Flow and Added Resistance of Slender Ships at High Froude Numbers and Low Wave Lengths', Journal of Ship Research, vol. 27, 1983.
[21] Chen C.Y. and Noblesse F., ‘Preliminary Numerical Study of a New Slender Ship Theory of Wave Resistance', Journal of Ship Research, vol. 27, 1983.
[22] Chen C.Y. and Noblesse F., ‘Comparison Between Theoretical Predictions of Wave Resistance and Experimental Data for the Wigley Hull', Journal of Ship Research, vol. 27, 1983.
[23] Sclavounos P.D., ‘The Diffraction of Free Surface Waves by a Slender Ship', Journal of Ship Research, vol. 28, 1984.
[24] Sclavounos P.D., ‘The Unified Slender Body Theory: Ship Motions in Waves', Proceedings of the 15th Symposium on Naval Hydrodynamics, 1984.
[25] Breit S.R. and Sclavounos P.D., ‘Wave Interactions Between Adjacent Slender Bodies', Journal of Fluid Mechanics, vol. 165, 1986.
[26] Song W., Ikehata M. and Suzuki K., ‘Computation of Wave Resistance and Ship Wave Pattern by the Slender Body Approximation', Journal of the Kamsai Society of Naval Architects of Japan, vol. 209, 1988, in Japanese.
[27] Choi H.S. and Mei C.C., ‘Wave Resistance

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
and Squat of a Slender Ship Moving Near the Critical Speed in Restricted Waters', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[28] Kashiwagi M., ‘Theoretical prediction of Tank Wall Effects on Hydrodynamic Forces Acting on an Oscillating and Translating Slender Ship', Extended Abstracts of the 4th International Workshop on Water Waves and Floating Bodies, 1989.
[29] Kashiwagi M., ‘Side Wall Effects on Hydrodynamic Forces Acting on a Ship with Forward and Oscillating Motions', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[30] Calisal S.M. and Chan J.L.K., ‘A Numerical Model of Ship Bow Waves', Journal of Ship Research, vol. 33, 1989.
[31] Ando S., ‘Forward Speed Effect in Head Wave Diffraction over the Forebody of A Slender Hull', Extended Abstracts of the 4th International Workshop on Water Waves and Floating Bodies, 1989.
[32] Maruo H. and Song W., ‘A Numerical Appraisal of the New Slender Ship Formulation in Steady Motion', Proceedings of the 18th Symposium on Naval Hydrodynamics, 1990.
[33] Yeung R.W., ‘Numerical Methods in Free Surface Flows', Annual Review of Fluid Mechanics, vol. 14, 1982.
[34] Longuet-Higgins M.S. and Cokelet E.D., ‘The Deformation of Steep Surface Waves on Water I: A Numerical Method of Computation'. Proceedings of the Royal Society of London, Series A 350, 1976.
[35] Vinje T. and Brevig P., ‘Nonlinear Ship Motion', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics , 1981.
[36] Lin W.M, Newman J.N. and Yue D.K., ‘Nonlinear Forced Motions of Floating Bodies', Proceedings of the 15th International Symposium on Naval Hydrodynamics , 1985.
[37] Baker G.R., Meiron D.J. and Orzag S.A., ‘Generalized Vortex Methods for Free-Surface Flow Problems', Journal of Fluid Mechanics, vol. 123, 1982.
[38] Dold J.W. and Peregrine D.H., ‘Steep Unsteady Water Waves: An Efficient Computational Scheme', Proceedings of the 19th International Conference on Coastal Engineering, ASCE, 1984.
[39] Suzuki K., ‘Calculation of Nonlinear Water Waves Around a Two-Dimensional Body in Uniform Flow by Means of a Boundary Element Method', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[40] Sen, D., Pawlowsky J.S., Lever J. and Hinchey M.J., ‘Two-Dimensional Numerical Modelling of Large Motions of Floating Bodies in Waves', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[41] Cointe R., ‘Nonlinear Simulation of Transient Free Surface Flows', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[42] Cointe R., Geyer P., King B., Molin B. and Tramoni M., ‘Nonlinear and Linear Motions of a Rectangular Barge in a Perfect Fluid', Proceedings of the 18th International Symposium on Naval Hydrodynamics , 1990.
[43] Ghia U., Shin C.T. and Ghia K.N., ‘Analysis of a Breaking Free Surface Wave Using Boundary Fitted Coordinates for Regions Including Reentrant Boundaries', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981.
[44] Haussling A.J., ‘Solution of Nonlinear Water Wave Problems Using Boundary Fitted Coordinate Systems', Numerical Grid Generation: J.Thompson, editor, 1982.
[45] Telste J.G., ‘Calculation of Fluid Motion from Large-Amplitude Forced Heave Motion of a Two-Dimensional Cylinder in a Free Surface', Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics, 1985.
[46] Yeung R.W. and Wu C., ‘Nonlinear Wave-Body Motion in a Closed Domain', Computers & Fluids, vol. 17, 1989.
[47] Yeung R.W. and Wu C. , ‘On Nonlinear Wave Motion in a Closed Domain', Jahrbuch

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
der Schiffbautechnischen Gesellschaft, vol. 83, 1989.
[48] Luke J.C., ‘A Variational Principle for a Fluid with a Free Surface', Journal of Fluid Mechanics, vol. 27, 1967.
[49] Whitman G.B., ‘Nonlinear Dispersion of Water Waves', Journal of Fluid Mechanics, vol. 27, 1967.
[50] Whitman G.B., ‘Two-timing, Variational Principles and Waves', Journal of Fluid Mechanics, vol. 44, 1970.
[51] K.J.Bai and R.W.Yeung, ‘Numerical Solutions to Free Surface Flow Problems', Proceedings of the 10th Symposium on Naval Hydrodynamics, 1974.
[52] Bai K.J., ‘Diffraction of Oblique Waves by an Infinite Cylinder', Journal of Fluid Mechanics, vol. 68, 1975.
[53] Yim B., ‘A Variational Principle Associated with a Localized Finite Element Technique for Steady Ship-Wave and Cavity Problems', Proceedings of the 1st International Conference on Numerical Ship Hydrodynamics, 1975.
[54] Bai K.J., ‘A Localized Finite Element Method for Steady, Two Dimensional Free Surface Flow Problems', Proceedings of the 1st International Conference on Numerical Ship Hydrodynamics, 1975.
[55] Bai K.J., ‘A Localized Finite Element Method for Two-Dimensional Steady Potential Flows with a Free Surface', Journal of Ship Research, vol. 22, 1978.
[56] Bai K.J., ‘Blockage Correction with a Free Surface', Journal of Fluid Mechanics, vol. 94, 1979.
[57] Wellford C.L. and Ganaba T., ‘Finite Element Procedures for Fluid Mechanics Problems Involving Large Free Surface Motion', Proceedings of the 3rd International Conference in Finite Elements in Flow Problems, 1980.
[58] Oomen A., ‘Free Surface Potential Flow Computation Using a Finite Element Method ', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981.
[59] Jami A., ‘Numerical Solving of Transient Linear Hydrodynamics Problems by Coupling Finite Elements and Integral Representation', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981.
[60] Euvard D., Jami A., Lenoir M. and Martin D., ‘Recent Progress Towards an Optimal Coupling between Finite Elements and Singularity Distribution Procedures', Proceedings of the 3rd International Conference on Numerical Ship Hydrodynamics, 1981.
[61] Allievi A. and Calisal S.M., ‘Finite Element Application to Grid Generation for Submerged and Floating Bodies', International Symposium on Hydro & Aero Dynamics In Marine Engineering, 1991.
[62] Allievi A. and Calisal S.M., ‘A Bubnov-Galerkin Formulation For Orthogonal Grid Generation', Journal of Computational Physics, vol. 98, 1992.
[63] Allievi A., ‘On Nonlinear Free Surface Potential Flow by a Bubnov-Galerkin Formulation in Space and a Semi-Lagrangian Semi-implicit scheme in Time', Ph.D. thesis, University of British Columbia, in preparation.
[64] Robert A.J., ‘A Stable Numerical Integration Scheme for the Primitive Meteorological Equations', Atmosphere-Ocean, vol. 19, 1981.
[65] Robert A.J., ‘A Semi-Lagrangian and Semi-Implicit Numerical Integration Scheme for the Primitive Meteorological Equations', Journal of the Meteorological Society of Japan, vol. 60, 1982.
[66] Staniforth A. and Temperton C., ‘Semi-Implicit Semi-Lagrangian Integration Schemes for a Barotropic Finite-Element Regional Model', Monthly Weather Review, vol. 114, 1986.
[67] Staniforth A. and Temperton C., ‘An Efficient Two-Time-Level Semi-Lagrangian Semi-Implicit Integration Scheme', Quarterly Journal of the Royal Meteorological Society, vol. 113, 1987.
[68] Bermejo R., ‘An Analysis of an Algorithm for the Galerkin-Characteristic Method ', Numerische Mathematik, vol. 60, 1991.
[69] Ryskin G. and Leal L.G., ‘Orthogonal Map-

OCR for page 153

Proceedings of the Sixth International Conference on Numerical Ship Hydrodynamics
ping', Journal of Computational Physics, vol. 50, 1983.
[70] Chang S.C. and Wu S.T., ‘On the Natural Frequencies of Standing Waves in a Canal of Arbitrary Shape', Journal of Applied Mathematics and Physics, vol. 23, 1972.
[71] Lamb H., ‘Hydrodynamics', Cambridge University Press, 1932.
[72] Aanesland V., ‘A Hybrid Model for Calculating Wave Making Resistance', Proceedings of the 5th International Conference on Numerical Ship Hydrodynamics, 1989.
[73] ‘Proceedings of the First Workshop on Ship Wave Resistance Computations ', David Taylor Naval Ship Research and Development Center, 1979.
[74] ‘Proceedings of the Second Workshop on Ship Wave Resistance Computations ', David Taylor Naval Ship Research and Development Center, 1983.