An investigation on HOBEM in evaluating ship wave of high speed displacement ship *

Xi Chen , Ren-chuan Zhu, Ya-lan Song, Ju Fan

1. Marine Design & Research Institute of China, Shanghai 200011, China

2. State Key Laboratory of Ocean Engineering, Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University,Shanghai 200240, China

Abstract: A practical numerical tool is developed to evaluate ship waves of high speed displacement ships on the basis of potential flow theory, in which high order boundary element method (HOBEM) based on biquadratic shape functions is applied to solve the boundary value problem. Since the sinkage and trim of ship at high speeds are notable, influences of ship attitude on wave drag are investigated and three kinds of models are used to evaluate them. To make the numerical approach highly efficient, an incomplete LU factorization preconditioner is adopted and incorporated with the restarted generalized minimal residual method GMRES (m) to solve the boundary integral equation. A corresponding Fortran code is developed and applied to evaluate ship waves of the Wigley hull and 4a model, a transom stern ship. Computations are performed for both monohulls and catamarans over a wide range of Froude numbers. Numerical issues including mesh convergence and computational efficiency are investigated at first. Computed results of the wave drag, sinkage and trim show generally good agreement with experimental data. Reasonable wave patterns are obtained and physical phenomena that wake angle where the largest waves occur, would become narrow at high speeds is also captured by the present computations. Numerical results indicate the proposed method would be accurate and efficient to evaluate resistance for hull design of high speed displacement ship.

Key words: High speed ship, ship waves, HOBEM, ILU-GMRES, sinkage and trim

Introduction

Resistance performance is usually the key factor in hull form design and line optimization of high speed ships. Over the past few years, great efforts have been made by researchers to study the problem of high speed ship advancing in water. Among these,Mizine et al.[1], Wang et al.[2], Wang and Lu[3]and Chen et al.[4]carried out experimental study on ship resistance of high speed monohull, catamaran and trimaran over various forms. With the instructions of experiments, numerical techniques have been developed a lot. Numerical method with high accuracy and efficiency is favorable in design stage.Computational fluid dynamics (CFD) is one of the numerical tools, which can usually provide accurate predictions[1,5-6]. However, due to its well-known disadvantage of low efficiency, CFD is hard to be widely used in engineering field especially for hull form design and bodyline optimization which require massive computations and lots of results.

An alternative choice is potential flow method. It has overall good accuracy and high efficiency.Potential flow theory can be classified into free surface Green's function method and Rankine source method. Neumann-Kelvin and Neumann-Michell(NM)[7-8]theory, in which Kelvin source Green's function is employed, are two classic free surface Green's function methods. Song et al.[9]derived analytical expressions of integral along line segment of wavelike disturbance in Kelvin source applied NK theory to predicted ship wave resistance. In addition,Kara and Vassalos[10]applied transient Green function to calculate wave-making resistance.

Rankine source method takes elementary solution of Laplace's equation as Green's function. Therefore,time cost in calculating complex Green's function is greatly shortened. And in the meanwhile, the problem of complex Green's function oscillating near free surface shall be successfully avoided. Another difference with Green's function method is that by applying Rankine source method, sources are required to be distributed on all domain boundaries. Accordingly, total number of mesh will increase and time cost in solving boundary integral equation increases.Since Dawson first applied Rankine source method to simulate ship waves, many researchers like Raven[11],Tarafder and Suzuki[12], Peng et al.[13]have done researches to improve numerical scheme of Dawson's method. Wyatt[14]came up with a nonlinear method to evaluate ship waves. In the above listed studies,constant panel approach was adopted. It assumes that velocity potential is discontinuous over boundaries and discretization on free surface would bring large numerical damping and dispersion. Therefore,applying high order boundary element method(HOBEM) is meaningful when Rankine source method is used. Chen et al.[15]applied a time domain HOBEM method to calculate Wigley hull wave drags.Gao and Zou[16]used NURBS based HOBEM to compute wave drags of submerged ellipsoid and showed that HOBEM has high accuracy. Although HOBEM was employed, Chen et al.[15], Gao and Zou[16]did not apply the methods to realistic ships with more complex hull forms for further investigations of their applicability. Chen et al.[17]proposed a nonlinear potential flow method incorporated with HOBEM to evaluate ship waves of different hulls.

The above mentioned researches are mainly carried out for ships under condition of low or moderate speed (Froude number smaller than 0.5).With respect to ship waves of high speed ships,Moraes et al.[18]used both Rankine source method and slender body theory to investigate wave drag of high speed catamaran of simple hull forms without transom sterns. Gao[19]studied the implementation of transom conditions on high speed ship and used HOBEM to compute wave drags of the NPL hull, while the sinkage and trim which are important for high speed ships were not taken into account. Tarafder and Suzuki[20]computed wave making resistance of the Wigley catamarans in sufficient high speeds by a weakly nonlinear Rankine source method on the basis of second order perturbation and the sinkage and trim were also neglected in computation. Wang and Lu[3]applied a fully nonlinear Rankine source method with constant panel discretization to compute trimaran resistance with consideration of ship attitude.Numerical issues of their method like convergence and efficiency were not discussed in detail.

The present study aims at a practical numerical method to evaluate ship waves of high speed displacement ships with and without transom stern.Rankine source method in conjunction with HOBEM is used in this study. The linearized free surface condition with higher stability is adopted in computation on consideration that nonlinear computations could not always converge at high Froude numbers.The sinkage and trim that are significant to high speed ships are also taken into account. Three kinds of methods, including iteration approaches, are used to evaluate the sinkage and trim. Nonlinear hydrostatic forces and moments computed by pressure integral on exact wetted ship hull are introduced in the present iteration methods. It differs from the linear method, by multiplying the ship displacements with hydrostatic coefficients, adopted in previous studies like Wang and Lu[3]and is more accurate for non-wall sided ships.

To improve computational efficiency, an iterative matrix solver, the restarted generalized minimal residual algorithm (GMRES (m), m denotes restarted number), is imposed to solve the linear system. And an incomplete LU factorization (ILU) preconditioner specific to the present problem is proposed to overcome the difficulty of poor convergence. A Fortran code is developed and the mathematical Wigley hull and 4a model, which is of round bilge type with transom stern, are used for validation and verification studies. Both the monohull and catamaran are considered. Wave drags, sinkage and trim of ships,and wave patterns are computed and compared with related model test data followed by detailed discussions. Results show the present method has generally good accuracy and high efficiency, and the numerical scheme is robust. Detailed discussions on physical phenomenon of ship waves of high speed ships are also carried out.

1. Mathematical and numerical formulation of ship waves

1.1 Mathematical model

1.1.1 Governing equation

For the ship wave problem, considering a ship advancing in calm water with constant speed u and a Cartesian coordinate o-xyz translating with ship is defined. As shown in Fig. 1, the origin o is on midship and xoy plane is fixed on still water level with x and y axes directing towards bow and portside respectively.

According to the potential φ flow theory, flow velocity potential satisfies the Laplace equation.

Fig. 1 Coordinate system

where total velocity potential can be decomposed into uniform stream potential -ux and disturbance potential φ. The boundary value problem (BVP) about φ can be constructed with the following boundary conditions

where subscripts xx , z and n represent partial derivatives.

The free surface condition is the Kelvin-Michell linearized boundary condition at the free surface.Dawson's linearization is not adopted here for its slow ship hypothesis. Radiation condition signifies no free surface waves propagate upstream. In computation,free surface is truncated at limited range. So in order to satisfy the radiation condition, Dawson's upstream differential scheme incorporated with technique of shifting source panels downstream a distance, usually 0.25 times of longitudinal panel length is employed in computations.

For ships with transom stern, specific conditions should be imposed. Raven[11]classified flows occurring after transom stern into three types: (1) regular stern flow, (2) transom stern flow, (3) dead water type of flow. Dead water type of flow, which is turbulent flow and is hard to be modeled with potential flow theory, exists at low speed. With regard to high speed ships, only the other two types of flow would occur.Regular stern flow is the same as flow after cruiser stern, so no extra formulas are needed in calculation.For transom stern flow type, fluid would leave hull along the lower edge of transom. To simulate this type of flow, the following conditions are adopted

1.1.2 Wave drag, sinkage and trim

Hydrodynamic forces and moments are computed by pressure integrals on ship wetted surface.

where Sbdenotes the wetted hull, Fjrepresents hydrodynamic force and moment in j-th direction,n is unit normal vector which points outside of domain boundaries, p is hydrodynamic pressure and is determined by Bernoulli's equation, which is

Wave drag coefficients are non-dimensional quantities of longitudinal forces.

For transom stern ships, the hydrostatic correction on wave drag as defined in the following equation should be considered in order to compensate the absence of hydrostatic pressure on transom.

where the integral is performed along transom edge.

Fig. 2 Three models used in sinkage and trim computation

In this paper, three different kinds of models are used to compute sinkage and trim of ships. As illustrated in Fig. 2, free 1, free 3 models employ stationary iterative approach, while free 2 model is a linear method in which ship sinkage and trim are obtained by solving the linear ship motion Eq. (9).

where cijare restoring coefficients, h and θ denote sinkage and trim respectively.

In free 1 and free 3 models, sinkage and trim of ships are solved iteratively by the following equations

where Δh and θΔ are increments of sinkage and trim,and Sbsignify the mean and exact wetted ship surface respectively.

Models 1, 3 are different. By using free 1 model,hydrodynamic forces would be updated in each iteration and hydrostatic forces are computed with ship wetted hull under still water level. By contrast,hydrodynamic forces are invariable in the free 3 model and nonlinear hydrostatic forces are considered by pressure integration on exact wetted hull under wavy free surface. The boundary value problem So the computing efficiency is free 2＞free 3＞free 1. The accuracy of different models are investigated in section 2.

1.2 Numerical method

1.2.1 Boundary integral equation

By using rankine source method, velocity potential and its spatial derivatives at field point p are expressed as

where p and q denote field and source point respectively, σ is source strength,=1/r is Rankine source Green's function.

Substituting Eq. (11) into hull boundary condition and free-surface condition as in Eq. (2),boundary integral equations (BIE) can be derived

where Nband Nfrepresent nodes number on hull and free surface respectively,bA and Afare influence coefficients. Note that transom flow conditions should be added in BIE for ships with transom stern.

1.2.2 High order boundary element method

HOBEM with 9-nodes isoparametric high-order panel as shown in Fig. 3 is adopted in the numerical computations. By applying this kind of high order panel, distributions of source strength and velocity potential are continuous over boundaries. And less curved panels can discretize curved ship hull better than 4 nodes constant panels. Therefore, HOBEM shall have higher accuracy and efficiency than constant panel method.

Position vector [ x , y , z] and source strength σ in panel can be approximated by sum of its nodal values weighted by biquadratic shape functions

where j is local nodal index.Applying Eq. (13) into free surface and body boundary conditions and choosing nodes of elements as collection points, influence coefficients A of BIE can be expressed by the following equations

where subscripts i and j are row and column number of coefficient matrix,iα stands for solid angle of i-th node, δ is the Kronecker delta function, Ine( e, k) represents global index of k-th node in e-th panel.

In Eq. (14), integrals of G ( p , q) and ? ( ,pG pon elements are generally computed by Gauss-Legendre algorithm. While if a source point coincides with a field point, singularities of influence coefficients would occur as integral kernels approach to infinite. In order to eliminate these singularities, the classical raised panel method is applied on freesurface. Refer to Raven's study[11], source panels are moved with a distance of 0.8 times of characteristic panel length above free surface in computations. On ship hull, singular integral ofis calculated by triangular polar-coordinate transformation[21]. And another singular kernelapproaches to zero on smooth hull boundary. As a result, diagonal values of influence coefficient sub-matrix with respect to ship hull are equivalent to solid angles that can be evaluated by formulas in the Ref. [21].

1.2.3 ILU Preconditioned GMRES (m) solver

Boundary element method would yield an un-symmetry and dense linear system. Since the direct solver with computation amount N3is inefficient if there are thousands of unknowns, the restarted generalized minimum residual method GMRES (m),which is applicable to un-symmetric linear system[22],is adopted as matrix solver in this paper. However,numerical tests indicate the iterative solver is hard to converge with respect to the present problem as the coefficient matrix is non-diagonally dominant. By this reason, a preconditioner is developed and imposed to improve matrix property. Algorithm of the preconditioned GMRES (m) solver can be found in numerous references such as Saad[22]. As simple Jacobi or SSOR preconditioner is futile in this work, a more sophisticated incomplete LU factorization (ILU) preconditioner is indispensable. The ILU-preconditioner adopted in this study uses a single threshold ε to filter out contributions of far fields. As the decomposed L and U matrixes are sparse, the compressed sparse row (CSR) storage format is employed for the precondition matrix to save storage and improve efficiency.

2. Numerical results and discussions

A Fortran program is developed for validation and verification studies. The Wigley and 4a model monohulls and catamarans are used as numerical examples. Principal dimensions of both demihulls are presented in Table 1. Note that ship models are non-dimensioned by the ship length in calculations.Computed results include ship wave drag, sinkage,trim and wave patterns over a wide range of Froude numbers (0.1-1.0).

The Wigley demihull is defined by the following formula

The 4a model is derived from NPL series and it is round bilge form high speed ship with transom stern.Body plan of the 4a model is shown in Fig. 4.

Fig. 4 Body lines of the 4a model

Fig. 5 (Color online) Panel arrangment of free surface

Panel arrangement of half of free surface is shown in Fig. 5, where Nlstands for number of nodes (With HOBEM, the quantity of mesh is measured by number of nodes) distributed in the longitudinal, Nwoand Nwiarenodesofouterandinner zones respectively, in lateral direction. Free surface mesh is the body fitted type and is extended to about 2.0 times of ship length downstream, 0.8 times of ship length upstream and one ship length portside.

Table 1 Comparison of computational times

Table 2 Mesh scheme of the Wigley monohull and catamarans

2.1 Numerical issues

Mesh convergence as well as efficiency and effectiveness of the matrix solver are studied in this section.

For high speed hulls, it is not suitable to adopt the criterion, number of mesh per transverse wave length, to discretize computational domain boundary.That is because transverse wave lengthis much larger than ship length at high speeds, and applying the criterion would result in oversize elements that cannot capture diverging waves. So number of nodes distributed per ship length(nodes/ L) is used as meshing standard in this work.Figure 6 depicts convergence of the Wigley monohull wave drags at =0.20-1.00. It is observed that the results can converge if nodes/L is greater than 31.In the following computations, mesh discretization of all cases will follow the standard of nodes/ L＞31.

Fig. 6 Wigley wave drag computations with different grid properties

2.2 Wave drag, sinkage and trim of the Wigley hull

Mesh properties in the computations are presented in Table 2.

2.2.1 Wigley monohull Wave drags, sinkage and trim as well as wave profiles of the Wigley monohull computed by the present method in comparison with model test data[7,17]are shown in Figs. 7-9.

Fig. 7 Wigley monohull wave drags of present computations and experiments in free and fixed conditions

Numerical results of the wave drags in both fixed and free conditions agree quite well with experiments as seen in Fig. 7. Evident effects of the sinkage and trim on wave drags are observed at. Free model results are much larger than that of fixed model.At, ship attitude would bring an increase of about 32% of the wave drag.wC increases with Froude numbers at the regime ofand humps and hollows on the Cwcurve, which indicates the constructive and destructive wave interferences,are obvious in the figure. A sharp growth is occurred nearon account of the approach of the transverse wave and ship lengths. As speeds go higher,would reach up to the main crest at ≈0.50 Fr .And after that, the curve experiences a smooth decrease with Fr, as transverse wave length becomes much larger than ship length and the interference effects for transverse waves vanishes.

Figures 8, 9 depict the sinkage and trim evaluated by the three models abovementioned. Numerical results are in overall good agreement with model test data except for slight deviations of sinkage at Fr=0.35-0.40. It indicates fluid viscosity has negligible effect on sinkage and trim. Slight differences of the results of the three models are only observed at. Besides, trim angles are close to zero at. After the transition point, 0.35 Fr ≈ , trim angle increases dramatically since transverse wave length gradually approaches to ship length. The trends of trim angle curves may explain why apparent discrepancies between free model and fixed model resistances are observed at.

Fig. 8 Wigley monohull sinkages of present computations and experiments

Fig. 10 Wigley monohull wave profiles of present computations and experiments at Fr=0.14

The computed wave profile are in good accordance with experimental data except for bow waves as plotted in Fig. 10. The underestimation of free surface elevation near bow may be due to wave spray that cannot be modeled in potential flow frame.

2.2.2 Wigley catamarans s/ L = 0.2, 0.4 and 1.0

Wave resistance, sinkage and trim of the Wigley catamarans with separation ratios, s/ L = 0.2, 0.4 and 1.0, are computed and results are shown in Figs. 11,12. The present computation agree well with results of Moraes et al.[18], which are computed by Shipflow software (a 3-D panel method), as seen from Fig. 12.The figure also shows a severe increase of Cwwith the reduction of separation ratio at =0.40-0.70.The peak value of wave drag coefficient of catamaran at s/ L=0.2 is about 2.5 times higher than that of monohull. Wave interferences between demihulls seem to be insignificant at higher speeds, i.e.,. This is because the wake angle of wave system decreases and thus significant wave interferences merely exist in far field behind stern.

Fig. 11 Wigley catamaran ( s/ L =0.2) wave drags of present computations in free and fixed conditions

Fig. 12 Wigley monohull and catamaran ( s/L = 0.2, 0.4 and 1.0) wave drags of present and Moraes et al.'s computaions[18] in free comdition

2.3 Wave drag, sinkage and trim of the 4a model

The second numerical example considered is the 4a model that is of round bilge type with transom stern. Towing tank experiments are conducted by Molland et al. Experimental data can be found in report ”Resistance experiments on a systematic series of high speed displacement catamaran forms: variation of length-displacement ratio and breadth-drought ratio, 1994”. Measured wave drags of 4a model in free condition are obtained by 1+k method. Mesh properties of the numerical computations are in Table 3.

Table 3 Mesh scheme of 4a model monohull and catamarans

2.3.1 4a model monohull

Figures 13, 14 depict sinkage and trim predicted by different methods in comparison with the measured data. Results of iterative computations, by the free 1, 3 models, generally agree well with the experimental data and have overall higher accuracy than the linear calculation. Some discrepancies between the measured and computed results exist at =0.60-1.00 Fr ,since ship model is in semi-planning state in this speed regime and nonlinear effects of flow field actually have impact on ship attitude. Results of the free 1, free 3 models are very close. Considering computational efficiency, the free 3 model is more approval.

Fig. 13 4a model monohull sinkages of present computations and experiments

Fig. 14 4a model monohull trims of present computations and experiments

Fig. 15 4a model monohull wave drags of present computations and experiment in free and fixed conditions

Fig. 16. 4a model monohull and catamaran ( s/ L =0.2, 0.3, 0.4 and 0.5) wave drags of present computaions and experiments in free comdition

2.3.2 4a model catamarans s/ L= 0.2, 0.3, 0.4 and 0.5

Figure 16 shows excellent agreements between calculated and measured wave drags of the 4a model catamarans at separation ratios 0.2-0.5. The computed sinkage and trim are also in fair accordance with experimental results as shown in Figs. 17, 18.Demihull wave interferences have a distinct tendency to arise the magnitudes of 4a model catamaran wave drags at0.40-0.70. And wave interferences gra- dually vanishes with the increase of speeds.Possible reasons for the physical phenomenon have been stated in Section 2.2.

Fig. 17 4a model monohull and catamaran ( s/ L =0.2, 0.3, 0.4 and 0.5) sinkages of present computaions and experiments

Fig. 18 4a model monohull and catamaran ( s/ L =0.2, 0.3, 0.4 and 0.5) trims of present computaions and experiments

2.4 Wave patterns of high speed ships

Figure 19 shows a comparison of the Wigley monohull wave patterns at Fr= 0.40, 0.80. Transverse wave length is equivalent to ship length at,and both transverse and diverging waves are easily observed in upper half of the figure. By contrast, attransverse wave length are far more larger than ship length, and diverging waves are dominant in the ship wave system as seen in the lower part of the figure. Besides, the result shows the wake angle φmax,where the maximum wave amplitude occurs, decreases at high speed. This observation has been studied recently by Darmon et al.[23], Miao and Liu[24],analytically with Gaussian pressure and an arbitrary surface pressure disturbance, and by Zhang et al.[25]with general hulls of a continuous source distribution.In Fig. 19, the wake angle=9.14°, labeled as the red dash line, is computed with method of Zhang et al.[25]and the present evaluation conforms to the result.

Fig. 19 (Color online) Wigley monohull wave patterns of present computations at Fr = 0.4, 0.8, as well as Kelvin angle k φ and wake angle max φ = 9.14° evaluated by method of Zhang et al.[25]

To further investigate the physical phenomenon,wave patterns of the Wigley and 4a model monohull atare computed and depicted in Fig. 20.Wake angles φmaxrepresented by the red dash lines,which are 6.43°, 6.33° for the Wigley and 4a model repectively, are solved by method of Zhang et al.[25].In their method, the stationary-phase approximation is applied to determine relations of ray angle and amplitudes of divergent waves. The present computed highest waves are in accordance with the wake angle obtained by Zhang's method as shown in the figure.The purple shot dash line in the figure represents=12.85° of Gaussian pressure field computed by Darmon et al.[23]. It is found that φmaxof real ships are narrower than that of pressure field, while results of different hulls just have slight difference.

Fig. 20 (Color online) Wigley and 4a model monohull wave patterns of present computations at Fr =1.00 , as well as Kelvin angle k φ and wake angle max φ =43°,6.33° of the Wigley and 4a model evaluated by method of Zhang et al.[25], and ,ma.85 G φ ° of Darmon et al.'s solution[23]

3. Conclusions

The present study applies HOBEM Rankine source method to simulate ship waves of high speed displacement ships. The sinkage and trim that are notable at high speeds are computed with three different methods and are taken into account in computation of wave drags. An ILU preconditioned GMRES (m) method is employed to solve the boundary integral equation. Numerical investigations are carried out for the Wigley and 4a model monohulls and catamarans. Ship waves are well evaluated by the present method, the following conclusions can be drawn by discussions:

(1) The computed wave drag, sinkage and trim of different monohulls and catamarans advancing with high speeds in comparison with related model test data indicates the numerical method are generally satisfactory. Convergence study shows HOBEM is stable and robust to solve ship wave problems. The developed iteration method, the free 3 model, is more accurate than linear method to compute sinkage and trim of realistic ships like the 4a model. Moreover, the proposed ILU preconditioned iterative matrix solver can effectively improve computing efficiency.

(2) For high speed ships, the sinkage and trim have significant influences on wave drags for both monohulls and catamarans. For example, ship attitude would bring an increase of about 30%-40% of the wave drag at. Differing from conventional ships, wave drag coefficients of high speed ships diminish smoothly with speed if Fr is greater than about 0.50, as divergent waves become dominant and the interactions with transverse waves vanishes. For catamarans, wave interference between demihulls would markedly magnify magnitudes of wave drag coefficients, the sinkage and trim at the speed regime of =0.40-0.70. In much higher Froude numbers,Fr is larger than about 0.70, the interference effects shall become unapparent and be moved downstream due to the reduction of wake angle φmax.

(3) Numerical results of wave patterns of both hulls are reasonable at different speeds. In high speed regime, physical phenomena of the decrease of wake anglewith Fr is observed by the present computation. The numerical results are also validated by some other recent studies. It is also found thatof actual ships vary little with hull forms and is narrower than that of Gaussian pressure field.