An Analysis and Comparison of the Hydrodynamic Behavior of Ships Using Mesh-Based and Meshless Computational Fluid Dynamics Simulations

Abstract

This paper presents a comparison of two turbulence models implemented in two different frameworks (Eulerian and Lagrangian) in order to simulate the motion in calm water of a displacement hull. The hydrodynamic resistance is calculated using two open-source Computational Fluid Dynamics (CFD) software packages: OpenFOAM and DualSPHysics. These two packages are employed with two different numerical treatments to introduce turbulence closure effects. The methodology includes rigorous validation using a Wigley hull with experimental data taken from the literature. Then, the validated frameworks are applied to model a ship hull with a 30 m length overall (LOA), and their results discussed, outlining the advantages and disadvantages of the two turbulence treatments. In conclusion, the resistance calculated with OpenFOAM offers the best compactness of results and a shorter simulation time, whereas DualSPHysics can better capture the free-surface deformations, preserving similar accuracy.

1. Introduction

The study of the performance and dynamics of ship hulls is a subject of interest to all sectors in the naval world, from commercial ships to high-performance sailing ships, and even from displacement to planing hulls [1]. In the past, this type of analysis was exclusively performed by means of empirical models, simplified formulas, and tank tests. Nowadays, thanks to the development of computer technology, Computational Fluid Dynamics (CFD) approaches, based on the numerical resolution of viscous flow differential Navier–Stokes (NS) equations, have become largely available. CFD is becoming cheaper and faster and has provided increasingly accurate results, using highly complex numerical tools to study the performance and dynamics of ship hulls [2,3,4]. Fluid dynamics investigations can be carried out using the mesh-based CFD approach, which has the ability to calculate any complicated flow, but has the difficulty of deforming meshes (and also their possible distortion), or the particle method, which is based on the Lagrangian description of flow and is applied where the mesh-based method has difficulty [5].

CFD makes it possible to describe the turbulent motion that is generated around solid–fluid interfaces, such as hulls or aircraft [6]. Since exact solutions for turbulent flows are very complex and difficult to obtain computationally owing to heavy resolution requirements [7], turbulence closure models have been used to achieve practical calculation times [8]. Naturally, the model approximates the exact equations of motion. By their nature, turbulence models are primarily based on dimensional analysis combined with intuition and experimental observations [9].

Eulerian-based descriptions have been used to investigate the hydrodynamic performance of hulls. Wu et al. [10] studied a trimaran using CFD via Reynolds-averaged Navier–Stokes (RANS) equations. The results of a seakeeping experiment for the high-speed trimaran carried out at the China Ship Scientific Research Center (CSSRC) towing tank were also studied. Rather good agreement was shown between the computational and experimental results. Zha et al. [11] applied a RANS solver and a shear-stress transport (SST) k-$\omega$ turbulence model [12] to study hydrodynamic resistance at different Froude numbers ($Fr$) by comparing different types of hulls. Sulisetyono and Alifrananda [13] used RANS equations and an SST k-$\omega$ model to compare the experimental data of the DTMB 5415 hull with the NUMECA numerical code [14] and observed that the wake profile generated by the ship model compares well with the real case and the resistance results are very close to it. Desouky and Elhenawy [15] utilized a RANS solver and a k-$ϵ$ turbulence model to evaluate the free-surface flow around the KSC hull using the STAR-CCM+ (software version 14.02.010-R8) [16] and determined the hydrodynamics of the ship in the Suez Canal with good accuracy to understand the effects of squat in order to allow them to prevent ships in transit from having accidents. Islam and Guedes Soares [17] applied RANS equations and an SST k-$\omega$ model to estimate resistance in calm water and in the presence of waves for a ferry model using the open-source OpenFOAM toolkit. The paper showed that the solver is good for estimating resistance in calm water, but still has some shortcomings when it simulates waves with free motion. Furthermore, in another piece of work [18], for uncertainty in ship resistance prediction with OpenFOAM, both the factor of safety ($F_s$) and the correction factor ($C_i$) were estimated. The study was conducted on four different hull shapes, two container ship models (KCS, DTC), a crude oil carrier (KVLCC2), and a bulk carrier (JBC), and they were simulated in calm sea conditions; from the results the total drag coefficient, sinking, and trim were predicted. Kim et al. [19] studied the KSC hull model via STAR CCM+. The resistance performance and flow structure around the ship were analyzed. This study confirms and shows that full-scale results can be obtained at model scale by applying virtual fluid instead of full-scale numerical simulations, which require more computational resources.

Likewise, Lagrangian-based numerical approaches have been used to study hull dynamics. Kawamura et al. [20] validated the accuracy of the Smoothed Hydrodynamic Particle (SPH) method by comparing it with free and constrained motion tests in storm conditions, using a simplified ship model. The results showed that the SPH simulation has good potential for the quantitative assessment of vessel safety in seagoing situations. Priyambada and Tarwidi [21] applied artificial viscosity to study the vertical displacement of a floating object on free-surface flow, which was evaluated by varying the frequency of the piston wavemaker, and Eriksson [22] chose the KSC hull as a case study and evaluated the change in total drag coefficient using the SPH and laminar-and-sub-particle scale (SPS) models. Cheng et al. [23] used a weakly compressible SPH implementation to investigate the water entry and the ensuing slamming loads on a ship cross-section, and then extended the validity of their analysis to a three-dimensional (3D) framework. The emphasis of this work was on the stability and robustness of their implementation, which substantiated into capturing the spray features and the hull pressure field. Tagliafierro et al. [24] studied the spray phenomena around a planing hull by applying DualSPHysics code [25]. They calculated the hydrodynamic drag and running attitudes for a wide $Fr$ range by comparing them with the experimental results and reproduced the spray around the planing hull. Capasso et al. (2023) [26,27] studied the modeling of free-surface flows interacting with a planing hull under regular wave conditions using DualSPHysics. Overall, the results are in good agreement with the experiments.

Throughout the reviewed literature, mesh-based and meshless methods have been shown to possess strengths and weaknesses for the hydrodynamic simulation of hulls in calm or rough water conditions [6,28,29]. With respect to the usability of each approach, code-to-code comparisons become very important to optimize their usage for engineering applications. In this regard, Gruwez et al. [30] compared the capability of OpenFOAM and DualSPHysics in reproducing a large-scale wave channel experiment of bichromatic wave transformations on a steeply sloping, very shallow shoreline, showing similar outcomes. The results showed that overall, OpenFOAM provides the highest model skill, but has the highest computational cost. DualSPHysics is shown to have reduced model performance, but still comparable to OpenFOAM and at a lower computational cost. Meringolo et al. [31] performed a first comparison between the SPH method and Finite Volume Method (FVM) for multiphase flows. The numerical results were compared with those obtained experimentally, showing good agreement. The results showed very good agreement in terms of potential and kinetic energy component evolution among the proposed model, the Large Eddy Simulation (LES)–FVM solver, and the single-phase $\delta$-SPH model.

By recognizing the high value of turbulence models for engineering applications, in this work, two open-source tools, OpenFOAM (v2012) (mesh-based) and DualSPHysics (v5.2) (mesh-free), largely used for most applications reported here, are used to study the fluid flow velocity, vorticity, and forces acting on a ship surface in calm water conditions.

The main objective of this work is to compare the two open-source models in studying ship hydrodynamics in terms of total resistance, free-surface prediction, and wake evolution, as well as perform an analysis of the simulation runtimes. In order to validate the adopted numerical models, the experimental results obtained for Wigley hull [32,33] were used for comparison. Validation was carried out at three different $Fr$ values: 0.25, 0.316, and 0.408, respectively. Furthermore, the validated frameworks were applied to model newly designed ship hulls with a 30 m length overall (LOA), and the results obtained with the SPH and FVM approaches at different $Fr$ values were compared and analyzed.

The rest of this paper is organized as follows. In Section 2, the numerical model and numerical methods are introduced; in Section 3, the numerical setup is described; in Section 4, the case studies are illustrated; in Section 5, the results are presented and discussed; Section 6 provides the concluding remarks.

2. Numerical Models

This section describes the numerical models used by OpenFOAM and DualSPHysics for our case study.

2.1. Eulerian Approach in OpenFOAM

We consider the Navier–Stokes equations in a domain characterized by the presence of different fluids. The continuity and momentum equations are written for the Eulerian solver, considering the incompressible fluid approximation. This system of equations can be expressed as

$${\displaystyle \frac{\partial u_i}{\partial x_i}}=0$$

(1)

$${\displaystyle \frac{\partial (\rho u_i)}{\partial t}}+u_j{\displaystyle \frac{\partial (\rho u_i)}{\partial x_j}}=-{\displaystyle \frac{\partial p}{\partial x_i}}+\mu {\displaystyle \frac{{\partial }^2{u}_i}{\partial x_i{x}_j}}+g_i$$

(2)

where $u_i$ (i = 1, 2, 3) denotes flow velocity on the x-, y-, and z-axes, respectively, $\rho$ represents the mixture density of the two phases separately considered, p denotes pressure, $\mu$ is the dynamic viscosity of the fluid, and $g_i$ represents the components of the gravitational acceleration.

The Reynolds averaging approach was adopted to eliminate the fluctuating components of the solution variables in the instantaneous NS equations to reduce the computational cost. The average of the continuity and NS equations yields

$${\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_i}}=0$$

(3)

$${\displaystyle \frac{\partial {(\rho \overline{u}}_i)}{\partial t}}+{\overline{u}}_j{\displaystyle \frac{\partial {(\rho \overline{u}}_i)}{\partial x_j}}=-{\displaystyle \frac{\partial \overline{p}}{\partial x_i}}+\mu {\displaystyle \frac{{\partial }^2{\overline{u}}_i}{\partial x_i{x}_j}}-{\displaystyle \frac{\partial \left(\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)}{\partial x_j}}+g_i$$

(4)

where ${\overline{u}}_i$ and $\overline{p}$ represent the ensemble-averaged velocity and pressure, respectively; the product $\rho \overline{u_i^{\prime }{u}_j^{\prime }}$ denotes the Reynolds stress component. Thereafter, the RANS equations govern the transport of flow. The Reynolds stress in Equation (4) is modeled as a function of turbulent viscosity ${\mu }_T$ and the mean velocity gradients according to the Boussinesq hypothesis:

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_T\left({\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_j}}+{\displaystyle \frac{\partial {\overline{u}}_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\rho k{\delta }_{ij},$$

(5)

where k is the turbulent kinetic energy (TKE) and ${\delta }_{ij}$ is the Kronecker delta. The latter is a function of two discrete variables, in particular, of two variables on the integers or naturals, which is worth 1 if their values coincide, and worth 0 otherwise.

To obtain the turbulent viscosity, we use the SST k-$\omega$ turbulence model, which is a two-equation eddy-viscosity model. It is a hybrid model combining the Wilcox k-$\omega$ and the k-$ϵ$ models. The use of a k-$\omega$ formulation in the inner parts of the boundary layer makes the model directly usable all the way down to the wall through the viscous sub-layer; hence, the SST k-$\omega$ model can be used as a Low-Re turbulence model without any extra damping functions. The SST formulation also switches to a k-$ϵ$ behavior in the free stream and thereby avoids the common k-$\omega$ problem whereby the model is too sensitive to the inlet free-stream turbulence properties. A blending function, $F_1$, activates the Wilcox model near the wall and the k-$ϵ$ model in the free stream. The model is proposed by Menter [34]. In the SST k-$\omega$ model, the viscosity term ${\mu }_T$ is defined as:

$${\mu }_T=\rho {\displaystyle \frac{k}{\omega }}$$

(6)

while the TKE per unit mass k and the dissipation rate $\omega$ obey the following equations:

$$\rho {\displaystyle \frac{\partial k}{\partial t}}+\rho {\overline{u}}_j{\displaystyle \frac{\partial k}{\partial x_j}}={\tau }_{ij}{\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_j}}-{\beta }^{\ast }\rho k\omega +{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\sigma }_k{\mu }_T\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]$$

(7)

$$\rho {\displaystyle \frac{\partial \omega }{\partial t}}+\rho {\overline{u}}_j{\displaystyle \frac{\partial \omega }{\partial x_j}}=\alpha {\displaystyle \frac{\omega }{k}}{\tau }_{ij}{\displaystyle \frac{\partial {\overline{u}}_i}{\partial x_j}}-\beta \rho {\omega }^2+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\sigma }_{\omega }{\mu }_T\right){\displaystyle \frac{\partial \omega }{\partial x_j}}\right]$$

(8)

In the above equations, $\alpha =0.52;{\sigma }_k=0.5;{\sigma }_{\omega }=0.5,$ whereas $\beta$ and ${\beta }^{\ast }$ are computed through auxiliary functions, as presented in Wilcox [35]. The second selected turbulence model, the SST k-$\omega$ [12], has several relatively minor variations with respect to the original SST version [34,36]. The TKE and the dissipation rate are computed using the following equations:

$$\rho {\displaystyle \frac{\partial k}{\partial t}}+\rho {\overline{u}}_j{\displaystyle \frac{\partial k}{\partial x_j}}={\tilde{P}}_k-{\beta }^{\ast }\rho k\omega +{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\sigma }_{k1}{\mu }_T\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]$$

(9)

$$\rho {\displaystyle \frac{\partial \left(\omega \right)}{\partial t}}+\rho {\overline{u}}_j{\displaystyle \frac{\partial \left(\omega \right)}{\partial x_i}}=\alpha \rho S^2-\beta \rho {\omega }^2+{\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\sigma }_{\omega }{\mu }_T\right){\displaystyle \frac{\partial \omega }{\partial x_j}}\right]+2\left(1-F_1\right)\rho {\sigma }_{\omega ,2}{\displaystyle \frac{1}{\omega }}{\displaystyle \frac{\partial k}{\partial x_j}}{\displaystyle \frac{\partial \omega }{\partial x_j}}$$

(10)

where ${\tilde{P}}_k$ represents a production limiter used in the model to prevent the build-up of turbulence in stagnation regions [12], and $F_1$ represents the blending function [12,34,36,37]. All the constants are predicted through a blend from the corresponding constants [e.g., $\alpha ={\alpha }_1{F}_1+{\alpha }_2(1-F_1)$]. This model contains the following closure coefficients: ${\beta }^{\ast }=9/100$, ${\sigma }_{k1}=0.85$, ${\sigma }_{\omega 1}=0.5$, ${\sigma }_{\omega 2}=0.856$, ${\alpha }_1=5/9$, ${\alpha }_2=0.44$, and ${\beta }_1=3/40$, as suggested in [37]. The choice of the SST k-$\omega$ model is due to its generally superior performance compared to the classical two-equation model in a variety of complex flow cases [12,37,38].

The volume of fluid (VOF) method with the artificial bounded compression technique proposed by Hirt and Nichols (1981) [39] was used to capture the free surface. The free surface was considered to be a mixture of water and air, and, accordingly, the VOF transport equation was expressed as

$${\displaystyle \frac{\partial \alpha }{\partial t}}+\nabla ·{\mathit{U}}_{\mathit{r}}(1-\alpha )\alpha =0,$$

(11)

where ${\mathit{U}}_{\mathit{r}}={\mathit{U}}_{\mathit{water}}-{\mathit{U}}_{\mathit{air}}$ is the relative velocity between the water and air, and is termed “compression velocity”; $\alpha$ represents the volume fraction, defined as the relative volume proportion of water in a cell; and ∇ represents the nabla operator. If the volume is completely filled with air, $\alpha$ is 0, but if the volume is filled with water, $\alpha$ is 1; $\alpha$ between 0 and 1 indicates the multiphase fluid interface. In our case, $\alpha$ equal to 0.5 indicates the free water surface.

The spatial variation in fluid density and dynamic viscosity can be then expressed as

$$\rho =\alpha {\rho }_{water}+\left(1-\alpha \right){\rho }_{air}$$

(12)

$$\mu =\alpha {\mu }_{water}+\left(1-\alpha \right){\mu }_{air}$$

(13)

The sequential solution implies that the velocity components are solved independently and only coupled through explicit terms. It follows a stability criterion that limits the duration of the time step, known as the Courant–Friedrichs–Lewy (CFL) condition [40]. The CFL number is given by

$$CFL={\displaystyle \frac{\left|u\right|\Delta t}{\Delta x}}$$

(14)

2.2. Lagrangian Approach in DualSPHysics

The DualSPHysics model uses the meshless SPH method to describe fluid mechanics and is used to discretize a fluid volume as a collection of particles, whose motion is dictated by the NS equations. Particles represent nodal points at which physical quantities (e.g., position, velocity, density, pressure) are approximated by interpolating the values of the NS equations. The interpolation of the values of neighboring particles takes place on a compact support at a short distance. The technique is ideal for studying violent flows owing to its inherent lack of mesh distortion and can easily handle multi-phase simulations because each particle stores its own properties. The SPH method has been used to describe a variety of free-surface flows (wave propagation on beaches, breakers, impact on structures and dams) [25].

The mathematical foundations of the SPH method are based on the approximation of any quantity by convolution integrals. Any function F can be defined by

$$F\left(\mathit{r}\right)=\int F\left({\mathit{r}}^{\prime }\right)W\left(\mathit{r}-{\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }$$

(15)

where W is the kernel function [41], $\mathit{r}$ is the position of the point at which the function is calculated, and ${\mathit{r}}^{\prime }$ is the position of a generic calculation point. This function F can be approximated by the contributions of the interpolating particles; a summation is performed on all particles within the compact kernel support:

$$F\left({\mathit{r}}_a\right)\approx \sum _{b}F\left({\mathit{r}}_b\right)W\left({\mathit{r}}_a-{\mathit{r}}_b,h\right){\displaystyle \frac{m_b}{{\rho }_b}}$$

(16)

where a is the interpolated particle, b is a neighboring particle, m and $\rho$ are the mass and density, respectively, ${\displaystyle \frac{m_b}{{\rho }_b}}$ is the volume associated with the neighboring particle, and h is the smoothing length. The kernel functions W satisfy several properties, such as positivity on compact support, normalization, and a monotonic decrease with distance [42]. Several methodologies are available; one option is the Quintic Wendland polynomial kernel [43]:

$$W\left(q\right)={\alpha }_D{\left(1-{\displaystyle \frac{q}{2}}\right)}^4\left(2q+1\right)0\le q\le 2$$

(17)

where ${\alpha }_D$ is a real number such that the kernel ensures the normalization property, $q={\displaystyle \frac{r}{h}}$ is the non-dimensional distance between the particles, and r is the distance between a certain particle a and another b. In this way, the Wendland kernel is used to calculate the interactions between the particles at a distance of up to 2h. In the Lagrangian framework, the differential form of the NS equations can be written as a discrete version using the kernel function:

$${\displaystyle \frac{d{\mathit{v}}_{\mathit{a}}}{dt}}=-\sum _b{m}_b{\displaystyle \frac{P_a+P_b}{\left({\rho }_a{\rho }_b\right)}}{\nabla }_a{W}_{ab}+{\mathbf{\Gamma }}_a+\mathit{g}$$

(18)

$${\displaystyle \frac{d{\rho }_a}{dt}}={\rho }_a\sum _b{\displaystyle \frac{m_b}{{\rho }_b}}{\mathit{v}}_{ab}{\nabla }_a{W}_{ab}+2\delta hc\sum _b({\rho }_b-{\rho }_a){\displaystyle \frac{{\mathit{v}}_{ab}{\nabla }_a{W}_{ab}}{\parallel {\mathbf{r}}_{ab}{\parallel }^2}}{\displaystyle \frac{m_b}{{\rho }_b}}$$

(19)

where t is the time, $\mathit{v}$ is the velocity, P pressure, $\mathit{g}$ is the gravitational acceleration, ${\nabla }_a$ is the gradient operator, $W_{ab}$ the kernel function, whose value depends on the distance between a and b, ${\mathit{r}}_{\mathit{ab}}$ = ${\mathit{r}}_a$−${\mathit{r}}_b$ with ${\mathit{r}}_k$ being the position of the kth particle, and c is the speed of sound. The code allows the use of an artificial viscosity model (not employed here) or a laminar viscosity model + SPS turbulence model [42]. An SPS model in combination with a laminar viscosity model resolves small unresolved motions. Vortices that can be solved by the grid are allowed to evolve according to the NS equations, and a model is used to represent turbulence at scales below the grid. The term ${\mathbf{\Gamma }}_a$ in Equation (18) introduces dissipation, as described by

$${\mathbf{\Gamma }}_a=\sum _b{m}_b{\displaystyle \frac{4{\nu }_0{\mathit{r}}_{ab}{\nabla }_a{W}_{ab}}{({\rho }_a+{\rho }_b)(\parallel {\mathbf{r}}_{ab}{\parallel }^2+0.01h^2)}}{v}_{ab}+\sum _b{m}_b\left({\displaystyle \frac{{\overrightarrow{\tau }}_a^{ij}+{\overrightarrow{\tau }}_b^{ij}}{{\rho }_a{\rho }_b}}\right){\nabla }^i{W}_{ab}$$

(20)

The first term in Equation (20) introduces the so-called laminar viscosity [44], where the term ${\nu }_0$ = $1.0\times {10}^{-6}{\mathrm{m}}^2{\mathrm{s}}^{-1}$ refers to the fluid kinematic viscosity. The second term describes the SPS model [45], corresponding to the variationally consistent form of the symmetric formulation proposed in [44]. The SPS stress tensor, ${\overrightarrow{\tau }}^{ij}$, in Einstein notation in the coordinate directions i and j, is defined according to the SPS strain tensor:

$${\overrightarrow{S}}^{ij}=-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\partial {\overrightarrow{v}}^i}{\partial x^j}}+{\displaystyle \frac{\partial {\overrightarrow{v}}^j}{\partial x^i}}\right)$$

(21)

where x is the position of the particle. A relationship between density and pressure binds the system of equations. DualSPHysics uses a weakly compressible SPH (WCSPH) formulation to model Newtonian fluids, and, for this formulation, the Tait equation of state is used to determine the fluid pressure, P, from the density of the particles.

$$P={\displaystyle \frac{c_0^2{\rho }_0}{\gamma }}\left({\left({\displaystyle \frac{\rho }{{\rho }_0}}\right)}^{\gamma }-1\right)$$

(22)

where ${\rho }_0$ is the density of the reference fluid, and $\gamma$ is the polytropic constant. The compressibility of the fluid is adjusted so that $c_0$ can be artificially lowered to ensure reasonable time values.

The explicit second-order accurate Symplectic [46] time integration scheme takes the following form, considering the weakly compressible formulation, for solving the position of particle a at step n + 1:

$${\mathit{r}}_a^{n+{\displaystyle \frac{1}{2}}}={\mathit{r}}_a^n+{\displaystyle \frac{\Delta t}{2}}{\mathit{v}}_a^n$$

(23)

$${\mathit{v}}_a^{n+{\displaystyle \frac{1}{2}}}={\mathit{v}}_a^n+{\displaystyle \frac{\Delta t}{2}}{\mathit{F}}_a^n$$

(24)

$${\mathit{v}}_a^{n+1}={\mathit{v}}_a^n+{\displaystyle \frac{\Delta t}{2}}{\mathit{F}}_a^{n+{\displaystyle \frac{1}{2}}}$$

(25)

$${\mathit{r}}_a^{n+1}={\mathit{r}}_a^n+\Delta t{\displaystyle \frac{{\mathit{v}}_a^{n+\frac{1}{2}}+{\mathit{v}}_a^n}{2}}$$

(26)

where ${\mathit{F}}_a={\displaystyle \frac{d{\mathit{v}}_a}{dt}}$, and ${\mathit{v}}_a={\displaystyle \frac{d{\mathit{r}}_a}{dt}}$. The density is updated at the next time step by means of a similar two-step strategy, which reads

$${\rho }_a^{n+{\displaystyle \frac{1}{2}}}={\rho }_a^n+{\displaystyle \frac{\Delta t}{2}}{R}_a^n$$

(27)

$${\rho }_a^{n+1}={\rho }_a^n+{\displaystyle \frac{2+\frac{R_a^{n+\frac{1}{2}}}{{\rho }_a^{n+\frac{1}{2}}}\Delta t}{2-\frac{R_a^{n+\frac{1}{2}}}{{\rho }_a^{n+\frac{1}{2}}}\Delta t}}$$

(28)

where $R_a={\displaystyle \frac{d{\rho }_a}{dt}}$.

As illustrated in Monaghan et al. (1999) [47], variable time stepping can be used and, considering the CFL condition, the force terms and the viscous diffusion term can be calculated as follows:

$$\Delta t=C_{CFL}\underset{a}{min}(\sqrt{{\displaystyle \frac{h}{|f_a|}},{\displaystyle \frac{h}{c+\underset{b}{max}\frac{|h{\mathit{v}}_a·{\mathit{r}}_a|}{\parallel {\mathbf{r}}_{ab}{\parallel }^2+0.01h^2}}}})$$

(29)

where $f_a$ is the force per unit mass.

For this work, the modified Dynamic Boundary Condition (mDBC) [48] is used. Its predecessor approach—Dynamic Boundary Condition (DBC) [49]—has been shown to produce some non-physical gaps when free-surface flows interact with solid particles. In [50], the authors show that this issue can be addressed through the use of “Dummy” particles, following the implementation proposed in [51]. Instead, mDBC, by sharing the same particle arrangement as DBC, locates the interacting boundary surface at the middle line between the outermost layer of body particles and the fluid domain. Normal vectors are calculated with respect to this latter boundary surface. The solid body adjusts with a number of boundary particle layers to ensure the completeness of the kernel of fluid particles approaching the boundary, avoiding inconsistencies due to truncation effects. Once the geometry has been characterized by the normal vectors, the position of the interface is used to reflect the ghost nodes in the fluid domain.

An isotropic particle spacing is an important stability problem in SPH, especially in violent flows, as the particles cannot maintain a uniform distribution. It is corrected by the shifting algorithm [52]. Finally, it is important to mention that the initial condition in DualSPHysics is generated using a pre-processing tool that creates particles with an initial inter-particle distance $D_p$. This value also defines the resolution used in the simulations. Fluid and solid particles, whether part of a floating object or other solids, are then created following this initial spacing at the initial time step.

3. Numerical Setup

This section shows how to obtain the setups of the two working models. For numerical feasibility, the motion of the hull is simulated by a free-surface flow whose velocity corresponds to the speed of the hull, and all of the motions of the hull are completely locked. The Wigley hull is a well-known case for which measurements have been made from several authors [32,33]. The Wigley hull form is conceived to provide data for both physics and CFD validation for a simplified ship hull. Three different $Fr$ values of 0.250, 0.316, and 0.408 are taken from the tests in order to compare the resistances with the two frameworks used. It is typical for the Wigley hull to have a slender shape (Figure 1a), which is generally accomplished by having a relationships between LOA, Breadth (B), and Draft (D) of L/B = 10 and B/D = 1.6 (which are the most common values observed in our literature research). The Wigley hull has mainly academic interest owing to its simplicity. Though it is rarely used at real scales, it is commonly found at different model sizes, as well as tested with several $Fr$ values. It is used regularly for tests and validation work for new CFD software, as well as for theoretical studies. Its basic design is also an advantage as it enables it to perform hull optimization studies: there are few parameters to control, and the modification of the ship surface can be performed without many restrictions or smoothness problems. In simulations, with both approaches, the domain must be created according to the following: Upstream Channel Length (UCL), Downstream Channel Length (DCL), Channel Width (CW), and Channel Depth (CD) (Figure 1b).

3.1. Validation of the Numerical Model in OpenFOAM

In the CFD using OpenFOAM, the steps of the simulation include designing the computation domain (

Table 1. Computational domain lengths.

UCL = 2.5 LOA = 2.5 m

DCL = 4.5 LOA = 4.5 m

CW = 3 LOA = 3 m

CD = 48 D = 3 m

), generating a mesh, and setting up the computed conditions. The most important steps for Eulerian models are those required to create the mesh and to generate the mesh; OpenFOAM uses the routines blockMesh, topoSet, and snappyHexMesh. The snappyHexMesh generate an unstructured mesh including hexahedral and tetrahedral cells.

A mesh sensitivity study is conducted consisting of running the same simulation using grids with different resolutions and analyzing how much the converged solution changes with each mesh. Three different simulations are performed before selecting the final configuration of the grid. The features of the tested grids are outlined in

Table 2. The grid refinement progression.

Cells	First Mesh	Second Mesh	Third Mesh
Min element size	0.69 mm	0.62 mm	0.57 mm
Max element size	75.0 mm	50.0 mm	37.5 mm
Count	3,091,284	4,464,400	6,180,196
Runtime	2.15 h	3.22 h	4.62 h

. The mesh element counts and runtime for each test are reported too. The simulations are run on a Central Processing Unit (CPU) with an Intel® Xeon® Silver 4116 Core: 48 processor.

The final three-dimensional (3D) finite-volume computational domain includes about 6.2 × 10⁶ grid points, with a minimum grid spacing of 0.57 mm and a maximum grid spacing of 37.5 mm. Figure 2 visualizes the mesh in three directions X [a], Y [b], and Z [c], and meshing around the hull [d].

The physical domain consists of water and air under specific and homogeneous conditions. The two-phase model uses the VOF phase-fraction interface capturing approach, in which the free surface is well defined over the entire domain. The interFoam solver (embedded in the OpenFOAM^® C++ libraries) is used to numerically solve governing Equations (3) and (4), together with the equations of the turbulence model. In the case of OpenFOAM, the simulation time step varies between 5 × 10⁻² and 1 × 10⁻⁵ s.

The independent convergence study is shown in Figure 3, where the total resistance of the Wigley hull for the three different mesh sizes is considered, and for the case with $Fr$ = 0.408. Considering both accuracy and computational efficiency, a minimum mesh size of 0.57 mm and a maximum mesh size of 37.5 mm are considered to be adequate for the current application, resulting in about 6.2 million mesh cells in the computational domain.

The total resistance coefficients of a ship in calm water, $C_T$, are defined by the following equation:

$$C_T={\displaystyle \frac{R_T}{0.5\rho S{V}^2}}$$

(30)

where $C_T$ is the total resistance coefficient; $R_T$ is the total resistance acting on the hull, N; S is the wetted surface area, m²; $\rho$ is the water density, kg m⁻³; and V is the velocity of the ship, m s⁻¹.

Table 3. A comparison of the resistance coefficients.

$\mathit{Fr}$	${\mathit{R}}_{\mathit{T}}$ [N]	${\mathit{C}}_{\mathit{T}}$ × 10−3	${\mathit{C}}_{\mathit{exp}}$ × 10−3	Error [%]
0.250	0.301	5.26	4.55	−15.71
0.316	0.488	5.34	5.16	−3.48
0.408	0.843	5.53	5.73	3.49

shows that the numerical values of $C_T$ are in agreement with respect to the experimental results for all $Fr$ values. The maximum relative error is −15.71% at the lowest $Fr$. For $Fr$ = 0.25 and 0.316, the value of the resistance coefficient is overestimated, whereas for $Fr$ = 0.408, it is underestimated.

3.2. Validation of the Numerical Model in DualSPHysics

Initially, the domain is dimensioned according to the requirements of the wake development (

Table 4. Computational domain lengths.

UCL = 0.8 LOA = 0.8 m

DCL = 1.6 LOA = 1.6 m

CW = 1 LOA = 1.0 m

CD = 6 D = 0.375 m

). Therefore, UCL and DCL are dimensioned to allow the wake to develop undisturbed. To avoid obstruction effects by side walls, the width of the domain was appropriately designed. For the definition of the CD, an adequate value is defined to avoid any shallow-water effect.

DualSPHysics allows boundary conditions to be included via buffer zones, which are then characterized as inflow/outflow areas [53]. These zones are made up of particles to guarantee the compactness of the kernel, thus ensuring a safe process of particle addition/removal. In these areas, physical quantities such as velocity and density can be extrapolated via ghost nodes in the neighboring fluid region, where fluid quantities are calculated by interpolation and applying the procedure proposed by Liu and Liu (2006) [54]. The free-surface flow is simulated with inlet and outlet contours. The side walls, which are parallel to the flow direction (y), and the bottom are considered fixed, as if the fluid were contained in a channel. In DualSPHysics, an important step is particle sizing, which is useful for good case reproduction. Three different initial particle sizes are considered, $D_p$ equal to 0.009, 0.007, and 0.005 m, and based on these, simulations are run taking into consideration $Fr$ = 0.408 (

Table 5. Model resolution and runtimes.

${\mathit{D}}_{\mathit{p}}$ [m]	Particle Number	Runtime [h]
0.009	1,576,534	1.33
0.007	3,149,391	3.20
0.005	8,198,736	13.42

). Figure 4 shows the total resistance with different particle sizes for the case with $Fr$ = 0.408. The simulations are run on a Graphics Processing Unit (GPU) and not on a CPU. The GPU allows a reduction in times compared to the CPU and is an NVIDIA GeForce RTX 2080. In the case of DualSPHysics, the simulation time step varies between 1 × 10⁻⁴ and 1 × 10⁻⁵ s.

The second particle size (0.007 m) is chosen because it has good convergence and is more time-saving than smaller diameters. Finally, the resistances and resistance coefficients for the other Froude numbers are calculated.

Table 6. A comparison of the resistance coefficients.

Fr	${\mathit{R}}_{\mathit{T}}$ [N]	${\mathit{C}}_{\mathit{T}}$ × 10−3	${\mathit{C}}_{\mathit{exp}}$ × 10−3	Error [%]
0.250	0.284	4.95	4.55	−8.87
0.316	0.493	5.39	5.16	−4.53
0.408	0.908	5.95	5.73	−3.87

shows that the numerical values of $C_T$ are in agreement with respect to the experimental results for all $Fr$ values. The maximum relative error is −8.87% at the lowest $Fr$. The value of the resistance coefficient is overestimated.

From the analysis of data reported in Table 3 and Table 6, it is evident how the two models are able to correctly evaluate the values of the resistance coefficients. For both models and for all the considered Froude numbers, the value of the resistance coefficient is underestimated. The Lagrangian model, in this case, demonstrates better overall behavior than the Eulerian model in calculating the resistance coefficients, with smaller errors.

4. Case Study

The validated frameworks are applied to model a ship hull with a 30 m LOA. The general view and the main characteristics of the new designed hull are reported in Figure 5 and

Table 7. Main geometric characteristics of ship.

LOA	30.00 m
B	6.00 m
D	2.00 m
S	168.96 m2

, respectively.

The two models have two different setups, as can be seen from the difference in the methods they use.

4.1. Simulations via OpenFOAM

For the construction of the case, the same design choices are adopted as in the validation for the Wigley case. The computational domain, inserted in the blockMeshdict, is proportional to that created in the validation. The number of cells is 6,245,338. The cell size ranges from a minimum size of 0.0018 m to a maximum size of 1.2 m. Also, here, the turbulence model is the SST k-$\omega$ model. It is solved by the routine interFoam.

4.2. Simulations via DualSPHysics

For the construction of the case, the same design choices are adopted as in the validation for the Wigley case, and we also use discretization proportional to the validation case. The particle size has $D_p$ = 0.20 m. The number of particles generated is 7,580,338. This value is sufficient to accurately discretize the external characteristics of the hull, as the surface is smooth and free of sudden variations. After the static convergence, the domain is dimensioned according to the requirements of the wake development (

Table 8. Domain simulation.

UCL	22.5 m
DCL	67.5 m
CW	50 m
CD	12 m

).

The boundary conditions and parameters are also the same as in the previous case. Also, here, the turbulence model is the SPS model.

5. Results and Discussion

This section summarizes the results obtained from the newly designed ship hull by elaborating data derived from the numerical simulations. Three $Fr$ values are considered, 0.233, 0.291, and 0.408, which correspond to speeds of 4.0, 5.0, and 7.0 m/s, respectively, highlighting the resistance of the two frameworks.

Table 9. Resistance coefficient comparison.

Framework	$\mathit{Fr}$	${\mathit{R}}_{\mathit{T}}$ [kN]	${\mathit{C}}_{\mathit{T}}\times {10}^{-3}$	Runtime [h]
OpenFOAM	0.233	26.573	19.683	2.08
OpenFOAM	0.291	51.715	24.515	2.1
OpenFOAM	0.408	114.338	27.654	2.22
DualSPHysics	0.233	29.622	21.941	16.8
DualSPHysics	0.291	49.966	23.686	14.62
DualSPHysics	0.408	127.649	30.873	12.45

shows the resistance coefficients obtained from OpenFOAM and DualSPHysics and at different values of $Fr$. They give similar results. The forces obtained from the Eulerian approach are slightly lower than those obtained from the Lagrangian approach except for $Fr$ = 0.291. However, the differences in values between the two models show an absolute error rate of no more than 11%. This leads to an awareness that the two models used show very similar force and drag coefficient results despite the difference in the basic application of the two methods.

Furthermore, Table 9 shows the execution time for the two models, suggesting very big computational time differences between the models. The simulations with OpenFOAM have runtimes that are, on average, an order of magnitude shorter than those with DualSPHysics, in addition to the fact that the two numerical methods are fundamentally different. Another difference between the two methods is the hardware used. OpenFOAM simulations are run on CPUs, which are very powerful hardware for these simulations. DualSPHysics simulations are run on GPUs, and it is well known that the SPH method has a high computational cost. Simulation time is an important factor in hull optimization. In addition to optimization time, a significant reduction in CFD time improves the timing of obtaining results.

Figure 6 shows the total resistance obtained in the two approaches. The box size indicates the range of variation in resistance values, and the line in the box indicates the average value. The compactness of values is defined by the difference in the two extremes (whiskers) extending 1.5 times from the interquartile range from the top and bottom of the box. In addition, all outliers (which are the values that go outside the range) are eliminated.

OpenFOAM shows greater compactness of the results in the simulation. Only the $Fr$ of 0.408 shows a small difference between the values, which could already be seen from the value of the resistances.

Flow Visualization

In Figure 7, Figure 8 and Figure 9, there are six different panels representing the free-surface elevation in two different views (OpenFOAM and DualSPHysics) for three different simulation instants (0.25, 0.5, and the end of the simulation time, T), all of them related to the simulation with $Fr$ = 0.291. In Figure 7, Figure 8 and Figure 9, on the left (a,c,e), there are those generated by OpenFOAM and, on the right (b,d,f), those generated by DualSPHysics.

Figure 7 shows the free-surface elevation scaled with the LOA ship length. It can be seen that quantitative and qualitative agreement is reached around the ship at the time of the final simulation for both wave position and amplitude. In the first quarter of the simulation, the two simulations are similar in the bow and stern areas. In the middle of the simulations, the free surface reproduced in OpenFOAM shows more uplift. At the end of the two simulations, the two models show a good representation of the free surface with small differences. OpenFOAM simulates the free surface around the hull well, but does not propagate it beyond the end of the first refinement box owing to cell variation. DualSPHysics, on the other hand, shows better free-surface quality and excellent propagation, although it has a small reflection at the end of the domain generated by the side wall.

Figure 8 shows the velocity field at three different simulation times. The velocity field expressed by the OpenFOAM numerical model has a large area of velocity reduction around the bow and a narrow area of velocity reduction in the wake, in contrast to the velocity field generated by the DualSPHysics model, which has a much smaller area of velocity reduction around the bow and a very large area of velocity reduction in the wake. At one quarter of the simulation, the two simulations are nearly identical. At half of the simulations, the velocity field generated by OpenFOAM shows a larger area of velocity reduction around the bow than the velocity field developed in DualSPHysics. Finally, at the last instant of the simulation, similarity of the velocity field can be seen, but with a greater velocity reduction in the bow area by the simulation in OpenFOAM, and in the stern area, a narrower area of velocity reduction than in DualSPHysics. Also, from these figures, one can see the defect generated by cell variation in the OpenFOAM simulations and a small reflection in the DualSPHysics simulation.

Figure 9 shows the vorticity field at three different simulation times. During the test, both models accurately resolved the vorticity field, as evidenced by the similar vorticity fields displayed. The images represent the section along the y-direction in the middle of the domain. In the first quarter of the simulation, the two simulations are similar in the aft area, with the difference that in the DualSPHysics simulation, some particles are lost. The particles decrease owing to the anisotropic spacing of the particles, which is a stability problem of the model. In the middle of the simulations, the vorticity reproduced behind the stern in OpenFOAM shows a larger area. At the end of the two simulations, the two models show a good representation of vorticity, with minimal differences. From the vorticity field section, we can also see the difference in the accuracy of the free surface and how the simulation in DualSPHysics captures the free surface better.

6. Conclusions

In this paper, two turbulence methods are presented for calculating the total resistances acting on the ship under study. The turbulence models selected for analysis are the SST k-$\omega$ model, which provides high-precision boundary layer modeling as a combination of k-$ϵ$ and k-$\omega$, and the Laminar-and-SPS model, which is mainly used to adequately present the viscosity and turbulent motions of the fluid flow, especially in the case of the evolution of breaking waves, but adapted to our case study. In this work, two different numerical schemes, one mesh-based and one meshless, each one employing a different turbulence method are presented for calculating the total resistances acting on the ship under study. For the mesh-based approach, the OpenFOAM model is employed. For the meshless approach, the Weaky Compressible SPH-based DualSPHysics model is used. To predict the resistance of ships with different running attitudes, numerical simulations of ships advancing at different speeds and $Fr$ values were carried out.

A comparison of the tools may be useful for possible use in the geometry optimization field and, as a consequence, considering that the naval world is increasingly approaching an electric propulsion system, good hull optimization and knowledge of fluid dynamics would improve the range of electric ships. In the present work, the accuracy of the calculations is ensured by preventive validation of the two methods described in Section 3. The numerical results from the two methods show good agreement compared to the experimental references for all $Fr$ values, with a maximum relative error of −15.71% at the lowest $Fr$ for the SST k-$\omega$ turbulence model. The drag coefficients obtained by OpenFOAM and DualSPHysics and studied at the different values of $Fr$ show very similar results, with an absolute error rate of no more than 11%. The free surface in the two simulations shows quantitative and qualitative agreement for both wave position and amplitude. Regarding the velocity field, the one expressed by the OpenFOAM numerical model has a large area of velocity reduction around the bow and a narrow area of velocity reduction in the wake, in contrast to the velocity field generated by the DualSPHysics model, which has a much smaller area of velocity reduction around the bow and a very large area of velocity reduction in the wake. As evidenced by the vorticity fields shown, both models accurately resolved the vorticity field. In general, the simulation results show that the availability of OpenFOAM and DualSPHysics as tools to predict ship resistance and their application to the development of ship design can potentially be leveraged by end-users. In this case study, the advantages of OpenFOAM are the reduced calculation time for simulations and good compactness of the results. On the other hand, DualSPHysics stands out for its simplicity of setting up the model and good reproduction of the free surface.