VIEWS: 42 PAGES: 36 POSTED ON: 5/8/2011
COMPARISON BETWEEN NUMERICAL COMPUTATIONS AND EXPERIMENTS FOR SEAKEEPING ON SHIP'S MODELS WITH FORWARD SPEED C.Maury*, G.Delhommeau*, M.Ba**, J.P.Boin*** and M.Guilbaud*** * L.M.F. (UMR CNRS n°6598), Ecole Centrale de Nantes, BP 92101, 44321 Nantes Cedex 3, France; Tel : 33-2-40-37-25-95 Fax : 33-2-40-74-74-06 ** L.E.A. (UMR CNRS n°6609), ENSMA, 1 rue C. Ader, BP 40109, 86960 Futuroscope Chasseneuil Cedex, France; Tel : 33-5-49-49-80-87 Fax : 33-5-49-49-80-88 *** L.E.A.-CEAT (UMR CNRS n°6609), Université de Poitiers, 43 rue de l'Aérodrome, 86036 Poitiers Cedex, France; Tel : 33-5-49-53-70-27 Fax : 33-5-49-53-70-01 Corresponding author: Guilbaud Michel, LEA-CEAT, 43 rue de l’Aérodrome - 86036 Poitiers Cedex, France ; tel. 33-5-49-53-70-27;fax 33-5-49-53-70-01; e-mail:guilbaud@univ- poitiers.fr 24/10/01 Abstract The recent progress of computers and numerical algorithms enables today to use the translating and pulsating Green function in panel methods for seakeeping calculations. Two different methods of calculations of this function, its derivatives and their integrations over panels or segments are briefly presented and have been introduced into the seakeeping codes Aquaplus developed at Ecole Centrale de Nantes and Poseidon at CEAT-LEA. Both codes interchange the Fourier and boundary integrals on panels or waterline segments, the last part being performed analytically. These methods have been used to compute flows around Wigley or series 60 model ships. To check the numerical results, an experimental set-up has been developed at the CEAT, which measure forces and moments on model in forced harmonic oscillations of pitch or heave. Tests have been performed in the recirculating water channel of Ecole Centrale de Nantes on two L=1.2m series 60 models of CB=0.6 and 0.8 block coefficients. Unsteady wave patterns have been recorded using a resistive wave probe. The experimental results are compared with the numerical ones. 1 Introduction The numerical difficulties to solve numerically the seakeeping problem with forward speed lead to the necessity to check numerical results with tests. Few experimental data concerning this problem are available and in some cases, not appropriate for accurate comparison. Moreover, the comparison of measured and calculated global forces or motions of ship running in waves is not an efficient check for the quality of numerical results because they involve hydrodynamics and mechanics, (Okhusu & Wen 1996, Okhusu 1998). So comparison with calculations does not enhanced the progress of the numerical methods and it becomes necessary to perform new tests with local measurements (pressure distribution, wave pattern) to improve the knowledge about the accuracy of the calculation methods. Some test results can be found in (Okhusu 1998), but also in (Iwashita et al 1993) concerning OHS (a bulk carrier's hull form of CB=0.8 and L/B=5.48 where L is the length and B the width), Series 60 or VLCC (very large crude carrier) ships. Two different techniques can be used in tests, free models in waves where forces and moments but also hull motions have to be measured or forced motion technique where the motion amplitude is known accurately and only the phase lag reference can be subject to errors. Panel methods are very efficient for seakeeping calculations. For the problem without forward speed, Green function is easy to compute. For the problem with forward speed, a lot of numerical and theoretical difficulties appear. So, for a long time, only frequency encounter approximation limited to thin ships (with the zero speed diffraction-radiation Green function) or Rankine methods available for high values of Brard's parameter τ=ωeU/g (ωe encounter circular frequency and U forward speed) have been used. Recently, using either steepest descent method (Iwashita & Okhusu 1989), (Brument & Delhommeau 1997), (Brument 1998), (Maury 2000), Super Green function method (Chen & Noblesse 1998) or Simpson adaptative method (Ba & Guilbaud 1995), (Nontakaew et al 1997), progress in the calculations of the Green's function have been performed. It is now possible to develop seakeeping codes using this Green function on small workstation with moderate computational time and good accuracy, see (Du et al 1999,2000), (Boin et al 2000). In spite of the numerical difficulties of its computation, the advantages of this Green function are an automatic fulfilment of the linearized free-surface and radiation conditions, the absence of mesh on the free-surface avoiding reflection on the boundary and the filtering of smaller wave lengths and the use of method, it is necessary to compute accurately the integrals of the Green function and its derivatives on elementary panel or waterline segment. Numerical quadrature 2 is not able to give correct results for a panel and a collocation point close to the free surface, (Guilbaud et al 2000). So it is needed to integrate analytically on panels or segments following formulas given by (Iwashita 1992) for the steepest descent method or (Bougis 1981) with the formulation of (Guevel & Bougis 1982). We present here two numerical methods to perform seakeeping calculations computing accurately this Green function and its derivatives and their summation on panels or segments. The aim is to obtain, in all cases and especially for high forward speed, the Green function and its derivatives with a given accuracy. The first method, Aquaplus, developed at Ecole Centrale de Nantes is founded on Bessho's formulation (Bessho 1977) and uses the steepest descent method. The second method, developed by the CEAT, uses the formulation of (Guevel and Bougis 1982) with computation of single integrals by Simpson adaptative method. We present also an experimental study using the forced motion technique in order to compare the results with those of the two numerical methods already mentioned. Forces, moments and also wave patterns have been measured in circulating water channel on Series 60 ship's models with length L=1.2m in forced oscillations of heave or pitch, using a development of the experimental set-up used in (Guyot 1995) or (Guyot & Guilbaud 1995) where a L=0.6m model have been used. The experimental set-up is a development of the one used in (Delhommeau et al 1992). Calculations of the Green function Computation of the Green function with steepest descent method Point source : Let us consider a source point at P ' (x ', y ', z ') and a field point at P (x , y, z ) . The source is moving at constant speed U and pulsating with an encounter pulsation ωe . The Green function satisfies the linearized free surface condition and is expressed under the form (1) derived in (Bessho 1977), which presents the main advantage to be written as a single integral over the complex variable θ, allowing to integrate over any available paths in the complex plane. Using this expression, one can evaluate at the same time the Green function and its derivatives of first or higher order. Development of this method, initially proposed in (Iwashita 1992), has been done at Ecole Centrale de Nantes (Brument & Delhommeau 1997, Brument 1998), or more recently (Maury 2000). g is the Froude dependent (or Kelvin) part of the Green function. 3 1 1 − 1 + g (P, P ') = 1 1 − 1 − iK0 θ2 (P ,P ') k2ek2ξ − sgnc.k1ek1ξ G (P, P ') = 4π R R1 4π R R1 2π ∫ θ1 1 + 4τ cos θ .d θ (1) R1 2 2 2 where = (x − x ') + (y − y ') + (z ± z ') R X = K 0 (x − x ') g k1 1 + 2τ cos θ ± 1 + 4τ cos θ K0 = = U2 Y = K0 y − y ' and k2 2 cos2 θ U ωe τ= Z = K 0 (z + z ') g ξ = Z + i (X cos θ + Y sin θ ) π θ2 (P, P ') = ϕ (P, P ') − − iε (P, P ') π + arccos 1 1 2 − if τ > 4τ X θ1 = 4 and ϕ (P, P ') = arccos 1 1 X 2 +Y 2 −π − i arg ch if τ < 4τ 4 Z ε (P, P ') = arg sh X 2 +Y 2 −1 if t < 0 sgn c = sign [ Re (cos θ )] , with sign (t ) = 0 if t =0 sgn s = sign [ Re (sin θ )] 1 if t > 0 Due to their behaviours, terms in k1 and k2 are computed differently. Terms in k2 are computed directly in the θ -space, whereas the steepest descent method is needed for terms in k1 due to strong oscillations in the vicinity of ±π / 2 . To apply the steepest descent method we separate the integral over θ in different parts following the sign of terms sgnc and sgn s : −π 2 π ∫θ f1 (θ).dθ if X > 0 θ2 (P,P ') − 2(P,P ') I =∫ −sgnc.f1 (θ).dθ = ∫ f1 (θ).dθ + 2 θ1 θ1 − π π π ∫ 2 f1 (θ).dθ − ∫ 2 f1 (θ).dθ + ∫ 2 f1 (θ).dθ 0 if X < 0 0 θ2 (P,P ') k1e k1ξ with f1(θ) = . The steepest descent method is applied in a space M defined 1 + 4τ cos θ by : M = sgn c.k1.cos θ = sgn c.m . Rewriting integrals in respect to variable M , we obtain (2) for the function and its derivatives: I M ( ) if X > 0 ( ) θ2 I = I Mθ + (2) 1 I (M 0 ) sgn s =−1 − I (M 0 ) sgn s =1 ( ) + I Mθ 2 if X < 0 4 1 i.m ∞ 2 e φ(M ) with I (β ) = sgn c.sgn s ∫ m dM β i.sgn s. (m − τ )2 1 − 2 ( m m − τ )2 1 − (m − τ )2 (m − τ )2 2 m and φ (M ) = (m − τ )2 Z + i.m.X + i.sgnsY (m − τ )2 1 − . m − τ )2 . ( It must be noticed that the term (m − τ )2 is out of the square root avoiding the existence of a branch cut in the M -plane, which simplifies quite significantly the research of steepest descent lines. Integration is done following a path on which the imaginary part of φ is constant, suppressing the oscillations of integrands, in a direction where the real part decreases. Steepest descent lines do not exist in all cases because they must be entirely included in a θ -space area where sgnc and sgn s are constant. Moreover, the steepest descent line is determined numerically with a given accuracy. In the case of a failure in the research of a steepest descent line, the starting point β is moved to β ' in order to find another line. A complementary integral between β and β ' is then introduced and computed in θ -space. β ' must be chosen carefully to avoid large oscillations of the complementary integral. Panel and segment source : The solution of 3-D seakeeping problem with forward speed by the boundary element method needs to compute an integral of a source distribution over panels on the body or segments on the waterline. Thus the computation of the Green function associated with a panel or a segment of constant source density is developed. The analytical quadrature of Kelvin part g of (1) over a panel or a segment (Iwashita 1992) gives us expressions I S = ∫∫ gds for a panel ∆S with N vertices and I L = ∫ gdl for a straight line: δC 1 k2ξl 1 e − sgnc ek1ξl i N θ2 (P ,Ql ) 2Sl k2 k1 i 1 N S IS = − 2πK 0 ∑∫ θ1 ψl (θ ) 1 + 4τ cos θ dθ + πK0τ 2 ∫ ∑ ψ (θ) 0 l νl (η)d η l =1 l =2 l θ=θ2 (P,Qη ) i β (Y1 −Y2 ) N θ2 (P ,Ql ) 1 e k2ξl − sgnc.e k1ξl IL = ∑∫ dθ (3) 2π l =1 θ1 (ξl +1 − ξl ) 1 + 4τ cos θ with :Ql vertices of panel (l = 1, N ) and Sl area of triangle (Ql ,Ql +1,Ql −1 ) 5 ( Xl = K 0 x P − xQ l ) Yl = K 0 yP − yQ (with the sign of (yP − yQ ) constant for l = 1, N l l ( Zl = K 0 z P + zQ l ) ψl (θ ) = (ξl +1 − ξl )(ξl −1 − ξl ) and ξl = ξ (P,Ql ) X ηYl '− Xl 'Yη Zl ' (X η + Yη2 ) − Z η (X η Xl '+ Yl 'Yη ) 2 νl (η ) = +I 2 X η + Yη2 (X η2 + Yη2 ) (X η2 + Yη2 + Z η2 ) X η = X 1 + η (X l − X 1 ) X l ' = Xl − X 1 where Yη = Y1 + η (Yl −Y1 ) and Yl ' = Yl −Y1 Z η = Z 1 + η (Zl − Z1 ) Zl ' = Zl − Z 1 Expressions of the derivatives may be found in (Maury 2000) or (Ba et al 2001a,b). Here again, the steepest descent method is useful to evaluate efficiently terms in k1 . For this reason and the fact that θ2 varies with Ql , computation of integrals corresponding to each vertex l , are done separately. Due to this decomposition of the sum, singular points θl appears in integrands when ψl (θ ) equals to zero. The values of these singular points are known and evaluated by the (4) below: Xel = XQ − XQ Ye ± Xe 2 + Ye 2 + Ze 2 (l +1) l l with θl = 2 arctan l l l Yel = YQ −YQ (4) Xel − i.Zel (l +1) l Zel = ZQ − ZQ (l +1) l Paths of integration must avoid these singular points, and the contribution of poles must be taken into account. It is necessary to evaluate this contribution in θl , referred to vertices Ql and Ql +1 , when the paths used for both integrals ( l and l + 1 ) surround the singular point considered. When the influenced point is far from the panel (or the segment), we can optimise the evaluation of the Green function using an averaged path of integration to integrate directly the whole sum from θ1 to θc . Only one integral has to be evaluated and there is no singular point to take care. Nevertheless, θ2(l ) being a function of the vertices Ql , we have to complete this expression by a sum of integrals from θc to each θ2(l ) as written in (5). Those integrals are done either in θ -space or with steepest descent method: 6 θc N N θ2 (P ,Ql ) I S (or I L ) = ∫ θ1 ∑ − sgnc .f1,l (θ ).d θ + ∑ ∫ l =1 l =1 θc − sgnc .f1,l (θ ).d θ (5) 1 k1ξl e −i 2S l k1 i β (Y1 − Y2 ) e k1ξl where f1,l (θ ) = or f1,l (θ ) = 2πK 0 ψl (θ ) 1 + 4 τ cos θ 2π (ξl +1 − ξl ) 1 + 4 τ cos θ Concerning the terms in k2 , we can easily use a common path for each integral associated with each vertex of the panel. Integrands have no oscillatory behaviour so we have no constraint on the path of integration. Expression (3) is decomposed as in (5). Here again, we have to complete the integral by a sum of integrals from θc to each θ2(l ) . It had been shown, (Maury 2000) that the difficulty to calculate increase with the order of derivation but decrease from the function itself to the segment integral and to the panel one. This is due to the fact that as the power in k1 is higher the integrands decreases more rapidly. Computations of the Green function with a Simpson Adaptative method Point source The Green function used is given in (Guevel & Bougis 1982) and developed in (Ba & Guilbaud 1995) with the dimensionless parameters ω = ωe 0 /g , F =U / g 0 , τ = ωF Kc = 4ω 2 where 0 is a reference length. It can be written as the sum of 3 functions G = G0 + g = G 0 + G1 + G2 , defined by: 1 1 1 G 0 (P , P ') = R R− , (6) 0 1 π 2 K g (K ξ ) + g (K ξ ') − K g (K ξ ) + g (K ξ ') 1 1 1 2 1 d θ , (7) ∫ 1 1 1 2 1 2 G 1 (P , P ') = π 0 0 1 + 4 τ cos θ θc' Z g (Z ξ ) + g 2 (Z 3 ξ ') − Z 4 g 2 (Z 4 ξ ) + g 2 (Z 4 ξ ') ∫ −i 3 2 3 dθ 0 4 τ cos θ − 1 θc −αc g (Z ξ ) + g (Z ξ ') − Z g (Z ξ ) + g (Z ξ ') Z3 1 3 + ∫ −i 1 3 4 3 4 3 4 d θ θ' 4 τ cos θ − 1 (8) 1 c G 2 (P , P ') = ( π 0 2πK c (1 − i ) e Kc ξc + e Kc ξc ) ' − αc τ sin θc π2 g (K ξ ) + g (K ξ ') − K g (K ξ ) + g (K ξ ') K3 3 3 + ∫ 3 3 4 1 4 1 4 dθ θc + αc 1 − 4 τ cos θ 7 0 ξ = z + z '+ i [(x − x ') cos θ + (y − y ') sin θ ] where with ξc=ξ(θ=θc), ξ'c=ξ'(θ=θc), ξ ' = z + z '+ i (x − x ') cos θ − (y − y ') sin θ 0 The integration limits are given by: if τ < 1/ 4 θc = θ 'c = αc = 0, if 1/ 4 < τ < 1/ 2 θ = arcos (1/4τ ), θ ' = 0, α = 10−6 , c c c if τ > 1/ 2 θc = arcos (1/4τ ), θ c = arcos (1/2τ ), αc = 10−6 ' The functions Kj(θ) and Zj(θ) are given by : 1 + 2τεj cos θ + (−1)j 1 + 4τ cos θ Kj = for j = 1 to 4, εj = +1 if j = 1,2; εj = −1 if j = 3, 4 2F 2 cos2 θ τ 1 − 2τ cos θ + (−1)j 4τ cos θ − 1 K5 = ; Zj = for j = 3, 4 F 2 cos θ 2F 2 cos2 θ The functions gi(ξ) are defined by : g1 (ξ ) = e ξ ε1 (ξ ) if 0<arg (ξ ) < 2π, g 2 (ξ ) = e ξ E1 (ξ ) if − π<arg (ξ ) < π; g 3 (ξ ) = e ξ ε1 (ξ ) + 2i π if 0<arg (ξ ) < 2π, with: ε1 (ξ ) = E1 (ξ ), if ℑ (ξ ) ≥ 0, or ε1 (ξ ) = E1 (ξ ) − 2i π, if ℑ (ξ ) < 0; where ℑ corresponds to the imaginary part and E1 is the complex-valued exponential integral function: ∞ ∞ e−t e −ξ t E1 (ξ ) = ∫ dt if − π < arg (ξ ) < π ; E1 (ξ ) = ∫ dt if ℜ (ξ ) > 0 ξ t 1 t The main difficulties concerns the Fourier integrals (G1 and G2) and the singular behaviour close to τ=1/4; furthermore, the integrands become very oscillating when z and z' are both close to 0. For the numerical calculations, we have replaced the 4th order Runge-Kutta method used in (Ba & Guilbaud 1995) with a variable step by an adaptative quadrature method, (Nontakaew et al 1997). The integration is done using a Simpson 1/3 method with five points to avoid to calculate more that once each value of the integrand needed. (Guttman 1983) has proposed an expression to relate the absolute error ε for an interval with length ln to the global error εtot for the total path with length L, ε=(15εtotln)/L. This method presents various advan- tages as the time reduction for a prescribed error, the possibility of computing simultaneously the real and imaginary parts. The behaviour of the integrands close to π/2 must be also taken 8 into account. Asymptotic developments have been made close to this value (Nontakaew et al 1997). Panel or segment source The boundary integrations are performed by an analytical method which has been shown to more accurate and efficient that the combined method mixing Gauss point method and analytical one of (Boin et al 2000). This analytical method is based on a Stokes theorem transforming the surface integral into a contour one (Bougis 1981): d2 m f (ζk +1 ) − f (ζk ) Is = ∫∫ 2 fds = ∑ C k (10). S dζ k =1 ζk +1 − ζk In this equation, use of the properties of the modified complex exponential integral functions have been done, concerning integration ∫g j (ξ ) ξ = g j (ξ ) + ln ξ d and ∫ (∫ g j ) (ξ ) ξ d ξ = g j (ξ ) + (ξ + 1) ln j ξ − ξ with ln ξ = ln ξ + i arg(ξ ) , and 0<arg(ξ)<2π, d for j=1,3 or -π<arg(ξ)<π if j=2. Concerning the derivatives, we have d 1 d2 1 1 g j (ξ ) = g j (ξ ) = g j (ξ ) − and g j (ξ ) = 2 g j (ξ ) = g j (ξ ) − + 2 . The integration of dξ ξ dξ ξ ξ the term G0 has been extensively studied in the Rankine methods. For G1, the integration on a panel with N vertices gives: 1 1 D − D D ' − D ' π 2 N 1k 1k 2k 2k ∂ / ∂x Gds = 1 A/ L dθ C K D − K D +C ' (−1)i+1(K D ' − K D ' ) (11) ∫∫ i 1 π 0 ∫ 1 + 4τ cos θ ∑ k 1 1k 2 2k k k=1 1 1k 2 2k S 2 ∂ / ∂x 2 2 2 0 A / L 2 2 K1 D1k − K2 D2k 2 2 K1 D '1k − K2 D '2k i 0 Ck for i=1,2,3 and = (q ∓ ir sin θ )(xk+1 − xk ) − (p − ir sin θ)(yk+1 − yk ) . M k (x k , y k , z k ) is C'k the node k of a panel with MN+1=M1; the outer unit normal to the panel is next = (p, q, r ) and: χk +1 (g1(ξ)d ξ )d ξ χk j ∫ ∫ Djk = k +1 k j j=1,2. ( )' are deduced from the ( ) ones by replacing χk by χj′k χj − χj j χk K j j k = z + z k + i (x − x k ) cos θ ± (y − y k ) sin θ , j=1,2 (12) χ ' j lo . Concerning the first derivatives, coefficients ( D jk ) are deduced from the previous ones 9 by replacing the two integrations by a unique one. For the second derivatives, the coefficients .. ( D jk ) involved only the functions g1. Furthermore, we have the following relations: x1=x, A=icosθ if i=1;x2=y, A=isinθ if i=2 and x3=z, A= 1 if i=3. Integrations for G2 involve terms very similar to the one of equation (11) and are given in appendix 1. For the integration on a segment, a similar expression to equation (10) can be obtained: df f (ξ2 ) − f (ξ1 ) IL = ∫ d ξ dl = ξ2 − ξ1 Cl Details can be found in (Boin et al. 2000). Numerical results Integration on an isolated panel Computations of the Green function and its first and second derivatives by the two methods developed in the previous paragraph have shown a very good agreement, (Brument & Delhommeau 1997), (Brument 1998). To test the new methods of integration, we have performed calculations of the various integrals involving the Green function and its first derivatives on an isolated panel with an unit side and with its centre of mass located at z=-.51; the unit outer normal is directed along the y-axis. Only cases where differences can be observed are presented. The calculations have been performed for a field point M located on the three different axes. Examples of calculations of the integrals for the function are plotted on Fig. 1 with the location of M on the y-axis. The two analytical integration methods give the same results ensuring the accuracy of the computations, and the results can be considered as a reference to check the results of the numerical integration methods. Far from the panel, the simple Gauss method induces very low errors with a very short computational time. Similar results are obtained for the first derivative with respect to x as shown on Fig. 2, also along the y-axis. Nevertheless, close to the peaks of the curves, numerical integration with a unique Gauss point leads to some errors and even the method with four Gauss points is unable to give accurate results close to y=0 for the integration of the imaginary part of the x–derivative. The results concerning the y-derivatives will not be presented as leading to null values. Results of the integration for M on the y-axis are presented on Fig. 3 and 4, respectively for the function and the x-derivative; they show greater difficulties when the field point is close to the free surface, the Gauss method requiring more and more points when the absolute value of zM decreases. This effect is more sensitive for the derivative than for the function. For example, four points are needed for zM<-.0.5, 16 points for zM<-.3, 36 are required for zM>-.1. Closer to 10 the free surface even 36 points are not enough to obtain a good accuracy and only analytical methods can give correct results. This last fact is very important mainly for free surface elevation calculations. The derivatives imply more difficulties to be integrated than the function, as observed in (Maury 2000). Figure 5 plots the real parts of the integrations on a waterline segment ( - 0.433≤x≤0.433; -0.25≤y≤0.25) of the function g (left graph) and of its derivative ∂ g ⁄∂x (right one) for a field point parallel to the axis x with z=-0.01 and y=2 for F=0.32 and ω = 2 . Downstream the source point, only numerical methods using more than 8 Gauss points can give accurate results. Figure 6 presents similar results for a field point describing the z axis. Even numerical integration with 16 Gauss points are unable to calculate accurately the integrals close to the free surface and only the analytical integration can be used. The agreement is the same as for a panel but integrations are more difficult to perform for a segment than for a panel. Computations in the seakeeping codes The use of the third Green's formula applied to a domain closed by the body surface, the free surface and a surface located at infinity leads to an integral equation to compute the velocity potential. Due to the radiation condition the contribution of the surface at infinity vanishes and N ∂G M ∂G Vn (x , y, z ) = ∑ σk ∫∫ ds(M ') + F 2 ∑ σl ∫ dl (13) k =1 ∆Sk ∂n M l =1 δC ∂nM l the integral on the free surface can be transformed into a line integral on the waterline. After discretisation on the hull surface in N panels with M waterline segments, for a symmetric flow (only sources σ are used), the body condition leads to the following equation after derivation: This previous integration technique has been used to compute the boundary integrals of (13). The first application presented is for a Wigley hull (with relative thickness B/L=0.1 and relative draft T/L=.0625) defined by : ( ) 4 y / B = 1 − (z /T ) 1 − (2x / L ) 1 + 0.2 (2x / L ) + (z /T ) 1 − (z /T ) 1 − (2x / L ) 2 2 2 2 8 2 in forced heave and pitch motion. Added mass and damping coefficients (Aij and Mij being respectively the damping coefficient and the added mass for force or moment i and motion j, V being the immersed volume), are plotted Fig. 7 versus the non dimensional frequency (with l0 = L , at Froude number F=0.2. Results of both codes for all coefficients are in good agreement in spite of the different number of panels, showing that the convergence has been reached : 522 for Poseidon and 750 for Aquaplus. Some differences are visible close for 11 values of τ close to the singular value of τ=ϖF=1/4. Theses results have been also compared with the experimental measurements of (Gerritsma 1988), and the computations of (Lin & Yue 1990) using 3D singularity method with a Green function satisfying the unsteady free surface boundary condition and (Wang et al 1997), using the Chapman approach with a two and half dimensional method, both in the time domain. Here also the agreement is reasonably fair; results are very good for CM33 and CM55, and the graphs show similar variation with the frequency but on nearly parallel curves for CA33 and CA55. For the cross coupling coefficients, the results are in relatively good agreement but in some cases, some discrepancy either with numerical methods or measurements can be observed. 4. Experimental set-up and test conditions An experimental device to measure forces and moments on ship model in forced oscillations of heave and/or pitch has been designed and developed at LEA-CEAT (Guyot 1995, Guyot & Guilbaud 1995), Fig. 8. This experimental set-up enables also to measure the unsteady wave pattern around this model. The system is driven by an electric motor (1) with a variable speed of rotation. The model (2) is moved through a set of counter gears and cranked belts. Pure heave is obtained through the heave rod (5) (after uncoupling the pitch rod (8) or inversely). With the 2 eccentrics and rods, a combined heave-pitch motion is possible. The model is fixed to the oscillator through a 3-component dynamometer (3), which can be replaced by a rigid modulus in order to avoid deformation of the set-up during the wave pattern measurements. The motion amplitudes are from 0 to 20 mm for heave, 0 to 6 ° in pitch, with frequencies ranging from 1 to 6 Hz. To keep a good accuracy of the added mass and damping coefficient measurements, we need to reduce the inertia forces and so an effort has been done to build models as light as possible by use of composite materials (carbon fibber); the masses are M=0.8 and 1 kg for the two L=1.2 m long series 60 models for respectively the block coefficient CB=0.6 and 0.8. The three-component dynamometer is composed of three force transducers. These transducers being sensible to forces normal to their axis, a system for uncoupling forces has been designed, so transducers are only in contact with compression or traction forces. Two ±25 daN and one ±10 daN force transducers are used. One displacement transducer records the model motion : (6) for heave and (7) for pitch. Free surface elevations are recorded using a resistive probe (4) consisting of 2 parallel chromel wires (0.2 mm diameter and 150 mm high, separated from 10 mm) held by a plexiglas 12 frame to avoid electric perturbations. The electric resistance between the two wires is function of the water elevation, the 2 wires forming a branch of a Wheatstone bridge supplied with alternative current (3.0 kHz) to prevent water electrolysis, so the output signal has to be demodulated. This probe is moved by a traversing mechanism driven by a stepper motor in the streamwise direction (x), and can be moved transversally (y) by a manual system. Measurements have been performed every 40 mm in both directions (80 mm in the y direction for y>240 mm), on only one side of the model, due to the symmetry of the flow, on a total length of 1.88 m long and 0.40 m wide. The bow of the model is located at y=0 and x=240 mm and the rear part at x=1440 mm. Close to the model, due to the size of the probe, the number of wave measurements has been reduced (for y≤80 mm for the CB=0.6 model and y≤120 mm for CB=0.8 one). The number of measurement locations is about 350 for each configuration, following Fig. 9, which plots the relative locations of the model and of the probe measurements. Force and free surface elevation measurements have been performed with a PC com- puter with a 12 bits analog to digital converter using a four channel sample-hold box; 6000 samples for global measurements and 1024 samples are recorded for each channel at data acquisition frequency of 100 Hz. Four channels are used for force measurements (motion and dynamometer), but only two (motion and wave probe) for the wave pattern. For these latter measurements, the data acquisition time is about 10s for each point. For each signal, a Fourier analysis first achieves a rough determination of the motion frequency by looking for the maximal energy spectrum. Around this first value, 100 points are distributed for ±0.1 Hz with a frequency step decreasing closer to the first determination of the frequency, and the fundamental response is obtained from the maximum amplitude of the signal with respect to the frequency. This method is more accurate than the FFT method used in (Guyot 1995), but the computational time is somewhat larger. This determination is applied to the various signals to determine the transfer function of each signal with respect to the motion, enabling to calculate the motion frequency, the signal spectra and more particularly the response at same frequency that the input. So amplitude and phase with an arbitrary origin of signals and by using a calibration table, force, moment or wave elevation amplitudes and phase lags with respect to the forced motion are computed. For added mass and damping measurements, tests are performed two times at the same frequency, once in air to record inertia forces and once in water to measure total forces, inertia plus hydrodynamic forces, the last ones being obtained by subtracting inertia forces from total forces taken into account the phase lags; then the added 13 mass and damping coefficients are calculated. For the wave height measurements, the calibration is achieved by static measurements, but takes into account the water velocity U, so 3 different calibrations (fifth order polynomial) have been determined. After a development on a small model (L=0.3 m) in the recirculating water channel of the CEAT with a 0.3 m wide test section, the experiments have been performed in the 2*1 m2 test section of the recirculating water channel of Ecole Centrale de Nantes, where the maximum water velocity is 1.7 m/s. Only pure motions have been studied with amplitudes of a3=10.8 mm (a/L=0.09) and a5=1.8 ° for U=0.13-0.7 and 1 m/s, corresponding to Froude numbers based on the model length 0.04≤F≤0.29 and frequencies of 2.5 to 6 Hz for the global measurements (Brard parameter 0.21≤τ≤3.84). For the wave measurements the lower velocity was increased till 0.4 m/s to reduce the wave reflection by the walls of the test-section. So, the frequencies were 0.87≤f(Hz)≤3.92, corresponding to Brard parameter between 0.22 and 2.5. The Reynolds numbers of the tests are between 0.2.106 and 1.2.106. Comparison of the numerical and experimental results on the series 60 models The first results ("initial calculation"), for the CB=0.6 hull at Froude number F=0.2 only for the heaving motion, are presented Fig. 10 versus the non dimensional frequency fL/U. They are obtained by neglecting the waterline integral and show large oscillations of the coefficients, probably due to the existence of irregular frequencies, as already mentioned in (Du et al 2000). By adding the waterline integral ("wi"), it is clear that the amplitude of oscillations have strongly decreased. We have then introduced a classical technique to suppress these frequencies ("tif") by adding an horizontal flat surface divided into panels inside the hull. On each of these inner panels, a condition of zero normal velocity have been written. The oscillations of the curve are nearly entirely suppressed but are still present with weak amplitudes if the waterline integral is neglected. It shows the relation between this integral and the existence of irregular frequencies. These results have been observed for both codes. The results of computations for the two CB=0.6 and 0.8 models have been compared with the measurements described in the previous paragraph. Both numerical results are for a grid of 245 panels on the half hull. Figures 11 and 12 plot the added mass and damping coefficients CMjj=Mjj/ρLn and CAjj =Ajj/ρωLn for j=3 (n=5) and 5 (with n=3) versus the non- dimensional frequency fL/U, respectively for heave and pitching forced motions at the Froude number F=0.2. On these figures, the dashed lines are for Poseidon code, the full ones for 14 Aquaplus code and the symbols correspond to the present test measurements. The two methods show very similar results with oscillations of the result for some values of the frequency. The curves are similar for both ship models (CB=0.6 or 0.8), with higher values of the coefficients as the block coefficient increases, except for CA55. Nevertheless, in some cases, some differences on the absolute values of maximum or minimum (for example, CM33 for the CB=0.8, figure 12) or on their locations can be observed (for CA55 for the CB=0.6, Fig. 11). The test results are generally in relatively good agreement with the results of the computations, except for CM55, where the test results are quite weak, may be due to problems in the measurements. Figures 13 present the different measured wave-elevation patterns for both CB=0.6 (upper graphs) and 0.8 (lower ones) models for the same test configuration (F=0.2; f=3Hz; τ=1.35); graph a is for heave and graph b for pitch. The recorded wave patterns have the same characteristics whatever are the test configurations and the model motion : two zones in V- shape with the tip in the upstream direction are visible with strong amplitude values, at the bow and at the stern (these amplitude values being stronger at the stern). The whole wave field is contained in this V-shape pattern with an opening angle and amplitudes which increase with the block coefficient. For instance, for the CB=0.8 model, just after the bow wave, a small area with low wave-amplitude appears, and three zones with strong amplitude are clearly separated along the hull. For the pitch, the back V-shape pattern presents globally higher wave amplitudes in comparison with heave, but the areas of strong wave amplitudes are reduced. It must be noticed than for heave fore and aft waves have the same phase but there is a phase lag of about 180 ° for the pitch motion. A more accurate approach of the wave-elevation contours is furnished by the Fig. 14 and 15 where respectively axial and transverse cuts of the unsteady wave pattern amplitudes for the CB=0.8 model at F=0.2 and f=4 Hz (τ=1.79), far from an irregular frequency, both for heave and pitch motions. Figure 14 plots the wave amplitude for y/L=0.2; upper graph is for heave and the lower one is for pitch. The agreement between the 2 codes is quite good, but it is necessary to have a larger number of panels (about 1000 while about 500 seems to be enough for the force calculations) than for computation of global forces. It can be observed on this figure than the amplitudes behind the stern (x/L>1.1) are greatly overestimated, effect which has not been yet explained. If values of the maxima and minima seems to be correctly predicted, their locations are not very accurate with regard to the tests. The transverse cuts, Fig. 15, are plotted for x/L=0.8 and 0.9 (close to the model stern). Upper graphs are for 15 x/L=0.9 and the lower ones, for x/L=0.8. Left ones are for heave motion and right ones for pitch. Here also, the 2 codes are in good agreement except some phase lags which may be probably due to a different number of panels. It must be noticed that the values of both components of wave elevation are weak and that the step between two measurement points are larger than some oscillations observed on the calculations. The results of computations have also been compared to present measurements. The general trend is predicted but with some discrepancies concerning the location of the maxima and minima of the curves. Conclusion We have presented a numerical and experimental study of the radiation flow about ships advancing in waves. The numerical codes have been developed with the assumption of linear theory. Both are first order panel methods using the diffraction-radiation with forward speed Green function in the frequency domain. By satisfying the body condition on the hull, we obtain a linear system of equations. Using a velocity based formulation with only a source distribution, this distribution can be calculated and then the whole flow. The main feature of these methods is an accurate computation of the coefficients of the linear system. These coefficients involve boundary integrals of the Green function and its first derivatives over panels on the hull or segments on the waterline. The boundary integrals are calculated accurately after permutation of the boundary and single Fourier integrals. The first code uses the Bessho formulation and integration of complex oscillating integrands are performed by steepest descent method. The second code uses the formulation of Guével and Bougis with the single integrals involving the complex exponential integral function calculated by a Simpson adaptative method where the local step varies with the oscillations of integrand. The accuracy of these boundary integrations have been shown, even close to the free-surface, by comparison of results of the two codes on elementary panels or waterline segments. This comparison shows that the derivatives are more difficult to integrate than the function and that for a field point close to the free-surface, it is impossible to use a numerical method of integration like the Gauss one with a reasonable accuracy. Integrations on a segment are more difficult to perform that on a panel but easier than the calculation of the function or of its derivatives. Finally these integrations have been introduced in seakeeping codes which have been applied to a Wigley hull and on two Series 60 ones with CB=0.6 and 0.8. In order to check the validity of the calculations, we have performed an experimental study in the recirculating water channel where not only global measurements as added-mass 16 and damping have been measured but also local ones. The wave pattern has been measured and amplitude and phase-lag with respect to the motion have been recorded on about 350 points on half the free surface for heave and pitch motions. The measurements have been performed for various values of the Froude number and of the motion frequency. The comparison of the numerical results with available numerical method or with the present test measurements have been done. The agreement between the two codes is excellent. For the Wigley hull, the agreement is also good with other results, particularly for the added- mass of the pure coefficients. Some discrepancies appear for the damping, but the curves versus the frequency are nearly parallel. A good prediction of the cross coupling coefficients has been obtained in spite of the low values of these coefficients and of the difficulty to measure them. Concerning the Series 60 hulls, results have shown the existence of irregular frequencies, but less abrupt than these observed at zero forward speed. Calculations performed by neglecting the waterline integral show an increase of the amplitudes of the oscillations showing that this integral has a strong damping effect on the irregular frequencies. So a technique to remove these frequencies has been used. The hulls are closed by a flat horizontal surface, slightly immersed, where a zero velocity condition has been written. These frequencies have been almost completely suppressed, at least the lower of them. Results concerning added-mass, damping and cross-coupling coefficients calculated by the two codes are in good agreement between them and with the other results available. Finally some comparison have been done concerning the unsteady wave field. Some difficulties appear for the comparison of the wave pattern amplitudes. The general shape is predicted but with some differences of maxima and minima along the hull. Downstream of the stern, the calculations overpredict the amplitude of the wave field, phenomena which has not be yet explained and will need further experiments. Appendix 1 By use of equation (10), the integration of G2 on a panel gives: 1 1 1 ∂ / ∂ x G ds = 1/ L I + I + I + I with: ∫∫ i 2 π 0 1 2 3 4 S 2 ∂ / ∂ x 2 1/ L2 i 0 17 1 E − E E ' − E ' N 3k 3k θc' 4k 4k −i C I 1 = A ∫ ∑ Z E − Z E + C ' (− 1)i +1 (Z E ' − Z E ' ) d θ 3 3k 4 4k k 4k 0 2 4 τ cos θ − 1 k =1 k 3 3k 4 A 2 2 2 Z 3 E 3k − Z 4 E 4k 2 Z 3 E ' 3k − Z 4 E '4 k 1 F − F F ' − F ' θc −αc N 1k 1k 2k 2k −i I 2 = A ∫ C ∑ Z F − Z F + C ' (−1) (Z F ' − Z F ' ) d θ i +1 ' 4 2k 2k 4τ cos θ − 1 k =1 k 2 3 1k k 2 3 1k 4 2 θc A 2 Z 3 F1k − Z 4 F2k 2 Z 3 F '1k − Z 4 F '2k 1 1 χc +1 χ 'c +1 χ 'c −2π (1 − i ) αc k k k k N e −e χc e −e I 3 = Ac Kc τ sin θc ∑ C k χk +1 − χk + (−1)i +1 C 'k χ 'k +1− χk c 2 2 k =1 c c 1 c Ac Kc 1 H − H H ' − H ' π2 N 1k 1k 2k 2k 1 I 4 = A ∫ ∑ C k K H − K H + C ' (−1) (K H ' − K H ' ) d θ , i +1 1 − 4τ cos θ k =1 3 1k 4 2k k 3 1k 4 2k 2 θc +αc 2 2 A 2 2 K 3 H 1k − K 4 H 2k K 3 H '1k − K 4 H '2k Coefficients Ejk, Fjk, F2k, H1k and H2k are given by: χk +1 χk +1 χk +1 (g2(ξ)d ξ )d ξ χk (g1(ξ)d ξ )d ξ χk (g 3 (ξ)d ξ )d ξ χk j ∫ ∫ ∫ ∫ 3 ∫ ∫ 4 E jk = k +1 j j=3,4, F1k = k +1 k 3 , F2k = k +1 4 ; χj − χj k χ3 − χ3 χ4 − χ4 k χ5 +1 k χ6 +1 k ∫ ∫ (g 3 (ξ)d ξ )d ξ χ k ∫ ∫ (g1(ξ)d ξ )d ξ χ k H 1k = 5 , H 2k = 6 . χ5 +1 − χ5 k k χ6 +1 − χ6 k k In the expressions (12) for χk and χj′k , Kj must be replaced by Zj if j=3,4 and by Kj-2 if j=5,6; j χck K furthermore k = c z + z k + i (x − x k ) cos θc ± (y − y k ) sin θc and Ac=icosθc if i=1; χ 'c Lo Ac=isinθc if i=2 and Ac=1 if i=3. Nomenclature Aij : damping 18 A,Ac : coefficients used in the integration over panels B : ship width Ck = (q ∓ ir sin θ ) xk+1 − xk − (p − ir sin θ)(yk+1 − yk ) : coefficient C'k ( ) CAij, CMij : damping and added-mass coefficients CB : block coefficient χk +1 (g1(ξ)d ξ )d ξ χk j ∫ ∫ D jk = k +1 k j : coefficients j=1,2 χj − χj Ei : complex-valued exponential integral function χk +1 (g2 (ξ)d ξ )d ξ χk j ∫ ∫ E jk = k +1 k j : coefficients j=3,4 χj − χj χ 3 +1 k ( g m (ξ )dξ ) dξ ∫ ∫ χ3k Flk = : coefficients with m=1 if j=1 and m=3 if j=2 χ 3k +1 − χ 3k k1e k1ξ f1(θ) = 1 + 4τ cos θ 1 k1ξl e −i 2Sl k1 i β (Y1 −Y2 ) e k1ξl f1,l (θ ) = or f (θ ) = 2πK 0 ψl (θ ) 1 + 4τ cos θ 2π (ξl +1 − ξl ) 1 + 4τ cos θ 1,l F =U / g 0 : Froude number g : gravity acceleration g : Froude number dependant part of the Green function gi (i=1,3) : modified complex-valued exponential integral function G=G0+G1+G2 : Green function θ2 (P ,P ') I = ∫ θ1 − sgn c.f1 (θ ).d θ : integral in the steepest descent method, θ-space 19 IS = ∫∫ gds , I L = ∫ gdl : integrals of the Green function or derivatives on a panel or a ∆S δC segment χk +1 (g3 (ξ)d ξ )d ξ χk 5 ∫ ∫ H 1k = k +1 k 5 : coefficients χ5 − χ5 k1,k2 : poles of the Green function in the steepest descent method Kc = 4ω 2 K0= g /U 2 : wave number Kj : poles of the Green function in the formulation of Guevel and Bougis, j=1,5 (k1=K1 and k2=K2) l0 : reference length ln : length of a single interval L : ship length M=k1cosθ M = sgn c.k1.cos θ = sgn c.m :space to calculate the steepest descent method Mij :added-mass N : number of vortices for a panel P(x,y,z),P'(x',y',z') : field and source points Ql, l=1,N : vortices of a panel R ( x − x ' ) + ( y − y ') + ( z ∓ z ') 2 2 2 = : distances to the source point and its symmetrical with R1 respect to the free surface sgnc=sign[Re(cosθ)] sgns=sign[Re(sinθ)] 20 sign(t)=-1 if t<0; 0 if t=0;1 if t>0 Sl : area of triangle (Ql,Ql+1,Ql-1) T : ship draft U : forward speed of the ship V : immersed volume of the hull Vn : normal velocity x,y,z : reference frame fixed to the body X = K 0 (x − x ');Y = K 0 y − y ' ; Z = K 0 (z + z ') ( l ) Xl = K 0 x P − xQ ;Yl = K 0 yP − yQ , Zl = K 0 z P + zQ l ( l ) X η = X1 + η (Xl − X1 ); Yη = Y1 + η (Yl −Y1 ); Z η = Z 1 + η (Zl − Z1 ) Xl ' = Xl − X1; Yl ' = Yl −Y1; Zl ' = Zl − Z1 Xel = XQ − XQ ; Yel = YQ −YQ ; Zel = ZQ − ZQ (l +1) l (l +1) l (l +1) l 1 − 2τ cos θ + (−1)j 4τ cos θ − 1 Zj = for j = 3, 4 : complementary poles 2F 2 cos2 θ αc=10-6 : coefficient used to calculate the Green function β,β' : starting points for the integration in the M-space Z ε (P, P ') = arg sh X 2 +Y 2 εtot,,ε : absolute error on total interval L or elementary interval with length ln ζ,,ζ' θ : complex variable θc,θ'c : limits of integration for the Fourier integrals 21 π + arccos 1 1 − if τ > θ1 = 4τ 4 1 1 −π − i arg ch if τ < 4τ 4 π θ2 (P , P ') = ϕ (P , P ') − − i ε (P, P ') 2 Ye ± Xe 2 + Ye 2 + Ze 2 l θl = 2 arctan l l l Xel − i.Zel X ηYl '− Xl 'Yη Zl ' (X η + Yη2 ) − Z η (X η Xl '+ Yl 'Yη ) 2 νl (η ) = +I 2 X η + Yη2 (X η2 + Yη2 ) (X η2 + Yη2 + Z η2 ) ξ = Z + i (X cos θ + Y sin θ ) or ξ (ξ ') = z + z '+ i [(x − x ') sin θ ∓ (y − y ') cos θ ] / l 0 ξl = ξ (P,Ql ) σ : source density τ=ωeU/g : Brard parameter X ϕ (P, P ') = arccos X 2 +Y 2 2 m φ (M ) = (m − τ Z + i.m.X + i.sgnsY (m − τ )2 . )2 1− ( m − τ )2 χck K = c z + z k + i (x − x k ) cos θ ± (y − y k ) sin θ k c χ 'c Lo c ψl (θ ) = (ξl +1 − ξl )(ξl −1 − ξl ) ωe : encounter frequency ω = ωe 0 /g 22 References Ba, M. and Guilbaud, M. 1995 A fast method of evaluation for the translating and pulsating Green's function. SHIP TECHNOLOGY RESEARCH, 42, 2, 68-80. Ba M., Boin J.P., Delhommeau G., Guilbaud M. et Maury C. 2001a Calcul de tenue à la mer avec la fonction de Green de diffraction-radiation avec vitesse d'avance. Proceedings: 8èmes Journées de l'Hydrodynamique, 125-140. Ba M., Boin J.P., Delhommeau G., Guilbaud M. et Maury C. 2001b Sur l'intégrale de ligne et les fréquences irrégulières dans les calculs de tenue à la mer avec la fonction de Green de diffraction-radiation avec vitesse d'avance. C.R. ACAD. SCI. PARIS, 329, IIb, 141-181. Bougis, J. 1981 Bessho, M. 1977 On the Fundamental Singularity in the Theory of Ship Motions in a Seaway. MEMOIRS OF THE DEFENCE ACADEMY JAPAN, XVII, 3, 95- 105. Boin, J.P., Guilbaud, M. and Ba, M. 2000 Sea-keeping computations using the ship motion Green's function. Proceedings: ISOPE2000 Conference, Seattle, 398-405. Bougis J. 1981 Etude de la diffraction-radiation dans le cas d'un flotteur indéformable animé d'une vitesse moyenne constante et sollicité par une houle sinusoïdale de faible amplitude. Thèse de doctorat, Université de Nantes. Brument, A. et Delhommeau, G. 1997 Evaluation numérique de la fonction de Green de la tenue à la mer avec vitesse d'avance. Proceedings: 6èmes Journées de l'Hydrodynamique, Nantes, 147-160. Brument, A. 1998 Evaluation numérique de la fonction de Green de la tenue à la mer. Thèse de Doctorat, École Centrale de Nantes. Chen, X.B. and Noblesse, F. 1998 Super Green function. Proceedings: 22nd Symposium on Naval Hydrodynamics, Washington, 860-874. 23 Du, S.X., Hudson, D.A., Price, W.G. and Temarel, P. 1999 Comparison of numerical evaluation techniques for the hydrodynamic analysis of a ship travelling in waves. Trans. RINA, 141, 236-258. Du, S.X., Hudson, D.A., Price, W.G. and Temarel, P. 2000 A validation study on mathematical models of speed and frequency dependence in seakeeping. Trans. RINA, 214, 181-202. Delhommeau, G., Ferrant, P. and Guilbaud, M. 1992, Calculation and measurement of forces on a high speed vehicle in forced pitch and heave. APPLIED OCEAN RESEARCH, 14, 119- 126. Gerritsma, J. 1988 Motions, wave loads and added resistance in waves of two Wigley hull forms. Delft University of Technology, Ship Hydromechanics laboratory, Report n° 804. Guével, P. and Bougis, J. 1982 Ship-Motions with Forward Speed in Infinite Depth. INT. SHIPBUILDING PROG., 29, 103-117. Guttman, C. 1983 Etude théorique et numérique du problème de Neumann-Kelvin pour un corps totalement immergé. Rapport de Recherche n°177, ENSTA. Guilbaud, M. Boin, J.P., and Ba, M. 2000 Frequency domain numerical and experimental investigation of forward speed radiation by ships. Proceedings: 23rd Symposium on Naval Hydrodynamics, Val de Reuil, Tuesday session, 110-125. Guyot, F. 1995 Étude expérimentale de la résistance ajoutée d'une maquette de navire soumise à des oscillations forcées harmoniques: étude du champ de vagues instationnaires associé. Thèse de Doctorat, Université de Poitiers. Guyot, F. and Guilbaud, M. 1995 Force and free surface elevation measurements on a series th 60 CB=0.6 ship model in forced oscillations. Proceedings: 5 ISOPE Conference, The Hague, 507-514. 24 Iwashita, H. and Okhusu, M. 1989 Hydrodynamic Forces on a Ship Moving at Forward Speed in Waves. J.S.N.A. Japan, 166, 87-109. Iwashita, H. 1992 Evaluation of the Added-Wave-Resistance Green Function Distributing on a Panel. Mem. Fac. Eng. Hiroshima Univ., 11, 2, 21-39. Ishawita, H., Ito, A., Okada, T., Okhusu, M., Takaki, M. and Mizoguchi, S. 1993 Waves forces acting on a blunt ship with forward speed in oblique sea (2). T. Soc. Naval Arch. Japan, 173, 195-208. Lin, W.-M., and Yue, D. 1990 Numerical solutions for Large-Amplitude ship motion in the time domain. Proceedings: 18th Symp. on Naval Hydrodynamics, Ann Arbour, 41-66. Maury, C. 2001, Etude du problème de tenue à la mer avec vitesse d'avance quelconque par une méthode de singularités de Kelvin. Thèse de Doctorat, École Centrale de Nantes. Nontakaew, U., Ba, M. and Guilbaud, M. 1997 Solving a radiation problem with forward speed using a lifting surface method with a Green's function. AEROSPACE SCIENCE AND TECHNOLOGY, 8, 533-543. Okhusu, M. and Wen, G. 1996 Radiation and diffraction waves of a ship at forward speed. Proceedings 21th Symposium on Naval Hydrodynamics, Trondheim, 29-44. Okhusu, M. 1998 Validation of theoretical methods for ship motions by means of experiment. Proceedings: 22nd Symposium on Naval Hydrodynamics, Washington, 341-358. Wang, C.-T., Horng, S-.J. and Chiu, F.-C. 1997 Hydrodynamic forces on the advancing slender body with speed effects. INT. SHIPBUILDING PROG., 44, 438, 105-126. List of figures Figure 1 Integrals of the Froude dependant part of the Green function for a field point describing the y-axis 25 Figure 2 Integrals of the first x-derivative of Froude dependant part of the Green function for a field point describing the y-axis Figure 3 Integrals of the Froude dependant part of the Green function for a field point describing the z-axis Figure 4 Integrals of the x-first derivative of the Froude dependant part of the Green function for a field point describing the z-axis Figure 5 Integrals of the Froude dependent part of the Green function and x-derivative for a field point describing a parallel to the x-axis (z=-0.01;y=2) Figure 6 Integrals of the Froude dependent part of the Green function and x-derivative for a field point describing the z-axis Figure 7 Added mass and damping coefficients on the Wigley hull (F=0.2) Figure 8 Experimental set-up and motion oscillator Figure 9 Locations of the wave measurements Figure 10 Influence of the waterline integral on the added mass and damping coefficients on the Series 60 CB=0.8 (F=0.2) Figure 11 Added mass and damping coefficients on the Series 60 CB=0.6 (F=0.2) 26 Figure 12 Added mass and damping coefficients on the Series 60 CB=0.8 (F=0.2) a) heave motion (CB=0.6: top; CB=0.8: low) b) pitch motion (CB=0.6: top; CB=0.8: low) Figure 13 Measured waves amplitudes (F=0.3, f=3Hz, τ=1.92) Figure 14 Comparison of longitudinal cut of wave elevation y/L=0.2 (Series 60 CB=0.8 hull; F=0.2; τ=1.79) Figure 15 Comparison of transverse cut of wave elevation (Series 60 CB=0.8 hull; F=0.2; τ=1.79) 27 z -0.5 0.5 x -0.01 -1.01 ℜ ∫∫ gds ℑ ∫∫ gds S S Aquaplus Aquaplus Poseidon 0.06 Poseidon 0.15 Numerical integration 1 Gauss point Numerical integration 1 Gauss point Numerical integration 4 Gauss points 0.05 Numerical integration 4 Gauss points Numerical integration 16 Gauss points 0.04 Numerical integration 16 Gauss points 0.03 0.02 0.1 0.01 0 -0.01 0.05 -0.02 -0.03 -0.04 -0.05 0 -0.06 -0.07 -0.08 -0.05 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 y y Figure 1 Integrals of the Froude dependant part of the Green function for a field point describing the y-axis ∂g ∂g ℜ ∫∫ ds ℑ ∫∫ ds S ∂x S ∂x Aquaplus Aquaplus Poseidon 0.14 Poseidon 0.12 Numerical integration 1 Gauss point Numerical integration 1 Gauss point Numerical integration 4 Gauss points 0.12 Numerical integration 4 Gauss points 0.1 Numerical integration 16 Gauss points Numerical integration 16 Gauss points 0.08 0.1 0.06 0.08 0.04 0.06 0.02 0.04 0 0.02 -0.02 -0.04 0 -0.06 -0.02 -0.08 -0.04 -0.1 -0.06 -0.12 -0.08 -0.14 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 -5 0 5 y y 28 Figure 2 Integrals of the first x-derivative of Froude dependant part of the Green function for a field point describing the y-axis ℜ ∫∫ gds ℑ ∫∫ gds S S 0.075 0.4 Aquaplus Aquaplus P oseidon P oseidon N um erical integration 1 G auss point 0.05 N um erical inte gration 1 G auss point N um erical integration 4 G auss points N um erical inte gration 4 G auss points N um erical integration 1 6 G auss points 0.3 N um erical integration 3 6 G auss points 0.025 N um erical inte gration 1 6 G auss points N um erical inte gration 3 6 G auss points 0 0.2 -0.025 -0.05 0.1 -0.075 -0.1 0 -1 -0.75 -0.5 -0.25 0 -1 -0.75 -0.5 -0.25 0 z z Figure 3 Integrals of the Froude dependant part of the Green function for a field point describing the z-axis ∂g ∂g ℜ ∫∫ ds ℑ ∫∫ ds S ∂x S ∂x 0.3 0.3 Aquaplus P oseidon N um erical inte gration 1 G auss point 0.2 0.2 N um erical inte gration 4 G auss points N um erical inte gration 1 6 G auss points 0.1 N um erical inte gration 3 6 G auss points 0.1 0 0 -0.1 -0.1 -0.2 Aquaplus -0.2 P oseidon N um erical integration 1 G auss point -0.3 N um erical integration 4 G auss points -0.3 N um erical integration 1 6 G auss points -0.4 N um erical integration 3 6 G auss points -0.4 -0.5 -0.5 -0.6 -1 -0.75 -0.5 -0.25 0 -1 -0.75 -0.5 -0.25 0 z z Figure 4 Integrals of the x-first derivative of the Froude dependant part of the Green function for a field point describing the z-axis 29 ∂g ℜ ∫ gdc ℜ∫ dc δCl δCl ∂x 0.5 5 0.4 4 0.3 3 0.2 2 1 0.1 0 0 -1 -0.1 Aquaplus Aquaplus Poseidon -2 Poseidon -0.2 Numerical method 1 Gauss point Numerical method 1 Gauss point Numerical method 4 Gauss points -3 Numerical method 4 Gauss points -0.3 Numerical method 8 Gauss points Numerical method 8 Gauss points Numerical method 16 Gauss points -4 Numerical method 16 G auss points -0.4 -5 -8 -6 -4 -2 0 2 -8 -6 -4 -2 0 2 x x Figure 5 Integrals of the Froude dependent part of the Green function and x-derivative for a field point describing a parallel to the x-axis (z=-0.01;y=2) ∂g ℜ ∫ gdc ℜ∫ dc δCl δCl ∂x 0.3 7 Aquaplus Poseidon Aquaplus Numerical integration 1 Gauss point 6 Poseidon Numerical integration 4 Gauss points Numerical integration 1 Gauss point Numerical integration 8 Gauss points Numerical integration 4 Gauss points 5 Numerical integration 16 Gauss points Numerical integration 8 Gauss points 0.2 Numerical integration 16 Gauss points 4 3 2 0.1 1 0 0 -1 -0.2 -0.1 0 -0.2 -0.1 0 z z Figure 6 Integrals of the Froude dependent part of the Green function and x-derivative for a field point describing the z-axis 30 2.25 2.25 Poseidon 522 panels Lin Yue (1990) 2 2 Tests Gerritsma ex (1988) Chun-Tsung (1997) 1.75 Aquaplus 750 panels 1.75 1.5 1.5 A33(L/g)0.5/ρV M 33/ρV 1.25 1.25 1 1 0.75 0.75 0.5 0.5 0.25 0.25 0 0 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 ϖ ϖ 0.16 0.09 0.14 0.08 0.07 0.12 0.06 A55(L/g) 0.5 /ρVL2 0.1 2 M 55 /ρVL 0.05 0.08 0.04 0.06 0.03 0.04 0.02 0.02 0.01 0 0 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 ϖ ϖ 0.2 0.05 0 0.15 -0.05 A53(L/g)0.5/ρVL 0.1 M 53 /ρVL -0.1 -0.15 0.05 -0.2 0 -0.25 -0.05 -0.3 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 ϖ ϖ 0.05 0.6 0 -0.05 0.4 A35/(L/g) 0.5 /ρVL M 35/ρVL -0.1 -0.15 0.2 -0.2 -0.25 -0.3 0 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 ϖ ϖ Figure 7 Added mass and damping coefficients on the Wigley hull (F=0.2) 31 1 Motor – 2 Ship model 3 Dynamometer or rigid modulus 4 Wave probe 5 Heave rod 6 Heave transducer 7 Pitch transducer 8 Pitch rod Figure 8 Experimental set-up and motion oscillator Figure 9 Locations of the wave measurements 32 0.012 initial calculation (without wi and tfi) 0.008 calculation with wi without tif initial calculation (without wi and tif) 0.011 calculation with wi and tif 0.007 calculation with wi without tif calculation without wi with tif calculation with wi and tif 0.01 0.006 calculation without wi wit tif 0.009 0.005 0.008 0.004 0.007 0.003 CM 33 CA33 0.006 0.002 0.005 0.001 0.004 0 0.003 0.002 -0.001 0.001 -0.002 0 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 fL/U fL/U Figure 10 Influence of the waterline integral on the added mass and damping coefficients on the Series 60 CB=0.8 (F=0.2) 0.008 0.012 0.007 Aquaplus Tests 0.01 0.006 Poseidon 0.005 0.008 0.004 0.003 CM 33 CA33 0.006 0.002 0.001 0.004 0 -0.001 0.002 -0.002 0 -0.003 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 fL/U fL/U 0.0007 0.0006 0.0006 0.0005 0.0005 0.0004 0.0004 0.0003 CM 55 CA55 0.0003 0.0002 0.0002 0.0001 0.0001 0 0 -0.0001 -0.0001 -0.0002 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 fL/U fL/U Figure 11 Added mass and damping coefficients on the Series 60 CB=0.6 (F=0.2) 33 0.025 0.025 Aquaplus Tests 0.02 0.02 Poseidon 0.015 0.015 CM33 CA33 0.01 0.01 0.005 0.005 0 0 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 fL/U fL/U 0.0015 0.001 0.001 0.0005 CM 55 CA55 0.0005 0 0 -0.0005 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 fL/U fL/U Figure 12 Added mass and damping coefficients on the Series 60 CB=0.8 (F=0.2) C = 6 , 0 = C 0 6 , C = , 0 6 B C = 0 , 6 B C B = 0 .6 C 0 = 6 , C B 0 = , 6 h /L 0 .0 0 7 0 0 .0 0 6 3 0 .0 0 5 6 0 .0 0 4 9 0 .0 0 4 2 0 .0 0 3 5 0 .0 0 2 8 0 .0 0 2 1 0 .0 0 1 4 0 .0 0 0 7 C B= 0 .8 a) heave motion (CB=0.6: top; CB=0.8: low) 34 C = , 6 0 6 C =0. B C =0. 6 B C = 0 6 , C = 0 , 6 B C = 0 , 6 B C B = 0 .6 C 0 = 6 , C B = 0 , 6 C =0. B 6 h /L L h / 0 .0 0 7 0 0 .0 0 6 3 0 .0 0 5 6 0 .0 0 4 9 0 .0 0 4 2 0 .0 0 3 5 0 .0 0 2 8 0 .0 0 2 1 0 .0 0 1 4 0 .0 0 0 7 C B = 0 .8 b) pitch motion (CB=0.6: top; CB=0.8: low) Figure 13 Measured waves amplitudes (F=0.3, f=3Hz, τ=1.92) 0.015 0.0125 C B = 0 ; a 3 /L= 0 0 9 ; F = 0 .2 ; f= 4 H z ; τ= 1 .7 9 y/L= 0 .2 0.01 0.0075 A/L 0.005 0.0025 0 -0.0025 0 0.5 1 1.5 x/L 0.015 C B = 0 .8 ; a 5 = 1 .8 °; F = 0 .2 ; f= 4 H z ; τ= 1 .7 9 y/L= 0 .2 0.01 Aquaplus P oseidon T e sts A/L 0.005 0 -0.005 0 0.5 1 1.5 x/L Figure 14 Comparison of longitudinal cut of wave elevation y/L=0.2 35 (Series 60 CB=0.8 hull; F=0.2; τ=1.79) X/L=0.9 Heave X/L=0.9 Pitch 0.014 0.018 0.016 0.012 0.014 0.01 0.012 0.008 0.01 A/L A/L Aquaplus 0.006 0.008 Poseidon Tests 0.006 0.004 0.004 0.002 0.002 0 0 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 y/L y/L X/L=0.8 Heave X/L=0.8 Pitch 0.012 0.016 0.014 0.01 0.012 0.008 0.01 A/L A/L 0.006 0.008 0.006 0.004 0.004 0.002 0.002 0 0 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 y/L y/L Figure 15 Comparison of transverse cut of wave elevation (Series 60 CB=0.8 hull; F=0.2; τ=1.79) 36