Prediction of the Hydrodynamic Forces for a Ship Oscillating in Calm Water by an Improved Higher Order Rankine Panel Method

: This paper presents a frequency-domain Rankine source method based on a biquadratic B-spline scheme with an improved radiation mechanism. The improved radiation mechanism, based on the introduction of spatially varying Rayleigh artiﬁcial damping in addition to the simpliﬁed Seto’s radiation boundary conditions, is considered for modeling radiation of generated waves at various τ conditions, where τ = ω U / g including the undercritical condition ( τ < 0.25); this condition is present when a ship undergoes slow translation or low oscillatory frequency. In evaluations, the proposed method yields accurate solutions for unsteady ﬂows produced by an oscillating, translationally moving submerged singularity. The radiation problem induced by a RIOS bulker is solved to have the resultant added mass and damping coefﬁcients for further comparisons with the experimental data and the public numerical prediction by a simpliﬁed combined method at a wide τ region.


Introduction
Green function and Rankine source methods are commonly used to solve flow problems in marine engineering. The advantages of the Green function method are twofold: first, the boundary integral equation is constructed from only the wetted body surface, and second, the free surface boundary condition and the radiation boundary condition at infinity are automatically satisfied. The research studies in [1][2][3] have successfully applied the Green function method to seakeeping problems. However, a converged and efficient numerical way to evaluate the Green function singularity and the proper numerical treatment of the waterline integral are critical to Green function method; meanwhile they are thought of as limitations in usage of Green function implementation. To address these limitations, the Rankine source method, where the Rankine source is used as the singularity, is considered. The Rankine source method facilitates effective evaluation of the employed singularity and there is no need for treatment of the waterline integral; this method is thus thought to be promising in practical applications.
The pioneer work for use of the Rankine source method is from Dawson [4] for predicting the resistance of a body sailing in calm water. Dawson's success in solving this steady flow problem began the prevalent use of the Rankine source method for solving unsteady flow problems. Bai and Yeung [5] are the pioneers credited for application of the frequency-domain Rankine source method in solving unsteady flow problems, and achieving good prediction of two axi-symmetric three-dimensional bodies, such as sphere and ellipsoid, with the existing results. As illustrated by Bai and Yeung [5], in addition to the infinite-depth case, the Rankine source method can be applied in the finite-depth case by constructing the relevant boundary equation with the additional Rankine singularity distribution over the bottom and taking the finite-depth dispersion relation into account for determining the corresponding wave number and wave length at constant finite depth.
Since then, frequency-domain analyses of the realistic ship's unsteady flow problems using the Rankine source method have been reported by Nakos [6], Bertram [7], Yasukawa [8], Iwashita and Ito [9], Yang et al. [10], Söding [11,12], Lyu and Moctar [13], Kim and Kim [14] and others. Furthermore, the time-domain Rankine source method using the linear, weakly nonlinear or weak scatter assumption on free surface boundary conditions are considered by Kring [15], Kim and Kim [16], Huang [17], Singh and Sen [18] and Datta et al. [19], for attempts to deal with the flow problem induced by the ship in large amplitude motion. To handle the full nonlinearity existing in interaction between the body and water, the method of mixed Eulerian-Lagrangian, in which the evolution in the dynamic and kinematic free surface boundary conditions are considered, is introduced by Longuet-Higgins and Cokelet [20] to simulate the two-dimensional wave profile variation in time. Its succeeding application in the three-dimensional problems are investigated by Cao [21], Beck [22,23], Kashiwagi [24] and Abbasnia [25]. However, even though these variants of the time-domain Rankine source methods are proposed and have progressed in past decades, the frequency-domain method is considered to be the most mature and efficient way of dealing with various flow conditions in engineering practices. For this reason, the present paper aims to use the frequency-domain Rankine source method for wider applications in the marine field.
To satisfy the free surface boundary condition, the computation domain required for constructing the boundary integral equation covers the body surface and the extent of the free surface. Although the allocation of the free surface implies more computational cost, it provides good flexibility in satisfying the variant of the free surface boundary condition. However, the accompanying problems of wave distortion due to free surface discretization and avoidance of wave reflection must be addressed. Jensen et al. [26] and Sclavounos [27] have improved the method by reducing numerical dispersion and dissipation. To ensure that the radiation boundary condition is satisfied on the free surface, Dawson [4] employed the upwind differential, Raven [28] analyzed the source shift, Sclavounos [27] applied the rigid-lid condition and Yasukawa [8] suggested introduction of the Rayleigh artificial damping term. However, these numerical studies focused on solving the flow problem under overcritical conditions, in which τ(= ωU/g) > 0.25 and all the generated wave systems propagate backward; their proposals may not be appropriate for the flow problem under undercritical conditions, where τ < 0.25 and part of the generated wave system scatters ahead of the translating vessel. Herein, τ is the oscillatory frequency of the vessel, U is the forward speed, and τ is thought of as an index for interpreting the generated wave system. Das and Cheung [29] and Yuan et al. [30] have attempted to address the aforementioned limitation by formulating an extended Sommerfeld radiation boundary condition derived from the forward-speed-induced Doppler shift in circular waves at an additional control surface. In a recent study by Iwashita et al. [31,32], the joint conditions for describing wave radiation are proposed to be blended with the Rankine source method. These joint conditions, derived based on the concept of a combined domain, were employed at the control surfaces and defined as the induced velocity potential and its normal derivative simulated by the allocation of Green function singularity. Furthermore, the joint condition is suggested as the resultant ratios of the normal derivatives to velocity potential, and such a technique, which is called the simplified combined method, is used for seakeeping prediction of ship, e.g., the studies in [32,33]. However, implementation of the extended Sommerfeld radiation or joint conditions requires good control surface allocation and increases the number of unknowns and computational cost.
In order to have the continuous description of flow and pressure fields, the velocity potential modeled in spline form is considered in constructing the boundary integral equation; this lies at the heart of the so-called higher-order method. Hsin et al. [34] constructed a model using B-splines for a two-dimensional flow problem. Scholars such as Nakos [6], Maniar [35], Coaxley [36], and Kim and Kim [37] have extended this insight to threedimensional steady and unsteady flow problems. Furthermore, Nakos and Sclavounos [38] discuss the stability analysis, which is based on the dispersion relation in continuous and discretized problems, for the B-spline scheme and propose the stability criterion.
The present study presents a frequency-domain Rankine source method that is based on a B-spline scheme. This method can be used to solve unsteady flow problems regardless of the τ condition (making this method practical for use). This method contrasts with the unsteady Rankine source approaches in the literature that are applicable only to the overcritical condition. Specifically, the present study constructs the resultant boundary integral equation by employing simplified Seto's explicit conditions and Rayleigh artificial damping to accommodate wave radiation. In evaluation of its effectiveness, the proposed method was successfully used to solve for the flow field generated by a single submerged disturbance source under undercritical and overcritical conditions [39]. The present biquadratic B-spline scheme is taken into account for predicting the hydrodynamic forces and moments induced by forced motion of a realistic bulk carrier; its predicted coupled added mass and damping coefficients are compared with the experimental data and the numerical result of the simplified combined method to investigate applicability of the present method in a broader region of τ. Figure 1 illustrates the boundary problem of an oscillating vessel translating at steady forward speed U in deep calm water. The oscillatory frequency is represented as ω. The adopted right-handed Cartesian coordinate frame o − xyz undergoes translational motion with the vessel on the positive x-axis, pointing to bow; its positive z-axis points upward with the frame origin o located on the gravity center of the vessel at still water. As presented in Figure 1, the finite computation extents conducted for the Rankine source panel scheme includes the free surface and wetted surface of the vessel, and herein are denoted by Σ F and Σ H , respectively. constructed a model using B-splines for a two-dimensional flow problem. Scholars such as Nakos [6], Maniar [35], Coaxley [36], and Kim and Kim [37] have extended this insight to three-dimensional steady and unsteady flow problems. Furthermore, Nakos and Sclavounos [38] discuss the stability analysis, which is based on the dispersion relation in continuous and discretized problems, for the B-spline scheme and propose the stability criterion.

Boundary Value Problem
The present study presents a frequency-domain Rankine source method that is based on a B-spline scheme. This method can be used to solve unsteady flow problems regardless of the τ condition (making this method practical for use). This method contrasts with the unsteady Rankine source approaches in the literature that are applicable only to the overcritical condition. Specifically, the present study constructs the resultant boundary integral equation by employing simplified Seto's explicit conditions and Rayleigh artificial damping to accommodate wave radiation. In evaluation of its effectiveness, the proposed method was successfully used to solve for the flow field generated by a single submerged disturbance source under undercritical and overcritical conditions [39]. The present biquadratic B-spline scheme is taken into account for predicting the hydrodynamic forces and moments induced by forced motion of a realistic bulk carrier; its predicted coupled added mass and damping coefficients are compared with the experimental data and the numerical result of the simplified combined method to investigate applicability of the present method in a broader region of . Figure 1 illustrates the boundary problem of an oscillating vessel translating at steady forward speed U in deep calm water. The oscillatory frequency is represented as . The adopted right-handed Cartesian coordinate frame − undergoes translational motion with the vessel on the positive x-axis, pointing to bow; its positive z-axis points upward with the frame origin o located on the gravity center of the vessel at still water. As presented in Figure 1, the finite computation extents conducted for the Rankine source panel scheme includes the free surface and wetted surface of the vessel, and herein are denoted by Σ and Σ , respectively.  The unsteady motion response of the vessel in calm water and the associated ambient flow are considered to be periodical with the ship oscillating frequency given by ω. The forced motion responses of the vessel are described in the six degrees of freedom oscillatory amplitudes of the vessel in calm water and are denoted by ξ j (j = 1~6), which stands for surge, sway, heave, roll, pitch and yaw movement, respectively.

Boundary Value Problem
The fluid throughout the domain in deep water, in which the Laplace equation is satisfied, is inviscid, incompressible flow and irrotational. Thus, the flow field can be described in the framework of velocity potential in the present study. Let t denote time and → p = (x, y, z) means the position vector. The total induced flow potential Ψ → x , t by the translating vessel in calm water is expressed as where i = √ −1; Φ is the double-body flow and is considered as the basis flow resulting from the presence of the translating vessel in rigid water plane. The velocity potential φ and ϕ stand for the steady and unsteady flow. Furthermore, the radiation flow component of unsteady potential flow reads as follows in which ϕ j (j = 1~6) stands for the velocity potential per unit j-mode motion displacement.
In direct methodology of the Rankine source scheme, the unsteady velocity potential in fluid domain is expressed in terms of a normal dipole distribution of moment ϕ and a source distribution of strength ϕ/∂n distributed over the boundary surfaces of Σ F and Σ H , with the Rankine source as the kernel singularity, as in the following boundary integral equation: where ϕ stands for either one of ϕ j , (j = 1~6), → p = (x, y, z) and → q = (x , y , z ) represents the field point in the fluid and the source point over the boundary surfaces, respectively. that is a product of two univariate quadratic polynomials described in more detail in the later section. Therefore, the spline coefficients σ m are the unknowns in the present solving of the boundary value problem.
A series of boundary integral equations raised at various field points → p collocated on, and end conditions around, the boundary surfaces Σ F and Σ H are collected to construct a linear system of simultaneous equations for solving the spline coefficients for approximating velocity potential over the boundary surfaces. Once the solution of spline coefficients is obtained, we can have the flow potential, the spatial derivatives and the hydrodynamic forces acting on the vessel accordingly.

Boundary Conditions
To completely construct the unsteady wave flow problem, the flow disturbance due to the steady translation of the body is considered in the derivations of relevant boundary conditions and should be solved prior to solving the unsteady flow problem.
As mentioned in the preceding section, the velocity potentials Φ and φ denote the double-body flow and the steady flow, which are induced due to the presence of the body on free surface and the steady translation. Normally, the steady potential is disregarded in solving an unsteady flow problem because the double-body flow is considered to be the basis flow and is the dominant term. Hence, the relevant velocity vector due to the body steady translation is written as by ∇Φ, normalized by forward speed U. The phenomenon of no flow flux penetrating through the boundary surface gives the boundary conditions for the double-body flow as: where n 1 denotes the x-component of normal vector on the hull. On the exact free surface, the fluid flow is required to satisfy kinematic and dynamic boundary conditions-that is, the normal velocities of the fluid and of the boundary surface must be equal, and the pressure on the free surface must be equal to atmosphere. The single form of the boundary condition can be given by taking a substantial derivative of the dynamic boundary condition and combining the kinematic boundary condition to eliminate the wave elevation term. It is further linearized with respect to the basis double-body flow by neglecting the higher order terms and with Taylor series expanding at the still water plane z = 0, we can have resultant linearized free surface boundary condition for unsteady flow as follows where the parameters shown in Equation (6) are defined as κ = ω 2 g , κ 0 = g U 2 and τ = ωU g . With regard to the boundary condition at the body surface, the kinematic condition, in which the normal component of the fluid velocity is equal to the component of body velocity, is considered in aspect of the space-fixed coordinate. It is further considered that the unsteady force due to the oscillatory movement of the body is balanced by the unsteady flow. The linearization manipulation is taken by the Taylor series expanding with respect to the equilibrium position of the body and dropping higher order terms, and the linearized body boundary condition for the radiation problem is expressed in the form in which (n 1 , n 2 , n 3 ) = → n , (n 4 , n 5 , n 6 ) = → x × → n ; → n denotes the normal vector and points inwards of the hull. The so-called m-term represents the disturbance effects due to the steady flow in unsteady flow problem, and its components are defined by

Hydrodynamic Coefficients
Once the unknown double-body flow and the unsteady flow potentials, due to the oscillatory motion, are solved, the resultant wave elevation on free surface and the pressure on the vessel can be evaluated. The wave elevation and the resultant linear hydrodynamic pressure can be expressed as Equations (9) and (10), respectively.
in which ρ is the fluid density. The pressure integration over the mean wetted surface of the vessel gives the hydrodynamic forces and moments, which can be further classified as the radiation force and wave excitation force. The radiation force, associated with the monochromic oscillatory motions of vessel in six degrees of freedom, is derived from the radiation potential as where the terms a ij and b ij represent the added mass and damping coefficients, respectively (extracted from the real and image part in D i ). Here the subscript ij indicates that the coupled coefficient of added mass or damping acting at the i th degree of freedom is induced by the forced motion of the j th mode. Furthermore, a ij and b ij can be written in the form of the integral as Equations (12) and (13), respectively.

Radiation Boundary Condition
The present study is conducted based on the spline Rankine source scheme with collocation method, and the resultant linear algebraic system is an underdetermined system since the equations raised on one collocation point per panel are insufficient to adequately solve the unknowns of spline coefficients. Thus, the additional end conditions are required to make the system determined and be solvable.
Among end conditions imposed around the extents of the free surface and vessel, the upstream truncation of free surface extent is considered of most dominant and it should be treated carefully. As illustrated in Figure 2, the generated wave systems by the forwardmoving vessel at a two-dimensional aspect can be classified into two major wave groups, which are referred to as the A-wave and B-wave, respectively. A-wave is generated behind the moving vessel and propagates downward, which is grouped by the shorter a 1 -wave and the longer a 2 -wave. By contrast, B-wave propagates at the forward direction, which consists of b 1 -wave and b 2 -wave at different group velocities. in which ρ is the fluid density. The pressure integration over the mean wetted surface of the vessel gives the hydrodynamic forces and moments, which can be further classified as the radiation force and wave excitation force. The radiation force, associated with the monochromic oscillatory motions of vessel in six degrees of freedom, is derived from the radiation potential as (11) where the terms and represent the added mass and damping coefficients, respectively (extracted from the real and image part in Di). Here the subscript ij indicates that the coupled coefficient of added mass or damping acting at the i th degree of freedom is induced by the forced motion of the j th mode. Furthermore, and can be written in the form of the integral as Equations (12) and (13), respectively.

Radiation Boundary Condition
The present study is conducted based on the spline Rankine source scheme with collocation method, and the resultant linear algebraic system is an underdetermined system since the equations raised on one collocation point per panel are insufficient to adequately solve the unknowns of spline coefficients. Thus, the additional end conditions are required to make the system determined and be solvable.
Among end conditions imposed around the extents of the free surface and vessel, the upstream truncation of free surface extent is considered of most dominant and it should be treated carefully. As illustrated in Figure 2, the generated wave systems by the forwardmoving vessel at a two-dimensional aspect can be classified into two major wave groups, which are referred to as the A-wave and B-wave, respectively. A-wave is generated behind the moving vessel and propagates downward, which is grouped by the shorter a1wave and the longer a2-wave. By contrast, B-wave propagates at the forward direction, which consists of b1-wave and b2-wave at different group velocities. The associated group velocities can be determined by the following equations.
It is worth noting that > and concludes that the b2-wave can propagate ahead of the vessel only when undercritical condition < 0.25 . Therefore, this significant The associated group velocities can be determined by the following equations.
It is worth noting that C b 2 g > U and concludes that the b 2 -wave can propagate ahead of the vessel only when undercritical condition τ < 0.25. Therefore, this significant phenomenon of the b 2 -wave leads to demand for appropriate radiation boundary conditions as the end conditions upstream for both overcritical and undercritical conditions.
Based on the authors' previous study [39], the explicit radiation boundary conditions set, derived from the conditions proposed by Seto [40], at upstream truncation are found. However, to avoid the overdetermined system resulted from the triple conditions for undercritical condition in numerical implementation, a common formula is simplified to be as follows, In Equations (16) and (17), k 2 and δ stand for the local wave number and the propagation direction of b 2 -wave through the outward boundary, and both have variation with the local position change. It is noted that k 2 is directly regarded as zero for the overcritical condition.
Due to the induced effect by the body's forward translation, the generated b 2 -wave system consists of the shorter wave ahead of the body and the longer waves behind the body, as illustrated in Figure 3 for the undercritical condition. At the translating coordinate frame coinciding with the body at speed U, when the waves reach point B, the propagation direction rotates from the radial axis by a counterclockwise angle θ.
J. Mar. Sci. Eng. 2022, 10, x FOR PEER REVIEW 7 of 28 phenomenon of the b2-wave leads to demand for appropriate radiation boundary conditions as the end conditions upstream for both overcritical and undercritical conditions. Based on the authors' previous study [39], the explicit radiation boundary conditions set, derived from the conditions proposed by Seto [40], at upstream truncation are found. However, to avoid the overdetermined system resulted from the triple conditions for undercritical condition in numerical implementation, a common formula is simplified to be as follows, In Equations (16) and (17), and stand for the local wave number and the propagation direction of b2-wave through the outward boundary, and both have variation with the local position change. It is noted that is directly regarded as zero for the overcritical condition.
Due to the induced effect by the body's forward translation, the generated b2-wave system consists of the shorter wave ahead of the body and the longer waves behind the body, as illustrated in Figure 3 for the undercritical condition. At the translating coordinate frame coinciding with the body at speed , when the waves reach point B, the propagation direction rotates from the radial axis by a counterclockwise angle . The research of Das and Cheung [29] or Yuan, et al. [30] gives Equation (18) to link the local characteristics of the wave number 2 , the scattered frequency s and the propagation angle (= Θ − ) corresponding to the point in the outward boundary.
Equation (18) is a nonlinear equation set, in which the first equation is derived from The research of Das and Cheung [29] or Yuan, et al. [30] gives Equation (18) to link the local characteristics of the wave number k 2 , the scattered frequency ω s and the propagation angle δ(= Θ − θ) corresponding to the point in the outward boundary.
Equation (18) is a nonlinear equation set, in which the first equation is derived from the law of sines and the condition of AO/U = AB/V p , the second one calculates the local phase velocity V p at B, and the last one is the local dispersion relationship in deep water. An iterative method can provide the local solution of k 2 and δ at the outward boundary point from Equation (18), for satisfying the radiation boundary condition in Equations (16) and (17).

Rayleigh Artificial Damping
As proposed in the proceeding section, the simplified Seto's radiation boundary condition is imposed upstream for all the flow conditions. However, to supplement the insufficiency of imposed upstream boundary conditions, the Rayleigh artificial damping is considered to see its effectiveness, especially for the undercritical conditions (τ < 0.25).
As explained in Nakos [6], introducing Rayleigh viscosity results in an infinitesimal shift of the poles of the linearized free surface boundary condition in complex Fourier plan, and leads to energy dissipation, accordingly. It transforms the encountered oscillatory frequency ω to (ω − iε), in which ε is a positive and small value. Furthermore, we can have Here, the coefficient µ = 2ε is called the Rayleigh artificial damping. The linearized free surface boundary condition in Equation (6) is modified with Rayleigh artificial damping to have expression as follows In the present study, Rayleigh artificial damping µ is assumed to be the spatial varying of the free surface extent for avoiding excessive energy dissipation near the vessel. The formula proposed by Iwashita [32] for µ(τ, R) reads as follows: Here, µ c presents the critical damping level; and f 1 (τ) and f 2 (R) are formulated as the polynomials. The formula of f 1 (τ) is expressed as follows: where t = 1 − (τ − τ s )/(τ e − τ s ) and τ s = 0.4 and τ e = 0.5. The formula of f 2 (R) reads as follows: The coefficients in Equation (23) can be determined by applying the constraints of

Quadratic B-Spline Scheme
In the present study, all values and derivatives of flow velocity potential are expressed in the biquadratic B-spline scheme, in which the associated basis function B (2,2) j is designed in terms of coordinates with local reference frame at the j th quadrilateral panel as the product of two quadratic functions: where ξ = x − x j and η = y − y j , in which x j , y j is the horizontal coordinate of origin of the local coordinate system at the j th panel. In Equation (24), the definition of quadratic The symbol υ (whether for the variable or as the subscript index) in this equation represents either ξ or η. h v presents the panel size along the local ξ-axis or η-axis of the j th panel and provides the nominal length scale. Considering the characteristics of the quadratic basis function, the value and first derivative of the basis function are continuous over the span 3h v . Because the supports span over the three panels in the ξ and η directions, the variation in characteristics in the j th panel depends not only on the spline control coefficient of this panel but also on the spline control coefficients of the eight neighboring panels j k ( k = 1 ∼ 4, 6 ∼ 9), as illustrated in is designed in terms of coordinates with local reference frame at the th quadrilateral panel as the product of two quadratic functions: where = − and = − , in which ( , ) is the horizontal coordinate of origin of the local coordinate system at the th panel. In Equation (24), the definition of quadratic function ( ) ( ) or ( ) ( ) has the expression as follows, The symbol (whether for the variable or as the subscript index) in this equation represents either or . ℎ presents the panel size along the local -axis or -axis of the j th panel and provides the nominal length scale. Considering the characteristics of the quadratic basis function, the value and first derivative of the basis function are continuous over the span 3ℎ . Because the supports span over the three panels in the and directions, the variation in characteristics in the th panel depends not only on the spline control coefficient of this panel but also on the spline control coefficients of the eight neighboring panels ( = 1~4,6~9), as illustrated in Figure 4. Through the present numerical spline scheme, the velocity potential at the point ⃗ on the th panel can be approximated in terms of the highest degree of the two as follows: Through the present numerical spline scheme, the velocity potential at the point → q on the j th panel can be approximated in terms of the highest degree of the two as follows: Here, β  (6) or (20) can be written in the spline form as follows, where γ (n,m−n) j k stands for the spline coefficients of an intermediate ξ n η (m−n) for the normal dipole expression.

Linear System Construction
The linear complex system for solving unknown splined coefficients is conducted by collecting the resultant discrete boundary integral equations at the collocated points and the employment of end conditions at the extent of truncations and boundaries in the present study.
In practice, the computation extents are to be discretized into the quadrilateral panels for numerical implementation. Herein, those quadrilateral panels over the free surface and the vessel's wetted hull are modeled along the I-and J-directions, so we can write N F = N By discretizing the computation extents into (N F + N H ) quadrilateral panels, we can rewrite the boundary integral equation into a discrete form. Let one field point → p be collocated in the l th panel. The discrete boundary integral equation is formulated as follows: where the above-mentioned coefficients are defined below.
In Equations (30) and (31), P (m,n) j and Q (m,n) j are the influence coefficients induced by the Rankine source G and the normal dipole ∂G/∂n over the j th panel Γ j , and defined as follows: P In the evaluation of the influence of the m-terms in right-hand side of Equation (34), the relation holds as follows This is derived from Stokes' theorem with incorporation of the double-body flow as the basis flow in the present study. Details can be found in Nakos [6]. This approach suggests manipulation of the first-order spatial derivatives to replace the second-order spatial derivatives in m-terms to raise evaluation accuracy for the m-terms. Hence the R j is further recast by The collection of the equations of the discrete boundary integral raised at the panel centroids and corresponding end conditions at the boundaries of extents gives the set of simultaneous equations in matrix form as follows where [A] is a square coefficient matrix; [σ] and [B] are the column vectors containing all the unknowns of spline coefficients and force terms of the right-hand side in Equation (28). Once the solution is determined from Equation (37), the flow potential and its normal derivative can be obtained for the subsequent hydrodynamic pressure calculation.

End Conditions
To introduce the appropriate radiation condition for handling the problem of nonunique solution and remedy the underdetermined problem caused by the collocated points on panels, the end conditions are imposed at the truncation of finite extents in free surface and boundary of the vessel's hull.
In the allocation of end conditions, except for the longitudinal truncation at free surface extent, the remaining truncations are imposed with a symmetry condition that makes the flow variation at truncation smooth. As illustrated in Figure 5, the flow variation at truncation point p = 0.5 h ξ , which is the local coordinate relative to the (l + 1) th panel, is symmetric for both sides and therefore the relation of continuous derivative: ∂ϕ ∂ξ (p + ) = ∂ϕ ∂ξ (p − ) → ∂ 2 ϕ ∂ξ 2 = 0, is obtained. At the present spline scheme, we employ the single quadratic spline function, and have σ l − 2σ l+1 + σ l+2 = 0. face extent, the remaining truncations are imposed with a symmetry condition that makes the flow variation at truncation smooth. As illustrated in Figure 5, the flow variation at truncation point = 0.5 ℎ , which is the local coordinate relative to the (l+1) th panel, is symmetric for both sides and therefore the relation of continuous derivative: ( ) = ( ) → = 0, is obtained. At the present spline scheme, we employ the single quadratic spline function, and have σ − 2σ + σ = 0. Regarding the end conditions at the longitudinal truncation, only upstream truncation is considered due to its dominance over flow development. Because the technique of the upwind difference or the collocation point shifting is absent from the present B-spline scheme, the radiation conditions proposed in Equations (16) and (17) are employed upstream, and can be expressed in spline form as Regarding the end conditions at the longitudinal truncation, only upstream truncation is considered due to its dominance over flow development. Because the technique of the upwind difference or the collocation point shifting is absent from the present B-spline scheme, the radiation conditions proposed in Equations (16) and (17) are employed upstream, and can be expressed in spline form as where the associated coefficients L

Solving Characteristics of Local Scattered Waves
An efficient way to obtain local wave characteristics, including the wave number k 2 and the local propagating direction δ required in Equations (16) and (17), is proposed to replace the traditional iteration way of solving Equation (18). Stationary phase method on Bessho three-dimensional translating-pulsating Green function, which is a single integral form satisfying the linearized free-surface boundary condition and radiation condition at infinity, is considered to have the set of parametric equations for the constant-phase curves as Equation (42). The detailed formulation of the Bessho Green function can refer to the works in [41,42].
where ψ(k 2 , ϑ) is the phase function of the wave number k 2 and variable ϑ. Both k 2 and its derivative . k 2 are the functions of ϑ. Here ϑ is the integral variable used in the Bessho Green function for the integral of −π ≤ ϑ ≤ 0 for the undercritical condition (τ < 0.25).
Let (x, y) represents the point B at outward boundary in Figure 3, which is the generated wave pattern by a translating body at undercritical condition. As indicated in the top diagram of Figure 6, it presents the generated wave patterns, including the circular waves (b 2 -waves) and two sets of divergent and transverse waves (a 1 -and a 2 -waves), for the undercritical condition (τ < 0.25). When the circular wave approaches the point B, at which the local wave direction deflects from the radial axis with the angle of θ, due to the Doppler-shifting effect. We can correct the circular wave propagating direction to angle of δ B (= Θ B − θ B ). Θ B corresponds to the polar angle at B, and it can be obtained by the resultant relationship of (κ 0 y/κ 0 x): −π to 0 into the Equations (44) and (45) could have a relationship as the diagram in Figure  6. Once specific , which corresponds to Θ , is pointed out, the associated solution of 2 and succeeds. Compared to the iteration method, the proposed alternative is explicit, and apparently efficient to determine the radiation boundary condition at various upstream points.

Induced Flow by the Submerged Source
The flow field induced by an oscillatory source translating under water has been evaluated by the present Rankine source method with the proposed radiation technique to investigate its feasibility under overcritical and undercritical conditions. The associated computation configuration is presented in Figure 7, in which the source is moving at the speed U and the depth of d. The investigation arises by comparing the present evaluation against the analytic solution obtained using the Bessho form three-dimensional translating-pulsating-source Green function, in which the steepest descent numerical integration proposed by [42] is utilized to obtain the analytic solution in this paper. The propagating direction δ can be obtained from Equation (18) by The consequence from Equations (44) and (45) can be expressed in function of ϑ to construct the relationship between the wave number k 2 and propagating direction δ at local position Θ, as presented in bottom diagram of Figure 6. Substituting ϑ varying from −π to 0 into the Equations (44) and (45) could have a relationship as the diagram in Figure 6. Once specific ϑ, which corresponds to Θ B , is pointed out, the associated solution of k 2 and δ succeeds. Compared to the iteration method, the proposed alternative is explicit, and apparently efficient to determine the radiation boundary condition at various upstream points.

Induced Flow by the Submerged Source
The flow field induced by an oscillatory source translating under water has been evaluated by the present Rankine source method with the proposed radiation technique to investigate its feasibility under overcritical and undercritical conditions. The associated computation configuration is presented in Figure 7, in which the source is moving at the speed U and the depth of d. The investigation arises by comparing the present evaluation against the analytic solution obtained using the Bessho form three-dimensional translatingpulsating-source Green function, in which the steepest descent numerical integration proposed by [42] is utilized to obtain the analytic solution in this paper. Let field point ⃗ be allocated in the l th quadrilateral panel. The associated discrete boundary integral followed by Equation (28) is expressed as in which the right-hand side presents the velocity potential induced by the Rankine source 0 ⃗ at ⃗. Once the spline coefficients in Equation (46) are solved, the unsteady velocity potential over the free surface can be evaluated. The induced unsteady flow field by an oscillating and translating singularity at (0, 0, -L/10) with Froude number of Fr = 0.2 is simulated by the present method. For convenience, only real part of the unsteady velocity potential is considered in the following discussion.
The associated comparison for the evaluated flow field by present Rankine source scheme without Rayleigh damping ( = 0.0) against the analytical solution at the undercritical condition of = 0.2, attributing from the oscillatory frequency at = 1.0, is shown in Figure 8, and the counter solution, introducing the Rayleigh damping ( = 0.6), is presented in Figure 9. The improvement in comparing Figures 8 and 9 can be observed at the pattern of circular waves ahead of the disturbance, and it reveals that the Rayleigh damping supplies the energy dissipation that is insufficient by Seto's radiation boundary condition only. Furthermore, the employed technique of radiation condition with Rayleigh artificial damping yields a complete and accurate flow around the disturbed singularity at this undercritical condition, although the crests are less or more damped. Figure 10 presents the flow patterns by the Rankine source scheme and the analytical solution at = 0.347, which attributes from oscillatory frequency of √3 and Froude number of = 0.2. The associated comparison on the generated divergent wave and transverse wave is at high agreement. Therefore, the good results in Figure 9 ( = 0.2) and Figure 10 ( = 0.347) conclude that the proposed Rankine source scheme, based on Rayleigh artificial damping in additional to simplified Seto's radiation boundary condition, is acceptable for the undercritical and overcritical conditions and can be further applied to the real ship analysis. Let field point → p be allocated in the l th quadrilateral panel. The associated discrete boundary integral followed by Equation (28) is expressed as in which the right-hand side presents the velocity potential induced by the Rankine source → q 0 at → p . Once the spline coefficients in Equation (46) are solved, the unsteady velocity potential ϕ over the free surface can be evaluated. The induced unsteady flow field by an oscillating and translating singularity at (0, 0, -L/10) with Froude number of F r = 0.2 is simulated by the present method. For convenience, only real part of the unsteady velocity potential is considered in the following discussion.
The associated comparison for the evaluated flow field by present Rankine source scheme without Rayleigh damping (µ c = 0.0) against the analytical solution at the undercritical condition of τ = 0.2, attributing from the oscillatory frequency at ω = 1.0, is shown in Figure 8, and the counter solution, introducing the Rayleigh damping (µ c = 0.6), is presented in Figure 9. The improvement in comparing Figures 8 and 9 can be observed at the pattern of circular waves ahead of the disturbance, and it reveals that the Rayleigh damping supplies the energy dissipation that is insufficient by Seto's radiation boundary condition only. Furthermore, the employed technique of radiation condition with Rayleigh artificial damping yields a complete and accurate flow around the disturbed singularity at this undercritical condition, although the crests are less or more damped. Figure 10 presents the flow patterns by the Rankine source scheme and the analytical solution at τ = 0.347, which attributes from oscillatory frequency of √ 3 and Froude number of F r = 0.2. The associated comparison on the generated divergent wave and transverse wave is at high agreement. Therefore, the good results in Figure 9 (τ = 0.2) and Figure 10 (τ = 0.347) conclude that the proposed Rankine source scheme, based on Rayleigh artificial damping in additional to simplified Seto's radiation boundary condition, is acceptable for the undercritical and overcritical conditions and can be further applied to the real ship analysis.        Figure 11, in which the employed covers from 0. and ⁄ is given at 1.5, 2.0 and 2.5. As presented in Figure 11, the value of mean ence obviously decreases with the increasing Rayleigh damping and trends aro level of e = 30%~35%, even at different ⁄ = 1.5, 2.0 or 2.5. The optimal param in this investigation is found to have = 0.6 and ⁄ = 2.0, and the compar evaluated velocity potentials against the analytic one along several longitudinal y/L = 0.06, 0.18, 0.30 and 0.42 are considered. As shown in Figure 12, good agreem tween the evaluated and analytic solutions in both real and image parts is note cially at the area ahead of the disturbance. Therefore, this parametric set of ( then be further employed for solving a radiation flow problem resulting from a re Due to the attempt to have an appropriate R c for Rayleigh artificial damping distribution, we investigated the difference between the solutions by the present Rankine source scheme and analytical method. The mean difference between both solutions is denoted by e, and has the definition in Equation (47), in which φ and ϕ denote the corresponding evaluated and analytical flow potential at arbitrary point → p i in free surface and M points are considered in total. The superscript * stands for the conjugate of complex number. Because the surrounding flow field of disturbance singularity → q 0 is our major concern, the specified region taken into account of the flow solution difference investigation is at a rectangular extent, bounded by L to −2L in longitudinal and from 0 to L in transversal direction, with respect to The induced flow fields at the undercritical condition of τ = 0.2 are solved with various parametric sets of critical Rayleigh damping (µ c ) and critical radial distance (R c ), against the corresponding analytical solution. The resultants of mean difference at various (µ c , R c ) are summarized in Figure 11, in which the employed µ c covers from 0.0 to 0.9, and R c /L is given at 1.5, 2.0 and 2.5. As presented in Figure 11, the value of mean difference obviously decreases with the increasing Rayleigh damping and trends around the level of e = 30%~35%, even at different R c /L = 1.5, 2.0 or 2.5. The optimal parametric set in this investigation is found to have µ c = 0.6 and R c /L = 2.0, and the comparisons of evaluated velocity potentials against the analytic one along several longitudinal lines at y/L = 0.06, 0.18, 0.30 and 0.42 are considered. As shown in Figure 12, good agreement between the evaluated and analytic solutions in both real and image parts is noted, especially at the area ahead of the disturbance. Therefore, this parametric set of (µ c , R c ) can then be further employed for solving a radiation flow problem resulting from a real ship.

Results and Discussion
The attempt to investigate the applicability of present biquadratic B-spline Rankine source scheme on both undercritical and overcritical conditions, the RIOS (Research Initiative on Ocean Ships) bulker is selected as the numerical model for calculations and the results are compared with the experiment data and numerical prediction by Iwashita, et al. [32]. The ship model seakeeping tests of RIOS bulker were conducted in the Research

Results and Discussion
The attempt to investigate the applicability of present biquadratic B-spline Rankine source scheme on both undercritical and overcritical conditions, the RIOS (Research Initiative on Ocean Ships) bulker is selected as the numerical model for calculations and the results are compared with the experiment data and numerical prediction by Iwashita, et al. [32]. The ship model seakeeping tests of RIOS bulker were conducted in the Research

Results and Discussion
The attempt to investigate the applicability of present biquadratic B-spline Rankine source scheme on both undercritical and overcritical conditions, the RIOS (Research Initiative on Ocean Ships) bulker is selected as the numerical model for calculations and the results are compared with the experiment data and numerical prediction by Iwashita, et al. [32]. The ship model seakeeping tests of RIOS bulker were conducted in the Research Institute for Applied Mechanics, Kyushu University and its principal dimensions of the model are shown in Table 1.  (x-axis) x G 0.051 (m) Center of gravity (z-axis) y G −0.020 (m) Gyration radius in pitch K yy /L 0.250 In addition to the capacity of the facility (such as carrier, load cell and so on), the available experiment data range in a seakeeping test is mainly restricted by the size of the test basin, because the near wave field around the ship model might be interfered with by the wave reflection from the side wall of the basin, especially for the longer wave lengths. This wave interference causes inaccuracy in measurement, and it is difficult to avoid for cases at τ ≤ 0.25. Therefore, the existing associated numerical results are taken into account in validation of the author's proposed method. Figure 13 shows the modeled panels representing the hull surface under still water, and 1258 quadrilateral panels in half hull surface are employed in the computation. Under the satisfaction of stability requirement, the total 6498 quadrilateral panels allocated on the free surface computation domain are adopted, which is a rectangular extent bounded by −2.5L ≤ x ≤ 1.5L in x-direction and 0 ≤ y ≤ 4.5L in positive y-direction, as shown in Figure 14. The truncation distance in y-direction is chosen to allow for sufficient space for developing the wave system such that the edges of the sector do not intersect the transverse outward boundary.
The resulting comparison for the coupled hydrodynamic coefficients induced by the oscillating model in heave and pitch modes translating at F r = 0.18 are presented in Figures 15-18, where two axes of abscissa using κL (dimensionless wave number) and τ are adopted for convenient reading. In calculating the results, the oscillatory wave numbers of κL ranging from 0.4 to 40.0 are considered to cover the flow conditions corresponding to τ < 0.25 and τ > 0.25. The critical Rayleigh damping µ c = 0.6 and critical distance R c = 2.0 are adopted in the calculations at various flow conditions. In addition to the capacity of the facility (such as carrier, load cell and so on), the available experiment data range in a seakeeping test is mainly restricted by the size of the test basin, because the near wave field around the ship model might be interfered with by the wave reflection from the side wall of the basin, especially for the longer wave lengths. This wave interference causes inaccuracy in measurement, and it is difficult to avoid for cases at ≤ 0.25. Therefore, the existing associated numerical results are taken into account in validation of the author's proposed method. Figure 13 shows the modeled panels representing the hull surface under still water, and 1258 quadrilateral panels in half hull surface are employed in the computation. Under the satisfaction of stability requirement, the total 6498 quadrilateral panels allocated on the free surface computation domain are adopted, which is a rectangular extent bounded by −2.5 ≤ ≤ 1.5 in x-direction and 0 ≤ ≤ 4.5 in positive y-direction, as shown in Figure 14. The truncation distance in y-direction is chosen to allow for sufficient space for developing the wave system such that the edges of the sector do not intersect the transverse outward boundary. Figure 13. Panels modeled on the model hull [32]. Figure 13. Panels modeled on the model hull [32]. It is worth noting that the explicit radiation boundary condition imposed upstream is an essential implementation for producing wave radiation in the present scheme structure. The calculations with and without damping are shown to investigate the contributions due to Rayleigh artificial damping in additional to inherent radiation boundary conditions.
In Figures 15 and 16, due to the forced heave oscillating motion, the heave-inducedheave and heave-induced-pitch added mass and damping coefficients are presented, respectively. Generally good agreement of the present method in comparison with either the experiments or the other numerical solutions based on the simplified combined method [32] are found. Only some discrepancies occurring in heave-induced-heave added mass, i.e., A33 in Figure 15, at low wave number around = 3.0 and in heave-inducedpitch damping, i.e., B53 in Figure 16, at < 10.0, which indicates that numerical solutions are overestimated against the experiment. Furthermore, the high similarity of trend and quantity reveals that the present method is capable of providing a reasonable solution as the simplified combined method does. The apparent difference between the solutions with and without Rayleigh damping is only found at the region of low wave number, i.e., < 0.25, which indicates Rayleigh artificial damping indeed makes the results more stable than those without imposing Rayleigh artificial damping which generally show large fluctuation phenomena.
The numerical solutions and experimental data for the coupled pitch-induced-pitch and pitch-induced-heave added mass and damping coefficient by the forced pitch oscillating motion are presented in Figures 17 and 18, respectively. Excellent agreement is generally seen in comparison with the experiment, although some overestimated pitch-induced-heave, i.e., B35, at < 10.0 is noted. Because no experiment is available, validation of the present method at low is specifically compared with the existing results from the simplified combined method. As shown in Figure 17, the differences among the pitchinduced-pitch results by numerical methods are small. Furthermore, it is noted that the drop-down phenomenon near = 0.25 in pitch-induced-pitch damping coefficient, i.e., B55, by the present method without damping is improved as smooth one. The effectiveness of the damping is that the coupled pitch-induced-heave results ( Figure 18) is also found, since the oscillation at low is depressed in the result with damping. At = 0.25, a peak in the pitch-induced-heave added mass, i.e., A35, result by the simplified combined method is noted. It is also seen in the results by the present method, although its peak value is higher. Generally, the same conclusions can be drawn for the hydrodynamic coefficients for the forced pitch oscillating motion mode, like the previous observation in forced heave motion mode. From the above comparisons either with experiment and other numerical method [32] in vertical motion modes, the validity of the present method by constructing the resultant boundary integral equation by employing simplified Seto's explicit conditions and Rayleigh artificial damping to accommodate wave radiation has been confirmed and application on the rest motion modes i.e., surge, sway, roll and yaw, may also be regarded to be workable although no related data can be compared. Some benefits on the calculation results by applying the present technique for these modes are shown in Figures 19-26 for reference. The comparisons in Figures 19-22, in which the induced diagonal added mass and damping coefficient at motion modes of surge, sway, roll and yaw are presented, illustrate that introducing Rayleigh artificial damping suppresses the fluctuation phenomena of results without damping in low frequency regions and makes the solutions more stable, i.e., < 0.25. A similar conclusion can also be drawn for the coupled horizontal motion modes, i.e., roll-induced-sway, roll-induced-yaw, sway-induced-roll and yaw-induced-roll, as presented in Figures 23-26, respectively.         It is worth noting that the explicit radiation boundary condition imposed upstream is an essential implementation for producing wave radiation in the present scheme structure. The calculations with and without damping are shown to investigate the contributions due to Rayleigh artificial damping in additional to inherent radiation boundary conditions.
In Figures 15 and 16, due to the forced heave oscillating motion, the heave-inducedheave and heave-induced-pitch added mass and damping coefficients are presented, respectively. Generally good agreement of the present method in comparison with either the experiments or the other numerical solutions based on the simplified combined method [32] are found. Only some discrepancies occurring in heave-induced-heave added mass, i.e., A 33 in Figure 15, at low wave number around κL = 3.0 and in heave-induced-pitch damping, i.e., B 53 in Figure 16, at κL < 10.0, which indicates that numerical solutions are overestimated against the experiment. Furthermore, the high similarity of trend and quantity reveals that the present method is capable of providing a reasonable solution as the simplified combined method does. The apparent difference between the solutions with and without Rayleigh damping is only found at the region of low wave number, i.e., τ < 0.25, which indicates Rayleigh artificial damping indeed makes the results more stable than those without imposing Rayleigh artificial damping which generally show large fluctuation phenomena.
The numerical solutions and experimental data for the coupled pitch-induced-pitch and pitch-induced-heave added mass and damping coefficient by the forced pitch oscillating motion are presented in Figures 17 and 18, respectively. Excellent agreement is generally seen in comparison with the experiment, although some overestimated pitchinduced-heave, i.e., B 35 , at κL < 10.0 is noted. Because no experiment is available, validation of the present method at low τ is specifically compared with the existing results from the simplified combined method. As shown in Figure 17, the differences among the pitch-induced-pitch results by numerical methods are small. Furthermore, it is noted that the drop-down phenomenon near τ = 0.25 in pitch-induced-pitch damping coefficient, i.e., B 55 , by the present method without damping is improved as smooth one. The effectiveness of the damping is that the coupled pitch-induced-heave results ( Figure 18) is also found, since the oscillation at low τ is depressed in the result with damping. At τ = 0.25, a peak in the pitch-induced-heave added mass, i.e., A 35 , result by the simplified combined method is noted. It is also seen in the results by the present method, although its peak value is higher. Generally, the same conclusions can be drawn for the hydrodynamic coefficients for the forced pitch oscillating motion mode, like the previous observation in forced heave motion mode.
From the above comparisons either with experiment and other numerical method [32] in vertical motion modes, the validity of the present method by constructing the resultant boundary integral equation by employing simplified Seto's explicit conditions and Rayleigh artificial damping to accommodate wave radiation has been confirmed and application on the rest motion modes i.e., surge, sway, roll and yaw, may also be regarded to be workable although no related data can be compared. Some benefits on the calculation results by applying the present technique for these modes are shown in Figures 19-26 for reference. The comparisons in Figures 19-22, in which the induced diagonal added mass and damping coefficient at motion modes of surge, sway, roll and yaw are presented, illustrate that introducing Rayleigh artificial damping suppresses the fluctuation phenomena of results without damping in low frequency regions and makes the solutions more stable, i.e., τ < 0.25. A similar conclusion can also be drawn for the coupled horizontal motion modes, i.e., roll-induced-sway, roll-induced-yaw, sway-induced-roll and yaw-induced-roll, as presented in Figures 23-26, respectively.

Conclusions
The radiation of generated waves, which caused by an oscillating and translating body on a free surface, must be modeled accurately in the Rankine source model to avoid the flow field from the influence of wave reflection around the outward boundary. Specifically, under the overcritical condition ( > 0.25), the generated waves propagate backward behind the body, whereas under the undercritical condition ( < 0.25), the generated waves (especially the circular waves) scatter ahead of the body; this different feature indicates that proper numerical means must be chosen when modeling flow on a free surface.
In this paper, a frequency-domain Rankine source method based on a biquadratic Bspline scheme is presented; it involves an improved radiation mechanism and solves the flow problem with no restriction to τ, making it have broader applicability and practical usefulness. Regarding the improved radiation mechanism in the present method, the simplified Seto's radiation boundary conditions are employed upstream and the introduction of spatially varying Rayleigh artificial damping over the free surface is incorporated. Both means are explicit and efficient and do not require control surface allocation in the boundary integral equation.
In evaluations conducted for cases with an oscillating and translating submerged source singularity, the proposed method is compared against the analytic solution obtained from the translating-pulsating-source Green function. The improved radiation boundary mechanism derives accurate wave radiation modelling, especially for the undercritical condition. The hydrodynamic forces induced by the forced oscillating movement of a translating bulker vessel are solved by the present method, and the obtained results, based on with and without introducing Rayleigh artificial damping, have good agreement in comparison with experimental data and public numerical prediction. Further investigation on the results with and without damping illustrates that proposed improved radiation mechanism, i.e., introducing Rayleigh artificial damping in additional

Conclusions
The radiation of generated waves, which caused by an oscillating and translating body on a free surface, must be modeled accurately in the Rankine source model to avoid the flow field from the influence of wave reflection around the outward boundary. Specifically, under the overcritical condition (τ > 0.25), the generated waves propagate backward behind the body, whereas under the undercritical condition (τ < 0.25), the generated waves (especially the circular waves) scatter ahead of the body; this different feature indicates that proper numerical means must be chosen when modeling flow on a free surface.
In this paper, a frequency-domain Rankine source method based on a biquadratic B-spline scheme is presented; it involves an improved radiation mechanism and solves the flow problem with no restriction to τ, making it have broader applicability and practical usefulness. Regarding the improved radiation mechanism in the present method, the simplified Seto's radiation boundary conditions are employed upstream and the introduction of spatially varying Rayleigh artificial damping over the free surface is incorporated. Both means are explicit and efficient and do not require control surface allocation in the boundary integral equation.
In evaluations conducted for cases with an oscillating and translating submerged source singularity, the proposed method is compared against the analytic solution obtained from the translating-pulsating-source Green function. The improved radiation boundary mechanism derives accurate wave radiation modelling, especially for the undercritical condition. The hydrodynamic forces induced by the forced oscillating movement of a translating bulker vessel are solved by the present method, and the obtained results, based on with and without introducing Rayleigh artificial damping, have good agreement in comparison with experimental data and public numerical prediction. Further investigation on the results with and without damping illustrates that proposed improved radiation mechanism, i.e., introducing Rayleigh artificial damping in additional to the simplified Seto's radiation boundary condition, can provide wider feasibility in various flow conditions.
In conclusion, the present B-spline Rankine source method, incorporating simplified Seto's radiation boundary condition and Rayleigh artificial damping, provides accurate flow solutions under various τ conditions, and affords flexibility and computational efficiency. It can thus be useful to naval engineering in offshore and deep-water applications.