THE HYDRODYNAMIC COEFFICIENTS FOR OSCILLATING 2D RECTANGULAR BOX USING WEAKLY COMPRESSIBLE SMOOTHED PARTICLE HYDRODYNAMICS (WCSPH) METHOD

ABSTRACT: An implementation of the weakly compressible smoothed particle hydrodynamics (WCSPH) method is demonstrated to determine the hydrodynamics coefficients through radiation problem of an oscillating 2D rectangular box. Three possible modes of motion namely swaying, heaving, and rolling are carried out to establish the influence of oscillating motions in predicting the added mass and damping. Both solid boundary and fluid flow are modelled by WCSPH and validated by the potential flow and experimental results. Discrepancies observed at lower frequencies are further investigated using different particle resolutions, different time steps, and extending the domain with longer runtime to demonstrate the performance of WCSPH. Finally, flow separation and vortices are discussed and compared with experimental results. 
ABSTRAK: Bagi fenomena yang melibatkan radiasi dalam air, segiempat kotak 2D diosilasikan dengan menggunakan simulasi WCSPH untuk memperoleh pekali hidrodinamik. Mod osilasi terbahagi kepada 3 iaitu sway, heave dan roll. Osilasi dengan mengguna pakai kotak akan mempengaruhi pergerakan air dalam menentukan nilai penambahan jisim dan rendaman. Keseluruhan domain air dan sempadan telah dimodelkan dengan menggunakan WCSPH. Semua model tersebut kemudiannya akan dibandingkan melalui keputusan eksperimen dan teori. Jika keputusan melalui kaedah WCSPH ini berbeza, terutama pada frekuensi rendah, penyelidikan lanjut akan dilakukan dengan menggunakan zarah resolusi yang berbeza, langkah masa yang berbeza dan menambah masa domain ujikaji bagi menilai keputusan WCSPH. Akhirnya, kriteria aliran dan kadar pusaran yang terhasil di sekeliling kotak akan dibincang dan dibandingkan bersama keputusan eksperimen.


INTRODUCTION
Fluid-structure interaction (FSI) is the physical phenomenon that occurs when a fluid force acts on a structure, deforming of moving the structure and this, in turn, changes the boundary conditions of the fluid; which will affect the fluid motion.The interaction has been studied broadly in civil, marine, and coastal engineering.Most approaches are theoretical based on potential flow where the fluid viscosity and nonlinear effect are negligible.With the help of increasing computational power, in many circumstances the FSI can be efficiently described in a large variety of methods ranging from linear theories to fully nonlinear methods.Besides that, combination of partial linear and nonlinear methods are also developed to optimise the results involving nonlinear effects.However, models that are based on potential flow theory are unable to deal with extreme deformation of the free surface in the fluid domain and when viscous or turbulent effects are significant.In order to address the problem, there have been several attempts to solve nonlinear FSI and recent progress in RANS code using either the finite difference method or the finite volume method (FVM) in incorporating both viscous and rotational effects in the flow, making RANS solver methods widely popular [1,2].However, most of these methods are Eulerian and based on the boundary integral method with complicated algorithms that are ineffective in the case of extreme events of waves breaking and water spray.Therefore, Lagrangian meshless methods such as WCSPH method are viewed as alternative in providing accurate numerical solutions to improve inadequacy of mesh-based discretization.
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 dicretized by the interaction between particles rather than grid cells [3,4].Each particle carries an individual mass, velocity, position and any other physical quantity, which evolve over time governed by the governing equation.All particles have a kernel function to define their range of interaction, while the hydrodynamic variables are approximated by integral approximations.The applications of SPH in FSI were started by Monaghan, who performed 2D simulations of wave propagation onto a shallow beach and investigated the entry of a box travelling down a slope into a numerical wave tank (NWT) [5].SPH has also been enhanced with exact information of incompressibility (ISPH) in water-breakwater interaction which uses an implicit pressure update that allows a larger time step but requires more computational work per time step [6,7].Antuono [8] introduced a diffusive term to model the free surface, while a modification in the SPH solver by Colagrassi [9] can correctly simulate dissipation phenomena for viscous flow.GPU-based SPH were also developed by solving shallow water equations in simulating landslide deformations sliding down into NWT and recreating tsunami [10][11][12][13].
This paper extends the previous works of improving the WCSPH method to simulate a propagating free surface in NWT with varies nonlinearities [14].In this study, a 2D fully nonlinear interation between an oscillating rectangular box and a free surface in three modes of motion are simulated using the WCSPH method.

SPH Interpolation
In the basic formulation of SPH, the entire system motion in SPH is discretized into particles.These particles hold individual mass and other properties.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 [15][16], and h is the smoothing length.Integral representation in Eq.
(1) can be written in the form of particle approximation.
where N is the total number of particles and xj is the position vector of particle j within the support domain of x.
For a particle i, Eq. ( 2) can also be concluded as: where mj and ρj are the mass and density of particle j ,respectively, within the support domain of particle i, ij=(xi-xj,h).The derivative of a function for particle i can be written as: 1 ( ) ( ) where the equation follows similar integral representation and particle approximation.
and rij 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 [17].A smoothing length of h=1.3dx is used, where dx is the initial particle spacing.

Force Around Fixed and Floating Bodies
Similar to fluid particles, boundary particles can also be used for simulating rigid bodies for fluid-structure interaction problems [18][19][20].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 dynamic boundary conditions (DBC) [20], 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 [21][22].Hence, boundary particles experience a force per unit mass given by: where FPs denotes fluid particles and   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,   , exerted by the boundary particle k on fluid particle a is the only force computed, but from Eq. ( 9), the force exerted on the moving body can be calculated.By integrating Eq. ( 9) 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 [23].

Computational Setup
The NWT is 16.0 m long and 0.5 m deep, as shown in Fig. 1.A space-fixed Cartesian coordinate system  −  is used, with -axis coincident together with the bottom of the NWT.The origin  is at the intersection of the bottom line and the vertical line, with x positive rightwards and z positive upwards.Beaches are installed at the ends of the wave tank to damp out the generated waves.About halfway down the length of the tank, a rectangular box is placed which is 0.40 m long and 0.40 m wide.The box is placed at points (7.2, 0.5), slightly further away from the fixed wall on the left side of the tank.The box is homogeneous, and its centre of mass coincides with its geometric centre, so half of the box is immersed in water (breadth to draught ratio, B/T = 2).The numerical simulation is performed using approximately 50,000 particles with an initial particle spacing  of 0.010 m.A release time of 2 s allows particle distribution to settle before any movement commences.The motion amplitude is 0.02 m for sway and heave, and 0.10 radians (5.73 deg.) for roll over frequency range between  = 1 rad/s to 9 rad/s as shown in Table 1.

RESULTS AND DISCUSSION
The results of a rectangular section harmonically oscillating at a free surface are presented in this section.The relation between added mass and damping can be represented by Eq. (10).The computation of the hydrodynamic load due to the fluid-structure interaction is carried out using Fourier decomposition to express the total fluid force in terms of a nondimensional complex-valued hydrodynamic function, which real and imaginary parts identify added mass and damping coefficients, respectively.Total force is obtained from the SPH momentum equation by first computing the acceleration vector of each water particle in the first layer of fluid surrounding the rigid body.This value is then multiplied by the fluid particle's mass and reversed in sign so that the resulting vector is sought as the force exerted by the water on the boundary layer [24][25][26].For swaying, only the x component of the acceleration is considered in the computation of the hydrodynamic force.The unit of total force is N/m and is commonly referred to as the 'unit width of force'.Both results from WCSPH and experiment can be non-dimensionalized with (N/m).
Predictions using WCSPH are compared against findings from Vugts [27].Vugts has used 2 different approaches for its predictions of hydrodynamic coefficients, one of which is a theoretical estimation using potential theory and the other uses experiments.For this particular model, the most suitable comparison should be made using its theoretical estimation.The experiments show to what extent the predictions invalidate the theoretical estimation, which is of significant comparison to WCSPH when viscosity and small motion amplitudes are concerned.Figure 3 shows convergence studies of horizontal force,   in sway motion at ω = 7.00 rad/s for different particle refinement.The preliminary tests were done to show convergence of the force computation.All tests are allowed to run for sufficient time for the motion to stabilize and the forces computed to be reliable.As mentioned for swaying motion, hydrodynamic coefficients are evaluated directly from the force around the body in the x direction.While in the case of heaving and rolling motions, computation from WCSPH will result in the total force experienced by the body.Therefore, static force needs to be excluded for comparison with Vugts.Static condition is assumed to be the body in the initial buoyancy with initial vertical force.Static force is subtracted from total force to obtain the dynamic force.

The Sway Test
For each different frequency in sway simulation, the rectangular cross section experiences a horizontal forced oscillation with small influence of vertical forces.The simulation is run for 20 seconds for the fluid flow to become steadier before the calculation of body forces commences.The added mass and damping in sway, together with the The agreement in added mass, ayy and mass coupling coefficients, aφy of WCSPH agrees well with the best section fit line and the experiment at motion at high frequencies.However, few discrepancies are observed at comparatively low frequencies.The experiments done are observed to overestimate the value of byy and underestimate the value of bφy.Hence, the WCSPH results appear to be closer to the theoretical values.In order to overcome this issue, further investigations are carried out by i) extending the domain and ii) implementing smaller particle spacing.The rectangle mark in Fig. 4 denotes an extension of 20 m length while the diamond mark denotes particle spacing of 0.005 m.The implementation of the refined particle spacing and extended domain, however, do not improve the predictions significantly.
Fig. 4: Added mass and damping coefficient in swaying.

The Heave Test
In this test, the rectangular section is restricted so that it oscillates only in a vertical direction, and therefore the influence of a horizontal force exerted by the fluid should be minimal.In both cases, a coarse pressure field can be seen around the rigid boundary between fluid-boundary particles.Boundary near the free surface has asymmetrical integral domains resulting in the neighbour particles' deficiency.
Therefore, when the fluid particles move along the rigid boundary, the densities of the boundary particles fluctuate, resulting in an unstable pressure field around the rigid body.Although the dimension is extended to 20 m, a maximum of 30% discrepancy is estimated for azz.Deviations appear only in the low-frequency range, especially in azz (ω = 1.97 rad/s) where WCSPH share similar discrepancies with the experimental values of Vugts.The influence of breadth-draught (B/T) ratio to azz and bzz is also quite large where it determines the immersed area of the rectangular section for heaving.
Fig. 6: Added-mass and damping coefficient in heaving.

The Roll Test
In this analysis, the rectangular section is rolled at 6 different frequencies by the amplitude of 0.10 radians.Both roll moments and cross coupling forces of sway into roll are investigated.In Fig. 8 and Fig. 9, predictions from WCSPH agree fairly well with experiments and potential theory in added moment of inertia, aφφ and mass coupling, ayφ except for the lower frequency, ω√ (B/2g) = 0.28.Here, nearly all influence is focused on the energy dissipating terms that are coefficients bφφ and byφ.Damping in coupling coefficient, byφ also seems to fit the trend with both experimental and theoretical results.However, it can be observed that the roll damping coefficients, bφφ are over-predicting the potential theory beyond the peak value at ω√ (B/2g) = 0.2 at comparatively high frequencies.There is a large contribution of viscosity in bφφ by WCSPH though they fit the trend line.
It is known that viscous roll damping is one of the dominant damping modes in the roll motion of the forced motion box, representing the effect of vortex shedding of the box with sharp corners.Although sway, heave and particle velocities may also have effect in the vortex damping of the box, these contributions were found to be very small in the case of determining roll damping from forced-roll [28].The vortex damping term directly depends on roll velocity only, which justifies the gradual increase of intensity of nonlinear forces and moments experienced by the box at comparatively higher frequencies.These nonlinear forces are due to viscous effects that can lead to flow separation and generation of vortices.The magnitude of vortices near the immersed sharp corners of the box are presented in Fig. 7.
A similar scenario is also mentioned by Vugts in his experimental measurements.According to him, the viscosity had an effect on the measurement of bφφ, while aφφ also suffers from experimental errors, which underpredict the values from the theoretical counterpart.Hence, we can say that the roll added inertia using WCSPH is consistent with the experimental results for the case ω√ (B/2g) > 0.28.Moreover, both coupling coefficients from sway and rolling share satisfactory results for both added mass and damping coefficients.In this analysis, extended dimension and refinement of particle spacing are not carried out based on previous studies done on sway motion and to prevent any possibility of numerical instability in the fluid domain that could affect the prediction of hydrodynamic coefficients.

CONCLUSION
WCSPH code is developed and employed as a numerical tool to predict the hydrodynamic coefficients in sway, heave, and roll.This is a 2D hydrodynamic test to validate the capability of the WCSPH method in predicting the hydrodynamic coefficients of a radiation problem.Here in WCSPH, the flow is modelled and compared to potential flow and experimental results from Vugts.Some concluding observations from the investigation are given below.• WCSPH showed a good overall agreement in the prediction of hydrodynamic coefficients, namely added mass and damping when compared with theoretical and experimental data.• Convergence analysis was carried out considering most small discrepancies recorded were observed at comparatively low frequencies.However, the results do not improve significantly for both i) refine particle resolution and ii) extending the fluid domain.• Results from damping coefficients particularly in sway and roll damping suggests that viscosity effects do exist in both experiment and WCSPH, which creates flow separation and generation of vortices at sharp corners of the box.• The result might be improved with boundary treatment to smooth out some pressure fluctuation between fluid and rigid particles, which is not discussed in this study.Improvements could also be made using double precision for particle position or variable particle distribution near the vicinity of the rectangular section.

Fig. 2 :
Fig. 2: The cross-section of the rectangular box.
the roll are shown in Fig.4and Fig.5.Both graphs share the same legend.The solid line from Vugts represents calculations based on the best section fit (theoretical value) while the circle denotes the experimental results by Vugts.

Table 1 :
Frequencies used for added mass and damping coefficient analysis