Hydrodynamic Coefficients for a 3-D Uniform Flexible Barge Using Weakly Compressible Smoothed Particle Hydrodynamics

The numerical modelling of the interactions between water waves and floating structures is significant for different areas of the marine sector, especially seakeeping and prediction of wave-induced loads. Seakeeping analysis involving severe flow fluctuations is still quite challenging even for the conventional RANS method. Particle method has been viewed as alternative for such analysis especially those involving deformable boundary, wave breaking and fluid fragmentation around hull shapes. In this paper, the weakly compressible smoothed particle hydrodynamics (WCSPH), a fully Lagrangian particle method, is applied to simulate the symmetric radiation problem for a stationary barge treated as a flexible body. This is carried out by imposing prescribed forced simple harmonic oscillations in heave, pitch and the two- and three-node distortion modes. The resultant, radiation force predictions, namely added mass and fluid damping coefficients, are compared with results from 3-D potential flow boundary element method and 3-D RANS CFD predictions, in order to verify the adopted modelling techniques for WCSPH. WCSPH were found to be in agreement with most results and could predict the fluid actions equally well in most cases.


Introduction
During fluid-structure interaction (FSI), fluid forces acting on the structure result in the structure moving and deforming, which in turn affects the flow boundary conditions, hence fluid motion. This interaction, also known as two-way coupling, has been studied in marine, offshore, civil and coastal engineering. The majority of approaches used are based on potential flow where nonlinear effects are either ignored or allowed for through various assumptions on the body and free surface boundary conditions, leading to a range of partly to fully nonlinear methods (ISSC 2012). However, numerical models based on potential flow theory are unable to deal with extreme free surface deformations, as well as FSI problems where the effects of viscosity, turbulence and compressibility are significant. In order to address such shortcomings, there have been several attempts to solve nonlinear rigid body FSI and benefitting from rapid progress in RANS code development, using either the finite difference method or the finite volume method (FVM) (Weymouth et al. 2005;Wilson et al. 2006;Castiglione et al. 2011;Hochkirch and Mallol 2013;Tezdogan et al. 2015). However, most of these methods are Eulerian and, by and large, ineffective in the case of extreme events of wave breaking and water spray. Therefore, Lagrangian meshless methods are viewed as alternatives in providing accurate numerical solutions to improve inadequacy of mesh-based discretisation.
Smoothed particle hydrodynamics (SPH) is a mesh-free, Lagrangian method whereby the computational domain is represented by a set of interpolation points called particles where the fluid medium is discretised by the interaction between particles rather than grid cells (Shadloo et al. 2012;Chen et al., 2013). Each particle carries an individual mass, velocity, position and any other requisite physical characteristics, which evolve over time through the governing equations. All particles have a kernel function to define their range of interaction, while the hydrodynamic variables are defined by integral approximations.
There are limited number of studies performing seakeeping analysis using particle methods, mainly focusing on extreme events such as slamming and green water. For example, Shibata et al. (2009) used the moving particle semi-implicit (MPS) method to simulate shipping water on a moving ship and validate the impact force on the deck. Slamming events were also modelled where accurate slamming pressures can be estimated using the SPH algorithm (Veen 2010;Veen and Gourlay 2012). More recently, Kawamura et al. (2016) predicts 6-DoF ship motions in severe conditions using a GPU-accelerated SPH simulation.
There have been quite a few recent investigations accounting for structural deformations and two-way coupling using conventional mesh-based RANS methods. For example, El Moctar et al. performed two-way coupling between RANS code and Timoshenko beam model and investigated the effects of springing and whipping (El Moctar et al. 2011); Lakshmynarayanana investigated a containership in regular head waves coupling STAR-CCM+ with the ABAQUS finite element software (Lakshmynarayanana et al. 2015).
The work presented in this paper is the first step in extending the application of WCSPH (weakly compressible SPH) to simulate two-way coupling in FSI. Forced oscillation tests are performed on a uniform flexible barge using WCSPH for the 3-D radiation problem. Hydrodynamic coefficients, namely added mass and damping coefficients, are obtained for the rigid body motions of heave and pitch and the two-node (2VB) and three-node (3VB) symmetric distortion mode shapes, including coupling terms. These are compared with potential flow (using 3-D hydroelasticity) and RANS (using STAR-CCM+) predictions. Domain size, particle numbers and damping zones are modified based on different frequencies of oscillation, allowing the free surface to be well captured by WCSPH.

SPH Interpolation
In the basic formulation of SPH, the entire domain is discretised with particles. These particles hold individual mass, position, velocity, density and any other physical quantity which evolves over time. The approximate integral form of a function at any given position vector of a particle is where x' is another arbitrary position vector in the domain of integration Ω, ω(x-x', h) is a smoothing function (Wendland 2006;Liu 2010), and h is the smoothing length. The integral representation in Eq. (1) can be written in the form of particle approximation, namely where N is the total number of particles and x j is the position vector of particle j within the support domain of x defined by the smoothing function. For a particle i, Eq.
(2) takes the form where m j and ρ j are the mass and density of particle j, respectively, within the support domain of particle i, The derivative of a function for particle i can be written as where the equation follows similar integral representation and particle approximation. In this equation and r ij is the distance between particle i and j. Using Eq. (3) and Eq. (4), the continuity and pressure contribution to the momentum conservation equations can be written as where P is the pressure and F is the acceleration due to gravity. In this paper, a Wendland quantic kernel is used for all WCSPH interpolations (Crespo et al. 2013(Crespo et al. , 2015. A smoothing length of h = 1.3dx is used, where dx is the initial particle spacing. There are two methods of applying incompressibility in flow simulations. The first method, namely incompressible SPH (ISPH), enforced true incompressibility by solving the Poisson's equation for the pressure field. However, due to the complexity of solving the two-step semi-implicit solution process in 3-D, incompressible fluid with a small compressibility, namely WCSPH, is preferred for simpler implementation in flow applications. In WCSPH, the pressure is solved with a state equation which follows the expression, The parameter B is a constant related to the bulk modulus of elasticity of the fluid, ρ 0 is the reference density, usually taken as the density of the fluid at the free surface, and γ is the polytrophic constant, usually between 1 and 7. The second term in Eq. (8), the minus one term, is in order to obtain zero pressure at a surface (Chen et al., 2013). In this paper, γ = 7, B ¼ c 2 0 ρ 0 =γ, with c 0 being the speed of sound at the reference density and c 0 ¼ c ρ 0 ð Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ∂P=∂ρ p .

Implementation
The wall modelled by particles must be treated with care in order to prevent any penetration from inner particles. Two types of boundary conditions were considered in this work, using boundary force particles that repel those in fluid and having layers of fixed wall around the boundaries known as dynamic boundary condition (DBC) (Colagrossi and Landrini 2003;Crespo et al. 2007). In the repulsive boundary condition, the force can take various forms such as Lennard-Jones force, or an empirical function with a singularity so that the force increases as the particle moves nearer to the boundaries. Following this approach, the force can be calculated by where D = 5gH and ρ 1 and ρ 2 are set as 4 and 2 respectively according to Monaghan (1994); r 0 is the cutoff distance normally selected to be the initial particle spacing. The value of f(r) is set to be zero when r > r 0 so that the force is purely repulsive. In the DBC method, boundary particles are forced to have the same properties as corresponding fluid particles. Thus, they follow the same momentum and continuity equation, as well as the equation of state. Their density and pressure are updated but their positions and velocities remained unchanged during the numerical test. They remain fixed in position or move according to an externally imposed position, i.e. wave maker, or imposed velocity (Sun et al. 2016). Symplectic time integration algorithms are time reversible in the absence of friction or viscous effects. This is also often known as the kick-drift-kick scheme, where the kick is the velocity changing according to the force and drift is the coordinate changing with the initial velocity. Symplectic methods preserve geometric features, such as the energy time-reversal symmetry, which are present in Eq. (4), leading to improved resolution of long-term solution behaviour. First, the values of acceleration and density are calculated at the middle of the time step as where the superscript denotes time step and t = nδt. In the second stage, gives the velocity and, hence, position of the particles at the end of the time step, namely At the end of the time step, Dρ for computational cost (Shao and Lo 2003). A variable time step is employed according to (Monaghan 2005) where δt cv is the combination of Courant and the viscous timestep controls and f i is the internal or external force. The values of time step that meet these two criteria are computed, and the smallest one is used for the time step value.

Force Around Fixed and Floating Bodies
Boundary particles, similar to fluid particles, can also be used for simulating rigid body in fluid-structure interaction problems (Colagrossi and Landrini 2003;Adami et al. 2012). The body might drift freely on the free surface with given initial velocity, or it might have a constrained movement along the fluid domain. All boundary particles have similar properties with fluid particles. However, according to DBC (Crespo et al. 2007), a boundary particle is bound to repel approaching fluid particles using repulsive force to prevent any penetration from the fluid particle. Within the same kernel, the force on each boundary particle is computed by adding up the contribution from all the surrounding fluid particles (Lee et al. 2008;Monaghan and Kajtar 2009). Hence, boundary particles experience a force per unit mass given by where FPs denotes fluid particles and f ka is the force per unit mass exerted by fluid particle a on boundary particle k. The force exerted by a fluid particle on each boundary particle follows the principle of equal and opposite action and reaction which is, In the simulations, the repulsive force (DBC), f ak , exerted by the boundary particle k on fluid particle a is the only force computed, and from Eq. (16), the force exerted by fluid particle a on the moving body can be calculated. By integrating Eq. (16) in time, the position of each boundary particle can be determined and moved accordingly. It can be shown that this technique conserves both linear and angular momentum (Monaghan et al. 2003; Tafuni 2016).

Computational Setup
Using the boundary and fluid geometry as described by Kim et al. (2014), a uniform rectangular barge is modelled with main particulars presented in Table 1. Both boundary and fluid domains are discretised using particles whose size and number depend on the initial particle spacing. For the boundary, two layers of boundary particles are adopted to model the conditions of a numerical wave tank (NWT). The domain is modelled in three dimensions where y axis is in the athwartships and x and z axes are along the barge and in the vertical direction, respectively. Unlike FVM model, the WCSPH code in DualSPHysics (Crespo et al. 2015) does not have the option to model only half of the barge about its longitudinal axis to reduce the computational cost (Chen et al. 2013). Therefore, due to the large dimension, several considerations, such as domain size and requisite computational power, need to be made in determining the appropriate particle spacing and particle number. Among them is the number of particles needed particularly in the vicinity of the barge and flow field near the free surface at high frequencies of oscillation. Furthermore, the sensitivity analysis for particle size carried out by Ramli et al. (2015) was also taken into consideration. Particle spacing of 1.0 m with particle numbers shown in Table 2 is considered adequate in capturing details of forces along the oscillating barge. Figure 1 shows all the domain with red denoting the barge comprising of two layers of boundary particles while cross sections of the domain are presented in Fig. 2. All simulations are carried out for a range of forced oscillation frequencies shown in Table 2 with each wavelength, λ, and oscillation frequency. The size of the numerical domain varies with each wavelength.
The NWT is extended either side of the barge in x and y directions that includes two zones, namely wave zone and damping zone. Wave zone allows time for the radiated wave from the oscillating barge to travel before being slowly damped by the beach at the end of the domain. Ideally, the lengths of both these zones should be set equal to the wavelength of the radiated wave. However, considering the number of particles used, wave zone and damping zone for all frequencies are first determined based on the characteristics of the wave with a frequency of 1.0 rad/s. In this case, the wave zone is set double the corresponding wavelength and damping zone is nearly equal to the wavelength. The water depth is kept large enough to avoid shallow water influence to the radiated wave for each individual frequency. However, for the case ω ≤ 1.0 rad/s, the dimension in the longitudinal and athwartship directions could not extend to double the length of the corresponding wavelengths due to the large particle numbers involved. Therefore, the domain sizes are reduced for ω = 0.2 rad/ s and ω = 0.4 rad/s. The reductions are based on some preliminary simulations, and these domain sizes are believed to be adequate in providing accurate results and to avoid precision problems (Ramli et al. 2015). The excitation amplitude, z a , is set to H/2 = 1 m for frequencies ω ≤ 1.0 rad/s and 0.2 m for frequencies ω > 1.0 rad/s. This large value of amplitude is assigned to lower frequency case to ensure that the dynamic force components are not too small compared to the static components.

Computational Parameters
Simulations are carried out with damping beaches set in both the longitudinal and athwartship directions as shown in Fig. 2. The delta-SPH (δ) term is set to 0.1. The Wendland kernel is chosen for the interaction between neighbouring particles, as mentioned before, while the symplectic scheme is used for the step algorithm. Particles are then initialised on a regular grid of particle spacing. The initial SPH smoothing length is then set at h = 1.3dx. WCSPH in DualSPHysics provides for a prescribed motion to be assigned to a body using <mvrectsinu> option for sinusoidal rectilinear movement. However, this prescribed motion can only be applied to heave motion to estimate the value of generalised force along the oscillating barge. In order to implement the flexible mode model, modifications of motions for pitch and vertical two-and three-node of distortion modes are applied directly to the velocity in the source code (Crespo et al. 2015). Following the velocity of a body undergoing a simple harmonic motion, where ż is the velocity, ω is the oscillation frequency (rad/s), z a is the amplitude of motion and t is the time (s). The velocity of each barge particle is imposed as a vector by multiplying the velocity with the eigenvector for individual rigid body motions or distortions. For this uniform barge, the eigenvectors for the distortion modes are calculated using Euler beam theory with free-free end conditions and the corresponding analytical mode shapes. As only the radiation problem is of interest, the structural properties of the barge (i.e. mass and second moment of area) have no effect on the mode shapes. The normalised mode shapes are approximated using polynomials to facilitate their input into WCSPH (Lakshmynarayanana et al. 2015). It should also be noted that for the distortion-related hydrodynamic coefficients, the mode shape is oscillated at the selected frequency of oscillation and not the natural frequency. Implementing Euler beam in this simulation ensures that the displacements change in the z direction. Displacement of each barge particle is updated with new velocity at the end of every time step. These steps can be simplified as follows: where [R y , 0, T z ] is equal to [0, 0, 1] in heave and [1/60, 0, (x − 60)/60] in pitch. Here, R y and T z are the rotational and translational vectors, respectively, and u, v, and w the three velocity components. Time step size for all simulations is kept to 1.0e −4 s, following the highest frequencies tested in Table 2. Time duration of simulation varies depending on individual wavelengths, i.e. 15 s for 2 rad/s and 60 s for 0.2 rad/s with similar output files of every 0.05 s.

Hydrodynamic Coefficients
The computation of the hydrodynamic force, F(t), is carried out using the recorded information on the particle properties at each time step. The force is determined by first computing the acceleration vector of each fluid particle in the vicinity of the oscillating barge and, making use of Newton's second law, by multiplying this value with the mass of each fluid particle. Referring back to Eq. (16), the resulting vector is assumed as the force exerted by the fluid on the barge with opposite sign. Since changes in x and y directions are negligible, only the z component of the acceleration is considered in obtaining the total hydrodynamic force, F(t).  However, in order to compute the dynamic force, the static force has to be subtracted from the total hydrodynamic force. This is similar to the procedure followed by Kim et al. (2014) when using STAR-CCM+. The static forces on the body for rigid body motions and distortions are computed on the mean free surface. In order to obtain the generalised force, F rs (t), the dynamic force for each particle is multiplied with the sth eigenvector where r denotes the index of the motion, including both the rigid body motions and distortions. Using this relation, one can also extract generalised forces for the cross-coupling motion, i.e. diagonal terms of heave-2VB, 2VB-heave, pitch-3VB and 3VBpitch. For example, considering the case of forced pitch motion, Instantaneous values of the hydrodynamic coefficient are obtained by Fourier analysis of time history of the generalised force F rs (t) using discrete window approach for one period, T = 2π/ω, of oscillation, namely

Results and Discussion
This section presents results for the stationary uniform barge oscillating on the free surface at different mode shapes, namely heave, pitch and two-node and three-node vertical bending modes. The predictions from WCSPH are compared to those obtained from 3-D predictions using STAR-CCM+ and 3-D (hydroelasticity) potential flow using the Green's function for pulsating source (Bishop et al. 1986). Kim et al. (2014), in STAR-CCM+, imposed a boundary that oscillates simply harmonically in the shape of the selected modes. Both STAR-CCM+ and WCSPH simulations use inviscid flow, but account for nonlinearities. Each eigenvector relevant to the motion of the barge is determined by approximating the Euler beam solution by a polynomial. It is worthwhile mentioning that results from STAR-CCM+ used for the comparisons can be categorised into two sets. The first set is the initial setup of the numerical model, with mesh size ranging between 1.1 and 4.4 M (lowest frequency), while the second set is obtained by refinement of the mesh around the body and the free surface, resulting in quadrupling the size of the mesh. It should be noted that these refinements were only used for a few low frequencies of oscillation. Figure 3 shows the generalised added mass and damping coefficients for the rigid body motions of heave and pitch. The variables are plotted as a function of the oscillation frequency and are not nondimensionalised. It should be noted that as a result of the symmetry, the cross-coupling coefficients in heave and pitch are zero. Overall, both added mass and fluid damping from WCSPH predictions agree well with STAR-CCM+ and potential flow theory with small discrepancies below 10% in the added mass of heave motion. The potential flow predictions in the vicinity of ω = 1.6 rad/s show the signs of an irregular frequency. Damping predictions show an increasing trend for ω < 0.8 rad/s before starting to decline gradually for shorter waves. Moreover, added mass by WCSPH, which includes up to eight million particles, is in closer agreement with a refined grid STAR-CCM+ at low frequencies in Fig. 3c which comprises approximately ten million cells. In the case of ω = 0.8 rad/s for pitch, small discrepancy is noted compared to the STARCCM+ refined grid in Fig. 3d which suggests that insufficient amount of force predicted around the oscillating barge. The generalised added mass and damping for the two-node (2VB) and three-node (3VB) distortion modes are represented in Fig. 4. For frequencies greater than 0.4 rad/s, 2VB and 3VB coefficients have good overall agreement with both STAR-CCM+ and potential flow. The influence of the irregular frequency can again be observed between ω = 1.6 rad/s and ω = 1.8 rad/s. For 2VB, differences between WCSPH and potential flow theory are smaller for fluid damping than for added mass. Added mass predicted by WCSPH shows a slightly higher value for ω = 0.2 rad/s and ω = 0.4 rad/s. Similar trends are noted in 3VB coefficients where there is a better agreement for added mass between WCSPH and STAR-CCM+. However, at the lowest frequency, the deviation in damping is recorded to be about ten times larger than STAR-CCM+ and potential flow. Primarily, this is due to the complexity of 3VB motion and insufficient length for the wave zone for the radiated wave in the NWT. Additional calculations are also added at ω = 0.9 rad/s. This is to make sure that the peak in fluid damping for 3VB motion is covered. Initially, the size of NWT was planned to be large enough to cover at least three periods of the radiated wave based on each individual wavelength. However, simulating such large domain with uniform particle spacing of 1.0 m would lead to total number of fluid particles in excess of 25 million, which is the limitation of the current graphic card.
The comparisons of hydrodynamic coefficients for the cross-coupling terms with the 2VB distortion modes are shown in Fig. 5. Apart from the predictions at lower frequencies, heave-2VB and 2VB-heave terms show a consistent trend with each other, agreeing well with both STAR-CCM+ and potential flow results. In heave-2VB, the first, relatively small, discrepancy is observed at ω = 0.4 rad/s in added mass and the errors continue to increase for ω = 0.2 rad/s. The error bound recorded is nearly 20% when compared to the potential flow results. While for the 2VB-heave term, large discrepancies are observed at ω = 0.2 rad/s, overpredicting the added mass similar to STAR-CCM+. In this frequency range, the agreement with potential flow results is better for damping compared to STAR-CCM+. Results with a refined grid using STAR-CCM+ are only available for ω = 0.4 rad/s, due to the fact that higher computational power is needed as mentioned before at ω = 0.2 rad/s for both WCSPH and STAR-CCM+, showing that this is an issue of grid refinement. As shown in Eq. (10), it is possible that the contribution of forces exerted Comparison between generalised added mass and damping coefficients. a, b 2VB motion. c, d 3VB motion. Dash line, predictions from potential flow; white circle, predictions from STAR-CCM+, white square, predictions from STAR-CCM+_ fine; and red triangle, predictions from WCSPH by fluid particles onto boundary particles is not accurately computed when the body is oscillating at relatively low frequencies. Figure 6 shows the comparison of hydrodynamic coefficients for the cross-coupling terms relating to the 3VB distortion modes. For relatively high frequencies (ω > 0.8 rad/s), pitch-3VB and 3VB-pitch terms are observed to follow a trend similar to STAR-CCM+ and potential flow although small discrepancies can be observed throughout the frequency range. The deviation in predictions between WCSPH and potential flow results, particularly for added mass in pitch-3VB and 3VB-pitch terms, becomes larger towards lower (a) Added mass-Heave 2VB (b)  (c) Added mass-2VB Heave (d) Damping-2VB Heave Overall, WCSPH overestimated added mass value for most frequencies. This overestimation may explain the complexity in obtaining the generalised force for 3VB motion in comparison to the rigid body motion due to a possibly higher nonlinearity involved when the body oscillates. Refined predictions from STAR-CCM+, denoted by STAR_CCM+_fine, show the sensitivity to mesh density of the force prediction around the oscillating body at relatively low frequencies. The results from the refinement of STAR-CCM+ at ω = 0.4 rad/s suggest that added mass for both terms can be closer to potential flow theory (Kim et al. 2014). Similar improvement may be achieved by WCSPH method if smaller particle spacing is used. However, use of double precision is required in order to use much smaller particle spacing, an additional limitation to the software in question. Furthermore, in Fig. 6b, d, both WCSPH and STAR-CCM+ can be seen to underestimate the value of damping at ω = 0.8 rad/s. The results from STAR-CCM+_fine suggest that refinement would not improve this prediction. The discrepancies computed in 3VB-pitch are observed to be larger than the one computed in pitch-3VB term. One of the reasons could be the difficulty in obtaining generalised force as dynamic components are smaller compared to static components; hence, an increased amplitude maybe needed at lower frequencies.
Overall, WCSPH is able to show consistency of trend in damping for both pitch-3VB and 3VB-pitch terms. This is because both cross-coupling terms should be the same for a stationary symmetrical body, like the barge used in this paper. Forced motions at different time instants by WCSPH for 2VB and 3VB modes are shown in Fig. 7. In Fig. 7, three figures at three different time instants show one complete cycle of the oscillatory forced motion, respectively. Colour scale bars range from 0 to 6 Pa are used to represent pressure distributions along the body. In Fig. 7a, pressure observed to be higher (red colour) where the boundary particles move at their maximum displacement (T = 5.0 s) along with the highest dynamic force, F rs (t). The interaction between boundary and fluid particles is very sensitive for the case where the magnitude of the dynamic force is small in lower frequencies. Figure 8c, d depicts the contour lines and motion of the waves radiated away from the oscillating body for pitch and 3VB respectively at t = 15.0 s. It should be noted that the barge has been excluded from Fig. 8c, d in order to better observe the actions exerted by the fluid particles onto boundary particles. Nevertheless, the position of the barge is the same as in Fig. 8e, f. Results show that fluid particles in WCSPH are radiating away from the oscillating barge in a similar manner obtained in STAR_CCM+, shown in Fig. 8a, b, respectively. In Fig. 8e, f, the colour contour are plotted via MATLAB routines with contribution of fluid particles near free surfaces where the colours show the pressure distributions of the flow field. The distributions are scattered throughout the fluid domain as pressure values on free surface particles do fluctuate. Looking at Fig. 8c, d, red to yellow colour indicates the moment where the barge slams onto the fluid particles, hence creating a high impact of pressure shown in Fig. 8e, f, also shown in the same colour. Observation on the behaviour of the flow field around the oscillating barge is then made by comparing the wave contour between pitch motions in Fig. 8a, c and 3VB motion in Fig. 8b, d by WCSPH and STAR-CCM+ _fine (Kim et al. 2014), respectively. In pitch motion, wave contours in WCSPH show similar wave pattern with STAR_CCM+_fine, except that fluctuations arise near each wall. These fluctuations are the results of computed forces between boundary and fluid particles. The effect could be due to the accuracy of the prediction of the dynamic force. For the 3VB motion, WCSPH shows significant wave pattern around the oscillating barge. Due to the movement in the threenode distortion mode, radiated waves appear to disperse in the vicinity of the barge and mostly cancelling each other before reaching the damping zones. The large domain used in the simulation plays a role in the 3-D effect on the flow field which

Conclusive Remarks
The accuracy of WCSPH in predicting the hydrodynamic coefficients of a stationary barge harmonically oscillating in still water is carried out and verified against predictions using 3-D potential flow and 3-D RANS CFD predictions. The WCSPH model has been successfully applied to the symmetric rigid body motions of heave and pitch and two-node and three-node distortion modes. The comparisons of the added mass and fluid damping coefficients show that WCSPH predictions agree well with STAR-CCM+ and BEM predictions for the vast majority of the cases investigated. It was observed that both STAR-CCM+ and WCSPH display similar features in the predicted hydrodynamic coefficients, including few discrepancies in the range of relatively low frequencies of oscillation. The most likely cause of the aforementioned discrepancies for WCSPH is the inadequacy of the domain size in the longitudinal directions, which requires to be larger (relative to the length of the radiated wave), as well as larger damping zones in preventing any wave reflections. Adopting such large domains may, however, be constrained by normally available computational power. Nevertheless, obtaining accurate predictions at such low frequencies is of academic interest, as they correspond to large wavelengths by comparison to the ship-wave matching region of practical interest for motions and wave-induced loads. Extracting the dynamic force from the total generalised force in such large time window in the case of low frequencies is likely to increase the errors within the calculation of added mass and damping. Moreover, poor performance of WCSPH particularly for the cross-coupling coefficients of pitch-3VB modes and vice versa similarly may be related to the insufficient particle refinement of the total fluid domain, although it has been shown when using STAR-CCM+ that extending the domain does not improve significantly the quality of the results at all frequencies. It has been noted that the STARCCM+ predictions show similar, sometimes worse, trends in the frequency variation of these hydrodynamic coefficients, e.g. in the vicinity of 0.8 rad/s. The pressure distributions and wave contours of the radiated waves obtained from WCSPH are comparable to the refined results by STAR-CCM+.
Simulation of fluid-structure interaction involving a flexible body using the WCSPH method employed in this paper has shown that it is a reliable numerical tool in predicting added mass and fluid damping coefficients. Although this investigation is limited to the radiation problem, the conclusions drawn for this work show that it is applicable to Fig. 8 Wave contours for different mode shapes at ω = 0.8 rad/s. a Pitch and b 3VB motions for STAR-CCM+ (Kim et al. 2014). c Pitch and d 3VB motions for WCSPH. Pressure contours for e pitch and f 3VB motions of WCSPH. These wave and pressure contours are obtained at time t = 15.0 s modelling the behaviour of the 3-D ship in waves. Future work should include double precision to increase accuracy of neighbouring particle tracking and variable resolution technique to improve particle refinement at large domains (Ramli et al. 2015;Dominuque et al. 2014).