Time domain simulation of coupled sloshing–seakeeping problems by SPH–FEM coupling.

<!-- metadata commented in wiki content


==Time domain simulation of coupled sloshing–seakeeping problems by SPH–FEM coupling.==

B. Serván-Camas<sup>1,*</sup>, J.L. Cercós-Pita<sup>2</sup>, J. Colom-Cobb<sup>1</sup>, J. García-Espinosa<sup>1,3</sup>, A. Souto-Iglesias<sup>2</sup>

''<sup>1</sup>''<span style="text-align: center; font-size: 75%;">'' Centre Internacional de Metodes Numerics en Enginyeria (CIMNE), ''[http://www.cimne.upc.edu/ ''http://www.cimne.upc.edu'']</span>

''<sup>2</sup>''<span style="text-align: center; font-size: 75%;">'' CEHINAV Res. Gr., DMFPA, ETSIN, Technical University of Madrid (UPM'')</span>

''<sup>3</sup>''<span style="text-align: center; font-size: 75%;">'' Universitat Politecnica de Catalunya, BarcelonaTech (UPC), Barcelona, Spain''</span>

''<sup>* </sup>''<span style="text-align: center; font-size: 75%;">''Corresponding author: ''[mailto:bservan@cimne.upc.edu,Phone ''bservan@cimne.upc.edu,Phone'']'': +34-934054627''</span>

'''Keywords: '''Seakeeping, sloshing, coupling, SPH, FEM, LNG
-->

==Abstract==

The aim of this work is to carry out numerical simulations in the time domain of seakeeping problems taking into account internal flow in tanks, including sloshing. To this aim, a Smooth Particle Hydrodynamics (SPH) solver for simulating internal flows in tanks is coupled in the time domain to a Finite Element Method (FEM) diffraction–radiation solver developed for seakeeping problems. Validations are carried out comparing against available experimental data. Good agreement between obtained numerical results and experimental data is found.

=1. Introduction=

There are a series of marine operations and/or navigation conditions in which the coupling between internal flow in tanks and seakeeping dynamics can be crucial in configuring the global dynamic response of the vessel. The assessment of such coupling effects may be thus extremely important in order to estimate the viability of certain operations as well as the risks associated to specific load conditions during navigation. Such navigation conditions and/or operations include: those carried out by offshore vessels, equipped with anti-roll tanks, commonly in the oil&gas industry, consisting of deploying piping systems, cables, drilling equipment, etc. (see e.g. Neves et al., 2009); those referred to the offloading of oil or Liquefied natural gas (LNG) from an [https://en.wikipedia.org/wiki/Floating_production_storage_and_offloading floating production storage and offloading] (FPSO) platform, bunkering vessel or a floating LNG (FLNG) to a shuttle, regular vessel or LNG carrier, respectively, in side-by-side configurations (see e.g. Watai et al., 2015); the transport in partially filled tanks of such LNG or its storage in FLNG tanks (see Zhao et al., 2011, for a review on FLNG vessels hydrodynamics).

Due to their large economic cost and much larger penalties in case of any unexpected problem, the planning of above mentioned operations or the definition of filling level thresholds of tanks, requires a deep understanding of the coupled dynamics between the vessel and the internal flows. In order to achieve such understanding it is important to notice that while the external dynamics can be usually resolved in relatively large time scales, the internal dynamics may incorporate violent sloshing flows, which requires resolving complex free-surface dynamics with much shorter time scales. In addition, when wave excitation is present, for low wave amplitudes, the motions are usually linear with amplitude and frequency while for moderate or large wave amplitudes, the response can be extremely nonlinear. For these reasons, in order to simultaneously resolve the seakeeping and internal flow dynamics, different types of solvers are typically used for each of them. For instance, Zhao et al. (2014a) use a frequency domain linear boundary element method (BEM) approach to obtain hydrodynamic coefficients and wave forces, which is coupled in the time–domain with a nonlinear potential solver for the internal dynamics. Another example is the commonly cited reference of Kim et al. (2007) who also resort to frequency domain calculations and resolve the internal flow modeled with Euler equations, with a finite difference scheme. With the same technique for ship motions Li et al. (2012) solved the internal flows with the widespread used open-source volume of fluid (VOF) based solver OpenFOAM.

Regarding the numerical method traditionally adopted for seakeeping simulation, the boundary element method (BEM) has dominated over others like finite element method (FEM). We might find the reason in the fact that most of the computational effort is spent in solving the Laplace equation. Then it might look like BEM offers a lower computational effort. However, Cai et al. (1998) carried out a heuristic study regarding the computational effort required for solving the Laplace equation by BEM and FEM. This study concluded that for a similar three-dimensional problem and a discretization size, the number of unknowns of BEM and FEM are O(N<sup>2</sup>) and O(N<sup>3</sup>), respectively, being N the number of unknowns needed in one dimension to achieve the desired discretization. But the computational costs are O(N<sup>4</sup>) and O(N<sup>3</sup>) respectively. Therefore, while BEM might be more efficient for problems with low number of unknowns, FEM becomes more efficient respect to BEM as the number of unknowns increases. Considering that computational capabilities continuously increases, so does the complexity of problems to be undertaken, and the number of unknowns required. Hence FEM become more efficient than BEM.

Time–domain diffraction–radiation solvers based on the Finite Element Method (FEM) have been previously used (see Ma and Yan, 2009, Serván-Camas and García-Espinosa, 2013). Their capability to naturally incorporate non-linearities is larger than those time–domain solvers based on frequency domain pre-calculations such as the ones used by Kim et al. (2007), Li et al. (2012), or Zhao et al. (2014a). Those solvers have been successfully used with coupled mooring models (Ma and Yan, 2009), in side-by-side (Yan and Ma, 2011) or even in multi-body configurations (Serván-Camas and García-Espinosa, 2013).

FEM based solvers have, however, to the authors’ knowledge, not yet been used to model coupled sloshing and ship motions. Due to the onset of violent sloshing, internal dynamics can be extremely complex and therefore computationally expensive to resolve. In this context, Smoothed Particle Hydrodynamics (SPH) solvers can be a competitive option since, on the one hand, they are able to cope with extremely non-linear and fragmented free-surface flows, and on the other hand, they are easily parallelized with cost competitive graphic card processing units.

<span id='_Ref420340472'></span><span id='_Ref420340474'></span><span id='_Ref420340518'></span>

SPH solvers have been successfully used in the past for solving coupled dynamics of single degree of freedom angular motion coupled sloshing problems (Bouscasse et al., 2014a, 2014b, Bulian et al., 2010) and their parallelization with GPUs has been subject of significant recent attention (see e.g. Crespo et al. 2011; Domingues et al. 2013; Herault et al. 2010; Nishiura et al. 2015). In this paper, the internal flow dynamics is solved using AquaGPUSPH (Cercós-Pita, 2015) and the seakeeping dynamics is solved with SeaFEM (Serván-Camas and García-Espinosa, 2013). Both solvers, AquaGPUSPH and SeaFEM, have been developed by the authors.

The paper is organized as follows: section 2 introduces the internal flow and wave diffraction–radiation solvers; section 3 describes the rigid body dynamics solver; section 4 focuses on the details of the coupling algorithm, as well as in the communication between solvers; section 5 presents a validation study comparing the results obtained in this work against experimental data available for three different case scenarios; finally, section 6 provides the conclusions of this work.

=2. Internal flow and seakeeping solvers=

<span id='_Ref420998627'></span>

Let’s mention that the seakeeping and internal flow solvers are independently developed, compiled, and distributed software. On one hand, the seakeeping solver is SeaFEM, a potential FEM time domain suite of tools for seakeeping simulation [<span id='cite-1'></span>[[#1|1]]] developed by the International Center of Numerical Methods in Engineering (CIMNE) [<span id='cite-2'></span>[[#2|2]]], in collaboration with CompassIS [<span id='cite-_Ref420998627'></span>[[#_Ref420998627|1]]]. On the other hand, the internal flow solver is AQUAgpusph, an open-source 3D SPH solver, accelerated with OpenCL in order to run on GPUs, developed by the Model Basin Research Group [<span id='cite-3'></span>[[#3|3]]] of the Technical University of Madrid.

==2.1. Internal flow solver: AQUAgpusph==

===2.1.1. General.===

AQUAgpusph (Cercós-Pita, 2015) is a recently released free 3D SPH solver, licensed under GPLv3, and accelerated with OpenCL. AQUAgpusph source code along with developers' documentation, a gallery of images, videos, validation cases and a news blog, is available for downloading at [<span id='cite-4'></span>[[#4|4]]]. AQUAgpusph validation is thoroughly documented in Cercós-Pita (2015), where it was used to model the sloshing flows inside a Lead-cooled Fast Nuclear Reactor during an earthquake, focusing on the evaluation of the loads caused by the fluid on the structure. In Cercós-Pita (2015) the solver is validated regarding fluid loads on walls (with emphasis on pressure loads) in dam-break flows and for a TLD case which works as a canonical representation of an anti-roll tank (Bouscasse et al., 2014a, 2014b; Bulian et al., 2010). The internal solver has been also used to simulate the influence of an anti-roll tank in a nonlinear ship motions solver based on pre-calculated (in frequency domain) diffraction forces (Cercos-Pita et al., 2016, 2016b).

===2.1.2. SPH method.===

SPH is a meshless numerical method that was developed in the seventies and first applied in the nineties to free-surface flows (Monaghan, 1994), and is based on the following convolution integral:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left\langle \left. \boldsymbol{f}\left( \boldsymbol{x}\right) \right\rangle \right. =</math><math>\frac{1}{\gamma (\boldsymbol{x})}\int_{\Omega }^{}\boldsymbol{f}\left( \boldsymbol{y}\right) W\left( \left| \boldsymbol{x}-\boldsymbol{y}\right| ,h\right) d\boldsymbol{y}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref431388484'></span>(1)
|}


where <math display="inline">\boldsymbol{f}</math> is a scalar or vectorial function, <math display="inline">W</math> is an even kernel function (In this specific application a Quintic Wendland function is considered), <math display="inline">\Omega</math>  is the fluid domain, and <math display="inline">\gamma (\boldsymbol{x})</math> is a renormalization factor:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\gamma (\boldsymbol{x})=\int_{\Omega }^{}W\left( \left| \boldsymbol{x}-\boldsymbol{y}\right| ,h\right) d\boldsymbol{y}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(2)
|}


The same convolution integral can be applied to a generic differential operator <math display="inline">D</math>,

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left\langle D\left. \boldsymbol{f}\left( \boldsymbol{x}\right) \right\rangle \right. =</math><math>\frac{1}{\gamma (\boldsymbol{x})}\int_{\Omega }^{}D\boldsymbol{f}\left( \boldsymbol{y}\right) W\left( \left| \boldsymbol{x}-\boldsymbol{y}\right| ,h\right) d\boldsymbol{y}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(3)
|}


where <math display="inline">\left\langle D\left. \boldsymbol{f}\left( \boldsymbol{x}\right) \right\rangle \right.</math>  is generally unknown. However, expanding the expression and applying the divergence theorem the differential operator can be conveniently moved to the kernel function:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left\langle D\left. \boldsymbol{f}\left( \boldsymbol{x}\right) \right\rangle \right. =</math><math>\frac{1}{\gamma \left( \boldsymbol{x}\right) }\left( \int_{\Omega }^{}\boldsymbol{f}\left( \boldsymbol{y}\right) \cdot \nabla W\left( \boldsymbol{y}-\boldsymbol{x},h\right) d\boldsymbol{y}+\right. </math><math>\left. \int_{\partial \Omega }^{}\boldsymbol{f}\left( \boldsymbol{y}\right) \cdot \boldsymbol{n}(\boldsymbol{y})W\left( \left| \boldsymbol{x}-\boldsymbol{y}\right| ,h\right) d\boldsymbol{y}\right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref431398521'></span>(4)
|}


In the previous equation <math display="inline">\boldsymbol{n}</math> is the outward oriented normal at the boundary. Traditionally the renormalization factor <math display="inline">\gamma \left( \boldsymbol{x}\right)</math>  and the first integral (so called the boundary integral) have been neglected by the usage of fluid extensions to enforce the boundary conditions. However, in recent years, the possibility of keeping these terms is attracting the interest of the researchers (see Cercos-Pita, 2015; Ferrand et al., 2013; Macià et al., 2012), and, as a matter of fact, they have been retained in the present implementation.

Hence, the SPH convolution (Eq.<span id='cite-_Ref431388484'></span>[[#_Ref431388484|(1)]]) allows the evaluation of differential operators using the already known field values. For practical purposes, the integrals in Eq. <span id='cite-_Ref431398521'></span>[[#_Ref431398521|(4)]] are numerically evaluated by discretizing the fluid domain in a set of particles, transporting the flow properties in a Lagrangian framework. To this end, for a generic particle <math display="inline">a</math>, the differential of volume is computed as <math display="inline">d\boldsymbol{y}=</math><math>{m}_{a}/{\rho }_{a}</math>, where <math display="inline">{m}_{a}</math> is the constant mass of the particle, and <math display="inline">{\rho }_{a}</math> its density. Therefore, the discrete version of the operators read:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\left\langle D\left. \boldsymbol{f}\right\rangle \right. }_{a}=-\frac{1}{{\gamma }_{a}}\left( \sum _{b\in \Omega }^{}{\boldsymbol{f}}_{b}\cdot \nabla W\left( {\boldsymbol{r}}_{b}-\right. \right. </math><math>\left. \left. {\boldsymbol{r}}_{a},h\right) \frac{{m}_{b}}{{\rho }_{b}}+\sum _{b\in \partial \Omega }^{}{\boldsymbol{f}}_{b}\cdot {\boldsymbol{n}}_{b}W\left( \left| {\boldsymbol{r}}_{a}-\right. \right. \right. </math><math>\left. \left. \left. {\boldsymbol{r}}_{b}\right| ,h\right) {s}_{b}\right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref431389593'></span>(5)
|}


with <math display="inline">{s}_{b}</math> being the area of a generic boundary element <math display="inline">b</math>.

===2.1.3. Governing equations.===

The methodology described above can be applied to the Navier–Stokes equations in a weakly compressible Lagrangian formulation:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{d{\rho }_{a}}{dt}=-{\rho }_{a}{\left\langle \nabla \left. \cdot \boldsymbol{u}\right\rangle \right. }_{a}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(6)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{d{\boldsymbol{u}}_{a}}{dt}=-\frac{{\left\langle \nabla \left. p\right\rangle \right. }_{a}}{{\rho }_{a}}+</math><math>\frac{1}{{\rho }_{a}}{\left\langle \nabla \left. \cdot \boldsymbol{T}\right\rangle \right. }_{a}+</math><math>\boldsymbol{g}</math>
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|(7)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{d{\boldsymbol{r}}_{a}}{dt}={\boldsymbol{u}}_{a}</math>
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|(8)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{p}_{a}=\frac{{c}_{0}^{2}{\rho }_{0}}{\Gamma }\left( {\left( \frac{{\rho }_{a}}{{\rho }_{0}}\right) }^{\Gamma }-\right. </math><math>\left. 1\right)</math> 
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|(9)
|}


where <math display="inline">{\boldsymbol{u}}_{a}</math> is the velocity of a generic material particle a, <math display="inline">\boldsymbol{g}</math> is the gravity force, <math display="inline">{\rho }_{0}</math> and <math display="inline">{c}_{0}</math> are the density and the speed of sound of the fluid respectively, <math display="inline">\Gamma</math>  is an incompressibility factor, and <math display="inline">\boldsymbol{T}</math> is the stress tensor.

===2.1.4. Discretization.===

For conservation and consistency reasons (see Colagrossi et al., 2009, 2011, 2013), the above differential operators are not directly approximated using equation Eq. <span id='cite-_Ref431389593'></span>[[#_Ref431389593|(5)]], but are conveniently symmetrized. More specifically, the approximations of the differential operators practised in this work read as follows:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{{\rho }_{a}\left\langle \nabla \left. \cdot \boldsymbol{u}\right\rangle \right. }_{a}=</math><math>\frac{1}{{\gamma }_{a}}\left( \sum _{b\in \Omega }^{}\left( {\boldsymbol{u}}_{b}-{\boldsymbol{u}}_{a}\right) \cdot \nabla {W}_{ab}{m}_{b}+\right. </math><math>\left. \sum _{b\in \partial \Omega }^{}\left( {\boldsymbol{u}}_{b}-{\boldsymbol{u}}_{a}\right) \cdot {\boldsymbol{n}}_{b}{W}_{ab}{\rho }_{b}{s}_{b}\right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref442856299'></span>(10)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{\left\langle \nabla \left. p\right\rangle \right. }_{a}}{{\rho }_{a}}=\frac{1}{{\gamma }_{a}}\left( \sum _{b\in \Omega }^{}\left( \frac{{p}_{b}}{{\rho }_{b}^{2}}+\frac{{p}_{a}}{{\rho }_{a}^{2}}\right) \nabla {W}_{ab}{m}_{b}+\right. </math><math>\left. \sum _{b\in \partial \Omega }^{}\left( \frac{{p}_{b}}{{\rho }_{b}^{2}}+\frac{{p}_{a}}{{\rho }_{a}^{2}}\right) {\boldsymbol{n}}_{b}{W}_{ab}{\rho }_{b}{s}_{b}\right)</math> 
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref442856301'></span>(11)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{1}{{\rho }_{a}}{\left\langle \nabla \left. \cdot \boldsymbol{T}\right\rangle \right. }_{a}=</math><math>\mu \sum _{b\in \Omega }^{}\frac{\left( {\boldsymbol{u}}_{b}-{\boldsymbol{u}}_{a}\right) \cdot \left( {\boldsymbol{r}}_{b}-{\boldsymbol{r}}_{a}\right) }{{\rho }_{a}{\rho }_{b}{\left| {\boldsymbol{r}}_{b}-{\boldsymbol{r}}_{a}\right| }^{2}}\, \nabla {W}_{ab}{m}_{b}</math>
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref442856328'></span>(12)
|}


Although turbulence modeling has not been used, Colagrossi et al. (2011, 2013, 2015) and Bouscasse et al. (2014 a,b) showed that SPH can model fairly accurately turbulent dissipation in this kind of free-surface flows, even in the presence of wave breaking. As a matter of fact, Bouscasse et al. (2014 a,b) used this capability to model with SPH a coupled sloshing and angular motion problem.

===2.1.5. Boundary conditions.===

Regarding boundary conditions, considering the negligible effect of boundary layers on this kind of sloshing flows, free-slip boundary conditions have been imposed. To this aim, boundary integrals, first introduced by Campbell (1989), and formalized later by Kulasegaram et al. (2003), De Leffe et al. (2009), Ferrand et al. (2013) and Macià et al. (2012) have been used. This method is based on the contour integral presented in eq. <span id='cite-_Ref431398521'></span>[[#_Ref431398521|(4)]] and discretized in equations <span id='cite-_Ref442856299'></span>[[#_Ref442856299|(10)]]-<span id='cite-_Ref442856301'></span>[[#_Ref442856301|(11)]]. Since free slip conditions are imposed, boundary does not affect the viscous term computation and is not included in eq. <span id='cite-_Ref442856328'></span>[[#_Ref442856328|(12)]].

To use this technique, the geometrical boundary is discretized in flat area elements whose typical size is that of the inter-particle distance. The main advantages of this technique are (see Cercos-Pita, 2015):<br/>1. It is efficient both in memory consumption and in computational time.

2. It is able to deal with complex geometries.

3. It is fully consistent in the many neighbors limit; contrary to other boundary techniques such as ghost fluids (see Macià et al., 2012 for details).

===2.1.6. Numerical aspects.===

The purely explicit numerical scheme presented above, corresponding to the Weakly-Compressible formulation (Monaghan, 1994), may lead to numerical instabilities unless either extremely small time steps, or additional diffusion, are used. Considering this, the <math display="inline">\delta</math> -SPH diffusive term has been used in this work (see Antuono et al. 2012, 2015), including it in equation <span id='cite-_Ref442856299'></span>[[#_Ref442856299|(10)]] as

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{{\rho }_{a}\left\langle \nabla \left. \cdot \boldsymbol{u}\right\rangle \right. }_{a}=</math><math>\frac{1}{{\gamma }_{a}}\left( \sum _{b\in \Omega }^{}\left( {\boldsymbol{u}}_{b}-{\boldsymbol{u}}_{a}\right) \cdot \nabla {W}_{ab}{m}_{b}+\right. </math><math>\left. \sum _{b\in \partial \Omega }^{}\left( {\boldsymbol{u}}_{b}-{\boldsymbol{u}}_{a}\right) \cdot {\boldsymbol{n}}_{b}{W}_{ab}{\rho }_{b}{s}_{b}\right) \, +</math><math>\delta h{c}_{0}{D\, }_{a}.</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(13)
|}


This diffusive term incorporates a non-dimensional coefficient &#x03b4; whose value is taken as 0.2, an optimal one as discussed in Antuono et al. (2012). The factor <math display="inline">{D\, }_{a}</math> is an approximation to the second derivative of the density field using the following expression

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{D\, }_{a}=\sum _{b\in \Omega }^{}{\psi }_{ba}\frac{\left( {\boldsymbol{r}}_{b}-{\boldsymbol{r}}_{a}\right) \cdot \nabla {W}_{ab}}{{\left( {\boldsymbol{r}}_{b}-{\boldsymbol{r}}_{a}\right) }^{2}}\cdot \frac{{m}_{b}}{{\rho }_{b}}.</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(14)
|}


Different formulations for the factor <math display="inline">{\psi }_{ba}</math> are available in the literature. The one chosen here is that of Antuono et al. (2012), which preserves mass conservation and is exact for hydrostatic fields.

Regarding the time integration scheme, Souto et al. (2006) discussed in detail a modified predictor-corrector, second order in space and first order in velocities and densities, which has been used in present study. The time step has been chosen so that considering a Courant condition on the pressure waves propagation (the most restrictive in this problem) is fulfilled, namely

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\Delta t\, }_{SPH}=0.25\frac{h}{{c}_{0}}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(15)
|}

===2.1.7. Code acceleration and parallelization.===

AQUAgpusph code adds up to existing free graphic processor unit (GPU) accelerated alternatives, namely the Compute Unified Device Architecture (CUDA) based GPUSPH (Herault et al., 2010) and DualSPHysics (Crespo et al., 2011). A performance comparison with these solvers is provided in Cercós-Pita (2015). Main differences with such implementations include: First, AQUAgpusph has been accelerated and parallelized with OpenCL. Second, commonly used boundary condition schemes have been implemented, including the boundary integrals used in present study (Macià, 2012). Third, AQUAgpusph is Python extensible, which allows to customize and couple solid body motions.

The key features that led to use the OpenCL instead of other alternatives like CUDA, OpenMP or MPI (all successfully used in other SPH implementations), can be summarized in the following list: (1) Hardware diversification: OpenCL is not restricted to GPUs, CPUs, or IBM Cell architecture, but is able to deal with all of them. This feature must be taken into consideration when a framework is selected, not only because the architecture may not be the best one, but also because OpenCL is likely easier to adapt to new architectures that will be created in the future. Other SPH developments require 3 different versions of the software: a serial CPU version, a parallel CPU version and a CUDA based version for GPUs. With OpenCL standard, a unique version of the code can cover all these versions, and moreover, it can do it simultaneously, allowing the use of GPUs and CPUs in the same simulation. (2) Highly moddable application: OpenCL is based on the GLSL working method, where the codes that will be executed in the computation device can be loaded and compiled in run time, allowing relying on external files.

==2.2. Potential flow solver==

The seakeeping solver includes a rigid body dynamics solver and a wave diffraction–radiation solver. The rigid body dynamics solver is described in the following section with details on the coupling with AQUAgpusph. Regarding the wave diffraction–radiation solver, we consider the wave problem governed by potential flow theory:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Delta \varphi =0</math> <math>in\, \Omega</math> incompressible and irrotational flow,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(16)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \xi }{\partial t}+\frac{\partial \varphi }{\partial x}\frac{\partial \xi }{\partial x}+</math><math>\frac{\partial \varphi }{\partial y}\frac{\partial \xi }{\partial y}-\frac{\partial \varphi }{\partial z}=</math><math>0</math> <math>on\, z=\xi</math> free surface kinematic boundary condition,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(17)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \varphi }{\partial t}+\frac{1}{2}\nabla \varphi \cdot \nabla \varphi +</math><math>g\xi =0</math> <math>on\, z=\xi</math> free surface dynamic boundary condition,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(18)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{v}}}_{p}\cdot {\mathit{\boldsymbol{n}}}_{p}+{\mathit{\boldsymbol{v}}}_{\varphi }\cdot {\mathit{\boldsymbol{n}}}_{p}=</math><math>0</math> <math>on\, {P\in \Gamma }_{B}</math> body boundary condition,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(19)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \varphi }{\partial z}=0</math> <math>on\, z=-H</math> bottom boundary condition,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(20)
|}


where <math display="inline">\varphi</math>  is the velocity potential, <math display="inline">\xi</math>  is the free surface elevation, <math display="inline">\Omega</math>  is the fluid domain, <math display="inline">{\Gamma }_{B}</math> represents the wetted surface of the present bodies, <math display="inline">{\mathit{\boldsymbol{v}}}_{p}</math> represents the local body velocity on the body surface, <math display="inline">{\mathit{\boldsymbol{v}}}_{\varphi }</math> is the velocity induce by <math display="inline">\varphi</math> , <math display="inline">{\mathit{\boldsymbol{n}}}_{p}</math> is the normal vector at <math display="inline">{\Gamma }_{B}</math>, and  <math display="inline">H</math> is the water depth. And the pressure at any given point P is

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}_{p}=-\rho \frac{\partial \varphi }{\partial t}-\frac{1}{2}\rho \nabla \varphi \cdot \nabla \varphi -</math><math>\rho g{z}_{p}</math> pressure at a point P.
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(21)
|}


Applying Taylor series expansion of the free surface boundary conditions around <math display="inline">z=</math><math>0</math>, introducing a Stokes’ perturbation solution <math display="inline">\varphi =</math><math>{\varphi }^{1}+{\varphi }^{2}+\ldots</math>  and <math display="inline">\xi ={\xi }^{1}+</math><math>{\xi }^{2}+\ldots</math> , and retaining the first order terms, the first order wave problem is governed by

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Delta {\varphi }^{1}=0</math> <math display="inline">in\, \Omega</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(22)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial {\xi }^{1}}{\partial t}-\frac{\partial {\varphi }^{1}}{\partial z}=</math><math>0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(23)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial {\varphi }^{1}}{\partial t}+g{\xi }^{1}=0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(24)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{v}}}_{p}^{1}\cdot {\mathit{\boldsymbol{n}}}_{p}^{0}+{\mathit{\boldsymbol{v}}}_{\varphi }^{1}\cdot {\mathit{\boldsymbol{n}}}_{p}^{0}=</math><math>0</math> <math display="inline">on\, P\in {S}_{B}^{0}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(25)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial {\varphi }^{1}}{\partial z}=0</math> <math display="inline">on\, z=-H</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(26)
|}


and the pressure at any given point P is

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}_{p}^{1}=-\rho \frac{\partial {\varphi }^{1}}{\partial t}-\rho g{z}_{p}</math> pressure at a point P.
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(27)
|}


From this point on, the superscript “1” will be omitted since only first order approach will be discussed. Next, the solution is decomposed as <math display="inline">\varphi =\psi +</math><math>\phi</math>  and <math display="inline">\xi =\zeta +\eta</math> , where <math display="inline">\psi</math>  and <math display="inline">\zeta</math>  are the incident wave velocity potential and free surface elevation respectively, and fulfill

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Delta \psi =0</math> <math display="inline">in\, \Omega</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(28)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \zeta }{\partial t}-\frac{\partial \psi }{\partial z}=0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(29)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \psi }{\partial t}+g\zeta =0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(30)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \psi }{\partial z}=0</math> <math display="inline">on\, z=-H</math>.
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(31)
|}


Hence the incident waves are described by the well solution to Eqs.  <span id='cite-_Ref409686535'></span>[[#_Ref409686535|(36]])-<span id='cite-_Ref441480721'></span>[[#_Ref441480721|(40]]) :

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\psi =\sum _{i}^{}\frac{{A}_{i}g}{{\omega }_{i}}\frac{\mathrm{cosh}\,\left( \left| {\mathit{\boldsymbol{k}}}_{i}\right| \left( H+z\right) \right) }{\mathrm{cosh}\,\left( \left| {\mathit{\boldsymbol{k}}}_{i}\right| H\right) }\mathrm{sin}\,\left( {\mathit{\boldsymbol{k}}}_{i}\mathit{\boldsymbol{x}}-{\omega }_{i}t+{\delta }_{i}\right)</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(32)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\zeta =\sum _{i}^{}{A}_{i}\mathrm{cos}\,\left( {\mathit{\boldsymbol{k}}}_{i}\mathit{\boldsymbol{x}}-{\omega }_{i}t+{\delta }_{i}\right)</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(33)
|}


where <math display="inline">{\mathit{\boldsymbol{k}}}_{\mathit{\boldsymbol{i}}}=</math><math>\frac{2\pi }{{L}_{i}}\left( \mathrm{cos}\,\left( {\gamma }_{i}\right) ,\mathrm{sin}\,\left( {\gamma }_{i}\right) \right)</math> , <math display="inline">{L}_{i}</math> is the wave length, <math display="inline">{\gamma }_{i}</math> is the wave direction of propagation, <math display="inline">{\omega }_{i}</math> is the wave frequency, and <math display="inline">{\delta }_{i}</math> is a phase delay. The following dispersion relation holds:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{\omega }_{i}^{2}=g\left| {\mathit{\boldsymbol{k}}}_{i}\right| \mathrm{tanh}\,\left( \left| {\mathit{\boldsymbol{k}}}_{i}\right| H\right)</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(34)
|}


and the fluid pressure induced by the incident waves at a point P is given by

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{P}_{p\psi }=\sum _{i}^{}\rho g{A}_{i}\frac{\mathrm{cosh}\,\left( \left| {\mathit{\boldsymbol{k}}}_{i}\right| \left( H+z\right) \right) }{\mathrm{cosh}\,\left( \left| {\mathit{\boldsymbol{k}}}_{i}\right| H\right) }\mathrm{cos\, }\,\left( {\mathit{\boldsymbol{k}}}_{i}\mathit{\boldsymbol{x}}-{\omega }_{i}t+{\delta }_{i}\right)</math> .
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(35)
|}


<math display="inline">\phi</math>  and <math display="inline">\eta</math>  are the diffraction–radiation velocity potential and free surface elevation respectively. Then, the  first order wave diffraction–radiation is govern by the following equations:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\nabla }^{2}\phi =0</math> in <math display="inline">\Omega</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref409686535'></span><span id='_Ref409686598'></span>(36)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\partial }_{tt}\phi +g{\partial }_{z}\phi =0</math> on <math display="inline">z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref409686582'></span>(37)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\partial }_{z}\phi =0</math> on <math display="inline">z=-H</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(38)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\nabla \phi ~\mathit{\boldsymbol{\cdot }}{\mathit{\boldsymbol{n}}}_{B}=\left( {\mathit{\boldsymbol{v}}}_{B}-\right. </math><math>\left. \nabla \psi \, \right) \mathit{\boldsymbol{\cdot }}{\mathit{\boldsymbol{n}}}_{B}</math> on  <math display="inline">{\Gamma }_{B}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(39)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\eta =-\frac{1}{g}\partial }_{t}\phi</math> on <math display="inline">z=0</math>.
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref441480721'></span>(40)
|}


And the pressure induce at a point P is

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}_{p}=-\rho \frac{\partial \phi }{\partial t}-\rho g{z}_{p}</math> pressure at a point P.
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(41)
|}


In the previous equations, the domain is assumed to be infinite in the horizontal directions. However, an infinite domain cannot be considered from a numerical point of view. Since waves represented by <math display="inline">\phi</math>  are born at the bodies and propagate in all directions away from the bodies, these waves have to either be dissipated or be let go out of the domain so they will not come back and interact with the bodies. Then, a Sommerfeld radiation condition at the edge of the computational domain will be used:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\partial }_{t}\phi +{c}_{R}\nabla \phi \mathit{\boldsymbol{\cdot }}{\mathit{\boldsymbol{n}}}_{R}=</math><math>0</math> on  <math display="inline">{\Gamma }_{R}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref441070810'></span><span id='_Ref429399623'></span>(42)
|}


where <math display="inline">{\Gamma }_{R}</math> is the surface limiting of the domain in the horizontal directions, <math display="inline">{c}_{R}</math> is a prescribed wave phase velocity, and <math display="inline">{\mathit{\boldsymbol{n}}}_{R}</math> is the normal vector at <math display="inline">{\Gamma }_{R}</math>. The numerical scheme used for Eq. <span id='cite-_Ref441070810'></span>[[#_Ref441070810|(42]]) is as follows:

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{\phi }^{n+1}-{\phi }^{n}}{\Delta t}=\, -{c}_{R}{\left( \nabla \phi \right) }^{\mathit{\boldsymbol{n+1}}}\mathit{\boldsymbol{\cdot }}{\mathit{\boldsymbol{n}}}_{R}</math> on  <math display="inline">{\Gamma }_{R}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(43)
|}


where <math display="inline">{c}_{R}</math> is calculated as the phase velocity of the waves with smaller frequency (larger waves). Eq. <span id='cite-_Ref429399623'></span>[[#_Ref429399623|(42)]] will let waves with phase velocity <math display="inline">{c}_{R}</math> to escape out of the domain. However, waves with different phase velocities (shorter waves) will not be leaving the domain. Then, wave dissipation is inserted into the dynamic free surface boundary condition by varying the pressure at the free surface such that

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}^{n+1}/\rho ={P}_{a}/\rho \, +\kappa \left( \boldsymbol{x}\right) {\left( {\partial }_{z}\phi \right) }^{n+1}</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(44)
|}


where <math display="inline">\kappa \left( \boldsymbol{x}\right)</math>  is a damping coefficient that is set to zero in the near field around the body (no dissipation), and increase gradually, starting at a given distance from the body, in order to dissipate the waves diffracted and radiated by the body.

The associated matrix form to the former governing equations using FEM is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\overline{\overline{\boldsymbol{L}}}\boldsymbol{\phi }={\boldsymbol{b}}^{B}+{\boldsymbol{b}}^{{Z}_{0}\boldsymbol{\, }}+</math><math>{\boldsymbol{b}}^{{Z}_{-H}}+{\boldsymbol{b}}^{R}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(45)
|}


where <math display="inline">\overline{\overline{\boldsymbol{L}}}</math> is the standard laplacian matrix, and <math display="inline">{\boldsymbol{b}}^{B}</math>, <math display="inline">{\boldsymbol{b}}^{{Z}_{0}\boldsymbol{\, }}</math>, <math display="inline">{\boldsymbol{b}}^{{Z}_{-H}}</math> and <math display="inline">{\boldsymbol{b}}^{R}</math> are the vectors obtained by integrating the corresponding boundary condition on  <math display="inline">{\Gamma }_{B}</math>, <math display="inline">z=</math><math>0</math>, <math display="inline">z=-H</math>, and <math display="inline">{\Gamma }_{R}</math>.

A fourth order Padé scheme is used for integrating the free surface boundary condition, and a fourth order finite difference is used for the kinematic boundary condition:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{\phi }^{n+1}-2{\phi }^{n}+{\phi }^{n-1}}{\Delta t}=-g{\partial }_{z}{\phi }^{n}-</math><math>\frac{1}{12}\left( {\partial }_{z}{\phi }^{n+1}-2{\partial }_{z}{\phi }^{n}+{\partial }_{z}{\phi }^{n-1}\right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref409686623'></span>(46)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\eta }^{n+1}=\frac{1}{g\Delta t}\left( \frac{25}{12}{\phi }^{n+1}-4{\phi }^{n}+\right. </math><math>\left. 3{\phi }^{n-1}-\frac{4}{3}{\phi }^{n-2}+\frac{1}{4}{\phi }^{n-3}\right)</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref409686631'></span>(47)
|}


The time step <math display="inline">\Delta t</math> is obtained based on the following stability criteria: <math display="inline">g\Delta {t}^{2}/h<\beta ;</math> where <math display="inline">h</math> is the minimum length in the vertical direction of the free surface discretization elements, and <math display="inline">\beta</math>  a stability limit which is order 1, depending on the discretization used.

More numerical details regarding the body dynamics and wave diffraction–radiation solvers can be found in Serván-Camas and García-Espinosa (2013).

=3. Body dynamics problem=

Under the assumption of small rotations, which is a usual assumption in seakeeping problems, the body dynamics of the floating device is described by the equation of motion (using a frame of reference with origin at the gravity center of the body):

{| class="formulaSCP" style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\overline{\overline{M}}\ddot{\boldsymbol{X}}={\boldsymbol{F}}_{ext}</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref429401588'></span>(48)
|}


where <math display="inline">\overline{\overline{M}}</math> is the body mass matrix, <math display="inline">\ddot{\boldsymbol{X}}</math> is the body acceleration vector, and <math display="inline">{\boldsymbol{F}}_{ext}</math> is the total external forces and moments vector. Inertias, All external loads, including hydrodynamic, hydrostatic, wave, and those coming from the internal flow calculations, are included in <math display="inline">{\boldsymbol{F}}_{ext}</math>. Then, Eq. <span id='cite-_Ref429401588'></span>[[#_Ref429401588|(48)]] is solved using n implicit Newark’s average acceleration method.

<math display="inline">{\boldsymbol{F}}_{ext}</math> is introduced straightforwardly into the body dynamic solver. Hence, the body dynamic solver will account for the non-linearities generated in the internal flow calculations, among others that could be considered.

The dynamics solver is integrated in the diffraction–radiation solver package used to calculate the action of wave loads on the floating body. Then, since the body dynamics solver is implicit and non-linear, it requires iterating in order to achieve convergence within each time step. As shown in Figure 1, there are two iterative loops within each time step. Where <math display="inline">\boldsymbol{X}</math> is the body position vector, <math display="inline">\overset{\cdot}{\boldsymbol{X}}</math> is the body velocity vector, <math display="inline">\ddot{\boldsymbol{X}}</math> is the body acceleration vector, <math display="inline">{\boldsymbol{F}}_{0}</math> are the external forces other than hydrodynamic loads, <math display="inline">{\boldsymbol{F}}_{H}</math> are hydrodynamic loads, and <math display="inline">{\boldsymbol{F}}_{SPH}</math> is the internal tanks vector of forces and moments computed by the SPH solver. With regards to superscripts, “n” represent the current time step, “k” represents the k-th iteration of the solver loop, and “l” represents the l-th iteration of the body dynamics loop.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image1.jpeg|600px]] </div>

<div id="_Ref408927245" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 1: Body dynamics solver algorithm including loads from the internal flow solver.</span></div>

On the one hand, the outer loop iterates over the hydrodynamics solver obtained from the diffraction–radiation problem. On the other hand, the inner loop iterates over the body dynamics solver updating any other loads acting on the body dynamics.

Looking at Figure 1, it can be observed that the diffraction–radiation problem and the body dynamics are solved within the time marching loop, and hence the same time step is used for both solvers. Regarding the SPH solver, it requires a smaller time step, and the exchange of data with the body dynamic solver will be carried out once per time step.

Regarding the body dynamics equations, the <math display="inline">\alpha</math>  Bossak-Newmark’s average acceleration method was implemented (Wood et al., 2005). The iterative procedure of the inner loop requires a relaxation schemes to ensure convergence. Aitkens method was used to relax body accelerations with the aim of accelerating the convergence of the iterative algorithm (Irons and Tuck, 1969). Let <math display="inline">\Delta</math>  be the standard difference operator, and <math display="inline">l</math> the iteration counter, the resulting algorithm of the body dynamics loop can be written as:

:(1)  Set iteration  <math display="inline">l=1</math>.

:(2)  Compute forces and moments acting on body:  <math display="inline">{\boldsymbol{F}}_{ext,l}=</math><math>{\boldsymbol{F}}_{ext}\left( {\ddot{\boldsymbol{X}}}_{r,l}^{t}\right)</math> '''.'''

:(3)  Solve:  ''' ''' <math display="inline">\left( 1-\alpha \right) \overline{\overline{M}}{\ddot{\boldsymbol{X}}}_{l+1}^{t}+</math><math>\alpha \overline{\overline{M}}{\ddot{\boldsymbol{X}}}^{t-\Delta t}={\boldsymbol{F}}_{ext,l}</math>'''.'''

:(4)  Compute relaxation factor:  ''' ''' <math display="inline">{\omega }_{l+1}=1-</math><math>{\omega }_{l}\left[ 1-\frac{\left( \Delta {\ddot{\boldsymbol{X}}}_{r,l}^{t}-\Delta {\ddot{\boldsymbol{X}}}_{r,l+1}^{t}\right) \cdot \Delta {\ddot{\boldsymbol{X}}}_{r,l+1}^{t}}{{\left( \Delta {\ddot{\boldsymbol{X}}}_{r,l}^{t}-\Delta {\ddot{\boldsymbol{X}}}_{r,l+1}^{t}\right) }^{2}}\right]</math> .

:(5)  Relaxation step:  ''' ''' <math display="inline">{\ddot{\boldsymbol{X}}}_{r,l+1}^{t}=</math><math>\left( \boldsymbol{1-}{\omega }_{l+1}\right) {\ddot{\boldsymbol{X}}}_{l}^{t}+\omega {\ddot{\boldsymbol{X}}}_{l+1}^{t}</math>.

:(6)  Check convergence. If not reached advance iteration <math display="inline">l=l+</math><math>1</math> and go to (b).

=4. Coupling scheme=

==4.1. Communication issues==

The seakeeping and internal flow solvers used in this work are independently developed, compiled, and distributed software. Therefore, the development and implementation of a communication strategy between them is a key point of the work.

The general idea of the designed coupling scheme is trivial: body movements calculated by the seakeeping solver are sent to the internal flow solver where forces and moments are evaluated and sent back to the seakeeping solver.

For the exchange of data, the seakeeping solver has a TCL interface, allowing to execute TCL scripts that can access different internal data structures of the solver during the calculation process. For its part, the internal flow solver includes a Python-based interface, developed to interchange data with external sources.

In order to implement the communication scheme, a TCL script [<span id='cite-5'></span>[[#5|5]]] was built, which is responsible for the communication between the two programs. This script is interpreted and executed by the seakeeping solver and communicates with different instances of the internal flow solver. Figure 2 explains the communication scheme used.

The exchange of data is carried out with respect to an instantaneous reference system. Tables 1 and 2 show the data sent from the seakeeping solver to the internal flow solver and vice versa.

The instantaneous reference system at time=t, <math display="inline">R</math>( <math display="inline">t</math>) <math display="inline">{X}_{t}{Y}_{t}{Z}_{t}</math>, is defined such that <math display="inline">R</math>( <math display="inline">0</math>) <math display="inline">{X}_{0}{Y}_{0}{Z}_{0}</math> coincides with the global reference system <math display="inline">OXYZ</math>, <math display="inline">R</math>( <math display="inline">t</math>) describes the trajectory of the point of  the rigid body located initially at the origin of the global system, and axes <math display="inline">{X}_{t}{Y}_{t}{Z}_{t}</math> are parallel to the global axes. To further understand the reference systems used, a simple rigid body movement is presented at Figure 3 showing the displacements and rotations.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 1: data sent from the seakeeping solver to the internal flow</span> <span style="text-align: center; font-size: 75%;">solver</span></div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"| <math>t</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Elapsed simulation time at each time-step</span>
|-
|  style="border: 1pt solid black;text-align: center;"| <math>R\left( t\right) =({X}_{R}\left( t\right) ,{Y}_{R}\left( t\right) ;{Z}_{R}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Coordinates of the body reference point </span> <math display="inline">R</math><span style="text-align: center; font-size: 75%;"> referred to the global system </span> <math display="inline">OXYZ</math>
|-
|  style="border: 1pt solid black;text-align: center;"| <math>\Theta \left( t\right) =({\Theta }_{x}\left( t\right) ,{\Theta }_{y}\left( t\right) ;{\Theta }_{z}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Rotation Angles referred to the global system </span> <math display="inline">OXYZ</math><span style="text-align: center; font-size: 75%;"> of the body</span>
|-
|  style="border: 1pt solid black;text-align: center;"| <math>{V}_{R}\left( t\right) =({V}_{x}\left( t\right) ,{V}_{y}\left( t\right) ;{V}_{z}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Linear velocities of the body reference point </span> <math display="inline">R</math><span style="text-align: center; font-size: 75%;"> referred to the global system </span> <math display="inline">OXYZ</math>
|-
|  style="border: 1pt solid black;text-align: center;"| <math>\omega \left( t\right) =({\omega }_{x}\left( t\right) ,{\omega }_{y}\left( t\right) ;{\omega }_{z}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Angular velocities of the body referred to the global system </span> <math display="inline">OXYZ</math>
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 2: data sent from the internal flow</span>''' '''<span style="text-align: center; font-size: 75%;">solver to the seakeeping solver</span></div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"| <math>F\left( t\right) =({F}_{x}\left( t\right) ,{F}_{y}\left( t\right) ;{F}_{z}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Internal flow forces referred to the global system </span> <math display="inline">OXYZ</math>
|-
|  style="border: 1pt solid black;text-align: center;"| <math>M\left( t\right) =({M}_{x}\left( t\right) ,{M}_{y}\left( t\right) ;{M}_{z}\left( t\right) )</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Internal flow moments referred to the global system </span> <math display="inline">OXYZ</math>
|}


 [[Image:Draft_Garcia_969903716-image2-c.jpeg|600px]]

<div id="_Ref409103438" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 2: Communication scheme.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image3.jpeg|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 3: Rigid body movement and reference systems.</span></div>

==4.2. Coupling algorithm==

The key point in this work is how the coupling between the two solvers is carried out. As aforementioned, both solvers run in the time domain. Therefore, the coupling is based on the exchange of data in Tables 1 and 2 at some specific time steps while both are running. The time step required by the SPH solver ( <math display="inline">\Delta {t}_{SPH}</math>) stability criteria is usually in the order of 10 to 1000 times smaller than the time step used by the seakeeping solver ( <math display="inline">\Delta {t}_{FEM}</math>), depending on the case under study, the grid resolution and number of particles. Then, an explicit staggered approach was selected as the most adequate to couple both solvers. This scheme is summarized next:

:1. Loads calculated by the internal flow solver at <math display="inline">time=</math><math>n\Delta {t}_{FEM}</math> are sent to the seakeeping solver.

:2. The seakeeping solver extrapolates loads from the internal flow solver at <math display="inline">time=(n+</math><math>1)\Delta {t}_{FEM}</math> using a five points Lagrange polynomial.

:3. The seakeeping solver calculates movements at <math display="inline">time=(n+</math><math>1)\Delta {t}_{FEM}</math>.

:4. The seakeeping solver sends movements to the internal flow solver at <math display="inline">time=(n+</math><math>1)\Delta {t}_{FEM}</math>.

:5. The internal flow solver runs from <math display="inline">time=n\Delta {t}_{FEM}</math> to  <math display="inline">time=(n+</math><math>1)\Delta {t}_{FEM}</math> interpolating the body movements sent by the seakeeping solver from <math display="inline">n\Delta {t}_{FEM}</math> to <math display="inline">(n+</math><math>1)\Delta {t}_{FEM}</math>.

Figure 4 provides the details of the coupling algorithm implemented between the seakeeping and internal flow solvers, where <math display="inline">\boldsymbol{F}</math> and <math display="inline">\boldsymbol{M}</math> are the forces and moments sent to the seakeeping solver respectively, <math display="inline">{\boldsymbol{F}}_{\boldsymbol{e}}</math> and <math display="inline">{\boldsymbol{M}}_{\boldsymbol{e}}</math> are the extrapolated forces and moments respectively, <math display="inline">\boldsymbol{R}</math> and <math display="inline">\boldsymbol{\Theta }</math> are the linear and angular movement of the internal tank respectively, and <math display="inline">\boldsymbol{VR}</math> and <math display="inline">\boldsymbol{W}</math> are the linear and angular velocities of the internal tank respectively. <span id='cite-_Ref408927245'></span>[[#_Ref408927245|Figure 1]] also shows how external forces and moments are introduced in the body dynamics solver at the beginning of the time step, and how the movements are sent to the SPH solver once the dynamics solver has converged.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image4.jpeg|498px]] </div>

<div id="_Ref409087331" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 4: Internal flow-seakeeping coupling algorithm.</span></div>

=5. Validation=

==5.1. Barge with water in tanks.==

===5.1.1. Case description===

In order to validate the present coupled internal flow-seakeeping solver, experimental results obtained by Molin et al. (2002) were used to compare with. The experiments consist on the study of the seakeeping response of a barge-like ship, where a tank is extending over a relatively large part of the model. There are two tanks next to each other at the mid-ship whose transverse dimension is close to the model breadth. Figure 5 shows the barge layout with the tanks and Figure 6 shows the unstructured mesh used in the seakeeping solver for the present case study. The mesh size in the floating line is 2.5 cm, and a total of 236,144 tetrahedral elements were used.

 [[Image:Draft_Garcia_969903716-image5.jpeg|600px]]

<div id="_Ref420506294" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 5: Barge with tanks.</span></div>

 [[Image:Draft_Garcia_969903716-image6.jpeg|600px]]

<div id="_Ref409429061" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 6: FEM mesh.</span></div>

===5.1.1 Case 1: same water level in tanks.===

The first case study to analyze is that where the water level is 19cm in both tanks. To achieve the target draft of 10.8cm, additional mass of 40kg was added on the deck of the barge. Table 3 provides the particulars of the barge including the additional mass. Only roll values were provided in Molin et al. (2002). Numerical simulations were carried out with free sway, heave and roll in order to better approximate the experiments. Since both tanks are filled with the same water level, and movements are occurring in the OYZ plane, only one tank was simulated and the second was assumed to behave equally.

Figure 7 represents the experimental results for the case study where both tanks are filled with the same water level (19 cm). The roll rao in the experiments was obtained via spectral analysis from irregular motion tests. The wave spectrum used was a JONSWAP with particulars provided in Table 4.

<div id="_Ref409082707" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 3: Barge main particulars including the additional mass of 40kg on deck.</span></div>

{| style="width: 42%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Breath</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Target draft</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.108 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Displacement</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">275kg</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">XG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">YG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.14 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radii of gyration respect to G: Rxx</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.3704 m</span>
|}


 [[Image:Draft_Garcia_969903716-image7.jpeg|528px]]

<div id="_Ref408996287" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 7: Experimental results obtained in Molin et al. (2002) [34].</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 4: Jonswap spectrum used in experiments.</span></div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Peak enhancement factor </span> <math display="inline">\gamma</math> 
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">2</span>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Significant wave height </span> <math display="inline">{H}_{S}</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">6.6 cm</span>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Peak period </span> <math display="inline">{T}_{p}</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">1.6 s</span>
|}

====5.1.2.1 Monochromatic wave test====

First, a calibration of the roll damping to account for viscous effects is carried out. A monochromatic test with wave amplitude of 2 cm and wave frequency <math display="inline">\omega =</math><math>4.026\, rad/s</math>, which is where the first resonance appears. The test was carried out with 10<sup>5</sup> particles and it was found that 11% of the critical roll damping was necessary to reproduce the RAO value of the experiments, larger than the 8% of the critical damping in Ludvigsen et al. (2013).

Secondly, several tests were carried out to check the effect of increasing the number of particles in the SPH solver. We select a wave frequency <math display="inline">\omega =</math><math>5.984\, rad/s</math>, which corresponds to the second resonance frequency in Figure 7. All cases were simulated including a damping factor of <math display="inline">C=</math><math>0.11\sqrt{4{I}_{xx}{K}_{44}}</math> , where <math display="inline">{I}_{xx}</math> is the roll inertia of the system barge plus weights, and <math display="inline">{K}_{44}</math> is the roll hydrostatic restoring coefficient.

Table 5 provides the roll amplitude obtained for each case. The amplification factor converges, slowly, to a value about 2.7.

<div id="_Ref409426447" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 5: Roll amplitude in rad/m for monochromatic test cases.</span></div>

{| style="width: 39%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Number of particles'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''RAO Roll [rad/m]'''</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.29</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">35,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.45</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">100,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.55</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">250,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.60</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">500,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.64</span>
|}


'''5.1.2.2 Irregular wave''' '' '' '''test'''

Next a RAOs analysis test is carried out. The wave conditions used in the numerical simulation correspond to a JONSWAP spectrum with similar characteristics of that used in the experiments (see Table 4).

The wave amplitude distribution used to mimic the irregular wave conditions of the experiments is provided in the Appendix A. The discretization of the wave spectrum was done using 51 frequencies from 2.5 to rad/s (this is about the frequency resolution used for the experiments). Also a roll damping <math display="inline">C=</math><math>0.11\sqrt{4{I}_{xx}{K}_{44}}</math> was used.

Figure 8 compares the results obtained in the present work using one hundred thousand particles against those from experiments and calculated by WADAM in Ludvigsen et al. (2013) using a potential flow approach to account for the fluid dynamics within the tank. It can be observed that while WADAM provides a smooth curve of results, experiments and the present work show a bit of noise in the obtained curve. We believe that this is because the SPH approach for the fluid flow in the tank is nonlinear, which contributes to energy transfer among frequencies.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image8.jpeg|522px]] </div>

<div id="_Ref409015359" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 8: Results comparison among FEM-SPH coupling (present work), WADAM Ludvigsen et al. (2013) [35], and experiments Molin et al. (2002) [34] for same water level in tanks (19cm+19cm).</span></div>

For frequencies around <math display="inline">\omega =6\, rad/s</math>, experiments show larger values than those obtained in the computational analyses. On the other hand, results obtained in this work are a bit larger than those obtained by WADAM, and they both compare quite well with the experiment for frequencies away from <math display="inline">\omega =</math><math>6\, rad/s</math>.

====5.1.2.3 Computational time====

Table 6 provides some data regarding the computational effort needed for the present approach. The computational time corresponds to those cases presented at Section 5.1.1. It is observed that the calculation of the fluid flow within the tank using the SPH with more than 3.5·10<sup>4</sup> particles is much more expensive than calculating the seakeeping behavior. Furthermore, if the number of particles is bumped up to 5·10<sup>5</sup> then this solver takes more than 98% of the calculation time.

<div id="_Ref409103816" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 6: Computational time</span></div>

{| style="width: 78%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Number'''</span>

<span style="text-align: center; font-size: 75%;">'''of particles'''</span>
|  style="border: 1pt solid black;text-align: center;"| <math display="inline">\boldsymbol{\Delta }{\boldsymbol{t}}_{\boldsymbol{SPH}}</math><span style="text-align: center; font-size: 75%;">''' (ms)'''</span>
|  style="border: 1pt solid black;text-align: center;"| <math display="inline">\boldsymbol{\Delta }{\boldsymbol{t}}_{\boldsymbol{FEM}}</math><span style="text-align: center; font-size: 75%;">''' (ms)'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Simulation'''</span>

<span style="text-align: center; font-size: 75%;">'''time (s)'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Computational'''</span>

<span style="text-align: center; font-size: 75%;">'''time (s)'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''SPH'''</span>

<span style="text-align: center; font-size: 75%;">'''(%)'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''FEM'''</span>

<span style="text-align: center; font-size: 75%;">'''(%)'''</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.17</span> <math display="inline">\boldsymbol{\, \cdot }</math><span style="text-align: center; font-size: 75%;">10<sup>-3</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10<sup>-2</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">30</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1938</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">33.5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">66.5</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">35,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">7.7</span> <math display="inline">\boldsymbol{\, \cdot }</math><span style="text-align: center; font-size: 75%;">10<sup>-4</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10<sup>-2</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">30</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">4109</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">68.4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">31.6</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">100,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">5.4</span> <math display="inline">\boldsymbol{\, \cdot }</math><span style="text-align: center; font-size: 75%;">10<sup>-4</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10<sup>-2</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">30</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">13110</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">88.5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">11.5</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">250,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">4.0</span> <math display="inline">\boldsymbol{\, \cdot }</math><span style="text-align: center; font-size: 75%;">10<sup>-4</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10<sup>-2</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">30</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">41577</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">95.5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">4.5</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">500,000</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.2</span> <math display="inline">\boldsymbol{\, \cdot }</math><span style="text-align: center; font-size: 75%;">10<sup>-4</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10<sup>-2</sup></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">30</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">98616</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">98.2</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.8</span>
|}


For all cases presented in Table 6, the artificial speed of sound for the SPH was set to 10m/s, and the maximum particle velocity during the simulation was in the order of 0.6m/s, which keeps the Mach number under 0.06.

===5.1.2 Case 2: different water level in tanks.===

The second case study to analyze is for the case where the water level in one tank is 19cm and in the other 39cm. No additional mass is necessary to achieve the target draft like in the previous case study. The wave amplitude distribution used to mimic the irregular wave conditions of the experiments is provided in Appendix A. The discretization of the wave spectrum was done using one 101 frequencies ranging from 2.5 to 8 rad/s (frequency resolution was increased respect to the experiments to be able to catch the middle resonance peak shown in the experiments). Also a roll damping <math display="inline">C=</math><math>0.11\sqrt{4{I}_{xx}{K}_{44}}</math> was used. The main particulars of the barge for this case study are given in Table 7.

<div id="_Ref414608973" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 7: Barge main particulars with no additional mass on deck.</span></div>

{| style="width: 41%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Target draft</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.108 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Displacement</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">275kg</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">XG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">YG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.132 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radii of gyration respect to G: Rxx</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.414 m</span>
|}


Figure 9 compares the results obtained in the present work against those from experiments and calculated by WADAM in Ludvigsen et al. (2013). 10<sup>5</sup> Particles were used in the 19cm filling tank, and two hundred thousand particles in the 39cm filling tank. In this case study, three resonance peaks appear. For frequencies around  <math display="inline">\omega =</math><math>4\, rad/s</math>  the largest deviation from numerical to experimental results is found. But, for frequencies that are around <math display="inline">\omega =</math><math>5.25\, rad/s</math> no apparent deviation of the peaks’ frequency is observed and the numerical results around the sharp peak achieve a very similar value to the experimental results. Finally, for frequencies around <math display="inline">\, \omega =</math><math>7\, rad/s</math>, no deviation in frequency is observed, but numerical results are below the predicted peak value.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image9.jpeg|480px]] </div>

<div id="_Ref414609014" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 9: Results comparison among FEM-SPH coupling (present work), WADAM Ludvigsen et al. (2013) [35], and experiments Molin et al. (2002) [34] for different water level in tanks (19cm+39cm)</span></div>

==5.2 Antiroll tank analysis==

Bai and Rhee (1987) provided experimental RAOs, obtained in a model basin, for a supply vessel equipped with an anti-roll tank (ART).  Later on, Kim et. al (2002) and Kim et al. (2007)  carried out numerical simulations of a modified S175 hull in order to compare the results with the experimental data provided in Bai and Rhee (1987). The modified S175 hull used by Kim et. al (2007) was an approximation to the supply vessel used by Bai and Rhee and, despite of it, results obtained by Kim et. al (2007) numerically showed the same trends than the experimental ones.

In this work, we followed Kim et. al (2007) approach using a modified S175 hull to approximate the supply hull used by Bai and Rhee (1987). We used the same aspect ratios of beam to length, and draft to length used by Kim et al. (2007), for the modified S175 hull and the ART. However, no information regarding the exact geometry was found, neither for the position of the gravity center nor for inertias. Therefore, we had to estimate some parameters as well as the hull form.

We started with the S175 hull form provided by the ITTC web site. Then, we proceeded to scale the ship to fulfill the aspect ratios provided by Kim et. al (2007). The resulting geometry is provided in [<span id='cite-6'></span>[[#6|6]]]. The longitudinal positions of the gravity center and ART was estimated as the longitudinal position of the buoyancy center. Transverse positions was set to Y=0 and the vertical position of the gravity center and roll radii of gyration were estimated to approximate roll RAOS results without ART provided by Kim et. al. (2007). Tables 8 and 9 provide the main particulars of the hull form and ART used in this work. Figure 10 shows the modified S175 hull with the ART in place.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 8: Modified S175 main particulars without ART</span></div>

{| style="width: 43%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length between perpendiculars (</span> <math display="inline">L</math><span style="text-align: center; font-size: 75%;">)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">47.6 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Breadth</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">13.7088 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Draft</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.9984 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Displacement</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1498.4 T</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">X<sub>FP</sub></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">23.8 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">X<sub>AP</sub></span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-23.8</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">XG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.73 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">YG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.25 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radius of roll gyration Rxx</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">6.26 m</span>
|}


<div id="_Ref410127170" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 9: Particulars of the ART</span></div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2. 8 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Breadth</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">13.699 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Draft</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.4 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Fresh water filling</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">50 %</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">XT (mid tank)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.73 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">YT (mid tank)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZT (base tank)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-1.8564 m</span>
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
''' [[Image:Draft_Garcia_969903716-image10.jpeg|396px]] '''</div>

<div id="_Ref420592835" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 10: Modified S175 hull with ART.</span></div>

Roll RAOs were obtained for a series of tests with monochromatic waves and were normalized with respect to the maximum wave slope <math display="inline">Ak</math>, where <math display="inline">A</math> is the wave amplitude, and <math display="inline">k</math> the wave number. Also a damping factor of <math display="inline">C=</math><math>\sqrt{4{I}_{xx}{K}_{44}}</math> was used to calibrate the RAO value near resonance.

Figure 11 compares the experimental results obtained by Bai and Rhee (1987) with the numerical results obtained in this work where the results are shown with and without ART. 104 Particles were used in the SPH solver for the ART SPH simulation. The modified S175 was calibrated to show a similar behavior to the supply vessel of Bai and Rhee without ART, as can be observed by comparing the RAOs results without ART. When inserting the ART the results tend to be the same to the ones obtained by Bai and Rhee with a very similar reduction of roll movements.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia_969903716-image11.jpeg|486px]] </div>

<div id="_Ref414610922" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 11: Experimental vs present work. Roll RAOs for monochromatic wave tests with </span> <math display="inline">\boldsymbol{A/L=0.005}</math><span style="text-align: center; font-size: 75%;">. Roll RAOs are obtained with maximum wave slope </span> <math display="inline">\boldsymbol{AK}</math><span style="text-align: center; font-size: 75%;">.</span></div>

Figure 12 presents snapshot of the modified S175 hull with the ART at the resonance wave frequency for the case when there is no ART.  It can be observed that the roll motion is very small, and the fluid flow inside the ART is very nonlinear, exhibiting wave breaking phenomena.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
''' [[Image:Draft_Garcia_969903716-image12.jpeg|600px]] '''</div>

<div id="_Ref414627315" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 12: Snapshots of modified S175 with SRT for </span> <math display="inline">\boldsymbol{\omega }\sqrt{\boldsymbol{L/g}}\boldsymbol{=1.55}</math><span style="text-align: center; font-size: 75%;"> and </span> <math display="inline">\boldsymbol{A/L=0.005}</math><span style="text-align: center; font-size: 75%;">.</span></div>

===5.3 2D vessel including internal flows===

Zhao et al. (2014b) carried out a series of two-dimensional model tests to study the hydrodynamic performance of a floating liquefied natural gas (FLNG) section including internal flow oscillations. The reference FLNG section was ballasted with fresh water and equivalent solid weights respectively, to clarify the coupling effects. RAOs of both motion responses and internal sloshing flows were calculated based on measured data. Particulars of the vessel section, internal tank, and complete model are provided in Tables 10-12, and Figure 13 shows the 3D FEM mesh model used for the simulation.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 10: Particulars of the internal tank</span></div>

{| style="width: 58%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Tank'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Full scale'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''model scale'''</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">40</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.8</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Breadth (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.2</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Water level (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">18</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.36</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mass (kg)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">7.2 E+06</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">57.6</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-4.72</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.0944</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Inertia respect to tank CG Ixx (kg m<sup>2</sup>)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.15 E+09</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.69</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radius of roll gyration Rxx (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">12.7</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.25</span>
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 11: Particulars of the vessel section</span></div>

{| style="width: 62%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Section'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Full scale'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Model scale'''</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">50</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Breadth (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">17.5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.35</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mass (kg)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">4.805 E+06</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">38</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-1.0473</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.021</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Inertia respect to section CG Ixx (kg m<sup>2</sup>)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.027 E+09</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.287</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radius of roll gyration Rxx (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">14.621</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.2924</span>
|}


<div id="_Ref410127171" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 12: Particulars of the complete model</span></div>

{| style="width: 68%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Complete model'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''Full scale'''</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">'''model scale'''</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Draft (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">13.72</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.2744</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mass (kg)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.20 E+07</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">96.04</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">ZG (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-3.25</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">-0.065</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Inertia respect to complete model CG Ixx (kg m<sup>2</sup>)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.22 E+09</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">7.1054</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Radius of roll gyration Rxx (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">13.6</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.272</span>
|}


<span id='_Ref418153621'></span> [[Image:Draft_Garcia_969903716-image13.jpeg|600px]]

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 13: FEM 3D mesh for 2D cross section.</span></div>

A white noise spectrum, with an energy density approximated to the one used in the experiments, was used to carry out the simulations. It consists of 51 frequencies equally spaced between <math display="inline">{\omega }_{min}=</math><math>0.3918\, rad/s</math> and <math display="inline">{\omega }_{max}=</math><math>16.6165\, rad/s</math>, with a wave amplitude of <math display="inline">A=</math><math>0.14\, m</math>. First of all, uncoupled simulations with the equivalent solid weight were carried out to calibrate damping coefficients. Comparing with the solid weight experimental results (see Figure 14, left), damping coefficients were estimated as  <math display="inline">{C}_{33}=</math><math>0.2\sqrt{4M{K}_{33}}</math> for heave, and <math display="inline">{C}_{44}=</math><math>0.055\sqrt{4{I}_{xx}{K}_{44}}</math> for roll. Secondly, coupled simulations were carried out using around 104 particles within the tank. Raos comparisons are provided in Figure 14 (right).

Overall, the results obtained in this work fit well the experimental results provided by Zhao et al. (2014b). Some differences appear regarding the resonance frequency in roll, which might be due to the experimental setup to enforce the 2D condition, which is not fully simulated in the computational model. Besides, the numerical results for the coupled problem shows a peak value in frequencies around the resonance frequency for the uncoupled problem. Although this effect is not present in the experimental results, it was observed in the previous case study with the antiroll tank. So the appearance of this peak is not surprising and the appearance or not is questionable.

==Conclusions==

In this work, a SPH flow solver and a FEM seakeeping solver were coupled to perform simulations of seakeeping dynamics and internal flows in tanks simultaneously. Both solvers run independently in the time–domain, hence a coupling algorithm, including a communication strategy, was developed.

The SPH-FEM coupling showed to be effective for solving seakeeping dynamics coupled with internal flows including sloshing. It was validated for several cases against available experimental data, providing good agreement. For the specific case study of the modified S175 hull with antiroll tank, where sloshing effects were observed, good agreement was found, demonstrating the capability of incorporating highly non-linear phenomena. As a result, the proposed coupling strategy has the capability of performing complex seakeeping analysis coupled with highly non-linear internal flows in reasonable computational times.

This coupled solver will be applied in a future work to solve realistic problems including LNG transportation analysis and complex sloshing phenomena under different sea conditions.

 [[Image:Draft_Garcia_969903716-image14.jpeg|600px]]

<div id="_Ref418154554" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 14: Experimental versus present work. RAOs comparison. Solid case on the left, and liquid (coupled) cases on the right.</span></div>

==Acknowledgments==

The research leading to these results has received funding from the Spanish  Ministry for Economy and Competitiveness under Grants TRA2013-41096-P (OPTIMIZACIÓN DEL TRANSPORTE DE GAS LICUADO EN BUQUES LNG MEDIANTE ESTUDIOS SOBRE INTERACCIÓN FLUIDO-ESTRUCTURA) and ENE2014-59194-C2-2-R (X-SHEAKS: CAMPAÑA EXPERIMENTAL Y VALIDACIÓN). The first author gratefully acknowledges the support of NVIDIA Corporation with the donation of the Tesla K20 GPU used for this research.

==References==

Antuono M., Colagrossi A. and Marrone S., “Numerical diffusive terms in weakly-compressible SPH schemes”, Comput. Phys. Commun. 183, 2012, 570–2580.

Antuono M., Marrone S., Colagrossi A. and Bouscasse B., “Energy balance in the δ-SPH scheme”, Comput. Methods Appl. Mech. Eng. 289, 2015, 209–226.

Bai K.J. and Rhee K.P., “Roll-Damping Tank Test”. Project report, 1987, Seoul National University.

Bouscasse B., Colagrossi A., Souto-Iglesias A. and Cercos-Pita J.L., “Mechanical energy dissipation induced by sloshing and wave breaking in a fully coupled angular motion system”. I. Theoretical formulation and numerical investigation, Phys. Fluids (1994–present) 26, 2014a, 3.

Bouscasse B., Colagrossi A., Souto-Iglesias A. and Cercos-Pita J.L., “Mechanical energy dissipation induced by sloshing and wave breaking in a fully coupled angular motion system”. II. Experimental investigation, Phys. Fluids (1994–present) 26, 2014b, 3.

Bulian G., Souto-Iglesias A., Delorme L. and Botia-Vera E., “SPH simulation of a tuned liquid damper with angular motion”, J. Hydraul. Res. 48 (Extra Issue), 2010, 28–39.

Cai X., Langtangen H.P., Nielsen B.F. and Tveito A., “A finite element method for fully nonlinear water waves”, J. Comput. Phys. 143, 1998, 544–568.

Campbell P.M., Some New Algorithms for Boundary Value Problems in Smooth Particle Hydrodynamics, 1989, Defense Nuclear Agency; USA.

Cercos-Pita J.L., “Aquagpusph, a new free 3D SPH solver accelerated with OpenCL”, Comput. Phys. Commun. 192 (0), 2015, 295–312.

Cercós-Pita, J. L., Bulian, G., Pérez-Rojas, L., and Francescutto, A. (2016). “Coupled simulation of nonlinear ship motions and a free surface tank”. Ocean Engineering, 120:281–288.

Cercós-Pita, J. L. and Bulian, G. (2016b). “Cosimulation of ship motions and sloshing in tanks”. (Under review)

Colagrossi A., Antuono M. and Le Touzé D., “Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model”, Phys. Rev. E (Stat. Nonlinear Soft MatterPhys.) 79 (5),2009, 056701.

Colagrossi A., Antuono M., Souto-Iglesias A. and Le Touzé D., “Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle hydrodynamics formulations in simulating freesurface flows”, Phys. Rev. E 84, 2011, 26705.

Colagrossi A., Souto-Iglesias A., Antuono M. and Marrone S., “Smoothed-particle-hydrodynamics modeling of dissipation mechanisms in gravity waves”, Phys. Rev. E 87, 2013, 023302.

Colagrossi A., Bouscasse B. and Marrone S., “Energy-decomposition analysis for viscous free-surface flows”, Phys. Rev. E 92, 2015, 053003.

Crespo A., Dominguez J., Barreiro A., Gómez-Gesteira M. and Rogers B., “GPUs, a new tool of acceleration in CFD: Efficiency and reliability on smoothed particle hydrodynamics methods”, PLoS ONE 6, 2011, 6.

De Leffe M., Le Touzé D. and Alessandrini B., “Normal flux method at the boundary for SPH”, In: Proceedings of the 4th ERCOFTAC SPHERIC Workshop SPH Appl. 2009, 149–156.

Dominguez J.M., Crespo A.J.C., Valdez-Balderas D., Rogers B.D. and Gómez-Gesteira M., “New multi-GPU implementation for smoothed particle hydrodynamics on heterogeneous clusters”, Comput. Phys. Commun. 184, 2013, 1848–1860.

Ferrand M., Laurence D.R., Rogers B.D., Violeau D. and Kassiotis C., “Unified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method” Int. J. Numer. Methods Fluids 71, 2013, 446–472.

Herault A., Bilotta G. and Dalrymple R.A., “SPH on GPU with CUDA”, J. Hydraul. Res. 48 (Extra Issue), 2010, 74–79.

Irons B.M. and Tuck R.C., “A version of the Aitken accelerator for computer iteration”, Int. J. Numer. Methods Eng. 1, 1969, 275–277.

Li Y., Zhu R., Miao R. and Fan, J G., “Simulation of tank sloshing based on OpenFOAM and coupling with ship motions in time domain”, J. Hydrodyn., Ser. B 24 (3), 2012, 450–457.

Kim Y., “A numerical study on sloshing flows coupled with ship motion—the anti-rolling tank problem”, J. Ship Res. 46 (1), 2002, 52–62.

Kim Y., Nam B.W., Kim D.W. and Kim Y.S., “Study on coupling effects of ship motion and sloshing”, Ocean Eng. 34, 2007, 2176–2187.

Kulasegaram S., Bonet J., Lewis W.R. and Profit M., “A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications”, Comput. Mech. 33 (4), 2003, 316–325.

Ludvigsen, A., Pan, Z.Y., Gou, P., Vada, T., 2013. “Adapting a linear potential theory solver for the outer hull to account for fluid dynamics in tanks”. In: Proceedings of the ASME 2013 32nd International Conference on Ocean, Offshore and Arctic Engineering OMAE2013 June 9–14, Nantes, France.

Ma Q.W. and Yan S., “QALE-FEM for numerical modelling of non-linear interaction between 3D moored floating bodies and steep waves”, Int. J. Numer. Methods Eng. 78 (6), 2009, 713–756.

Macià F., Antuono M., González L.M. and Colagrossi A., “Theoretical analysis of the no-slip boundary condition enforcement in SPH methods”, Progress. Theor. Phys. 125 (6), 2011, 1091–1121.

Macià F., González L.M., Cercos-Pita J.L. and Souto-Iglesias A., “A boundary integral SPH formulation”. Consistency and applications to ISPH and WCSPH, Progress. Theor. Phys. 128 (3), 2012, 439–462.

Molin, B., Remy, F., Rigaud, S., de Jouette Ch., 2002. “LNG-FPSO’s: frequency domain coupled analysis of support and liquid cargo motions”. In: Proceedings of INAM Conference, Rethymnon, Greece.

Monaghan J.J., “Simulating free surface flows with SPH”, J. Comput. Phys. 110, 1994, 399–406.

Neves M.A.S., Merino J.A. and Rodríguez C.A., “A nonlinear model of parametric rolling stabilization by anti-roll tanks”, Ocean Eng. 36 (14), 2009, 1048–1059.

Nishiura D., Furuichi M. and Sakaguchi H., “Computational performance of a smoothed particle hydrodynamics simulation for shared-memory parallel computing”, Comput. Phys. Commun. 194 (0), 2015, 18–32.

Serván-Camas B. and García-Espinosa J., “Accelerated 3D multi-body seakeeping simulations using unstructured finite elements”, J. Comput. Phys. 252, 2013, 382–403.

Souto-Iglesias A., Delorme L., Rojas L.P. and Abril S., “Liquid moment amplitude assessment in sloshing type problems with SPH”, Ocean Eng. 33, 2006, 11–12.

Watai R.A., Dinoi P., Ruggeri F., Souto-Iglesias A. and Simos A.N., “Rankine time–domain method with application to side-by-side gap flow modeling”, Appl. Ocean Res. 50 (0), 2015, 69–90.

Wood W.L., Bossak M. and Zienkiewicz O.C., “An alpha modification of Newmark's method”, Int. J. Numer. Methods Eng. 5 (10), 2005, 1562–1566.

Yan, S., Ma, Q.W., Cheng, X., 2011. “Fully nonlinear simulation of two floating structures in close proximity subjected to oblique waves”. In: Proceedings of the 21st International Offshore and Polar Engineering Conference (ISOPE). The International Society of Offshore and Polar Engineers (ISOPE).

Zhao W.H., Yang J., Hu Z.Q. and Wei Y.F., “Recent developments on the hydrodynamics of floating liquid natural gas (FNLG)”, Ocean Eng. 38, 2011, 1555–1567.

Zhao W., Yang J., Hu Z. and Tao L., “Coupled analysis of nonlinear sloshing and ship motions”, Appl. Ocean Res. 47 (0), 2014a, 85–97.

Zhao W., Yang J., Hu Z., Xiao L. and Tao L., “Hydrodynamics of a 2D vessel including internal sloshing flows”, Ocean Eng. 84, 2014b, 45–53.

==Web references==

<div id="1"></div>
[[#cite-1|[1]]] CompassIS (last accessed: January 20th 2016): [http://www.compassis.com/seafem http://www.compassis.com/seafem].

<div id="2"></div>
[[#cite-2|[2]]] Centre Internacional de Metodes Numerics en Enginyeria (last accessed: January 20th 2016):  [http://www.cimne.upc.edu http://www.cimne.upc.edu].

<div id="3"></div>
[[#cite-3|[3]]] CEHINAV Res. Gr., DMFPA, ETSIN, Technical University of Madrid (last accessed: January 20th 2016). [http://canal.etsin.upm.es/ http://canal.etsin.upm.es/].

<div id="4"></div>
[[#cite-4|[4]]] AQUAgpusph (last accessed: January 20th 2016): [http://canal.etsin.upm.es/aquagpusph/ http://canal.etsin.upm.es/aquagpusph/].

<div id="5"></div>
[[#cite-5|[5]]] Tcl (Tool Command Language) [on-line] (last accessed: January 20th 2016): [http://www.tcl.tk http://www.tcl.tk]

<div id="6"></div>
[[#cite-6|[6]]] (last accessed: January 20th 2016): ([http://canal.etsin.upm.es/papers/servanetal2016/ http://canal.etsin.upm.es/papers/servanetal2016/])