*Article* **Evaluation of the E**ff**ect of Container Ship Characteristics on Added Resistance in Waves**

#### **Ivana Marti´c 1, Nastia Degiuli 1,\*, Andrea Farkas <sup>1</sup> and Ivan Gospi´c <sup>2</sup>**


Received: 7 August 2020; Accepted: 6 September 2020; Published: 9 September 2020

**Abstract:** Added resistance in waves is one of the main causes of an increase in required power when a ship operates in actual service conditions. The assessment of added resistance in waves is important from both an economic and environmental point of view, owing to increasingly stringent rules set by the International Maritime Organization (IMO) with the aim to reduce CO2 emission by ships. For that reason, it is desirable to evaluate the added resistance in waves already in the preliminary ship design stage both in regular and irregular waves. Ships are traditionally designed and optimized with respect to calm water conditions. Within this research, the effect of prismatic coefficient, longitudinal position of the centre of buoyancy, trim, pitch radius of gyration, and ship speed on added resistance is investigated for the KCS (Kriso Container Ship) container ship in regular head waves and for different sea states. The calculations are performed using the 3D panel method based on Kelvin type Green function. The results for short waves are corrected to adequately take into account the diffraction component. The obtained results provide an insight into the effect of variation of ship characteristics on added resistance in waves.

**Keywords:** container ship; added resistance in waves; sea states; potential flow theory; variation of ship characteristics

#### **1. Introduction**

The ship hull is traditionally designed and optimized for calm water conditions, while an increase in the ship resistance owing to sailing in waves is being taken into account via sea margin. Added resistance in waves, as one of the additional loads acting on the ship in service, affects the attainable ship speed, and causes a change in the performance of the ship propulsion system. Its assessment is of great importance as early as in the preliminary design stage to adequately design and optimize the ship propulsion system. The assessment of added resistance ensures the ability of the ship to sail in severe sea states, as well as the possibility to estimate fuel consumption, with the aim of reducing operating costs and emissions of harmful gases. This is of particular importance from the environmental protection point of view and in accordance with regulations introduced by the International Maritime Organization (IMO). IMO has set a series of measures under the EEDI (Energy Efficiency Design Index) and EEOI (Energy Efficiency Operational Indicator) for new built ships and those already in service [1]. When sailing in waves, the total resistance can increase by 15–30% compared with calm water resistance, and still that increase could be even more pronounced [2]. For that reason, it is of particular importance to analyze the effect of variation of ship characteristics on the added resistance in waves. In that way, it is possible to evaluate a change in added resistance for different hull form characteristics rather than optimizing the hull for calm water conditions only. Moreover, it would enable the estimation of the effect of particular loading conditions on added resistance as early as in the preliminary design stage.

The additional power that a ship requires to attain speed when sailing in waves is related to three main components of added resistance. The first and largest one arises owing to the interference of waves generated by the ship response in waves (radiation waves) and incoming waves. It is often called the drift force, although the drift force in the longitudinal direction of the ship is equal to the added resistance in waves only in the case of zero forward speed. The largest contributors to radiation waves are heave and pitch motions and, to a certain extent, roll motion, depending on the phase shift regarding the incoming waves. In short waves, the component of added resistance caused by the ship motions tends to zero, but the second component, called the diffraction force, increases. The wave diffraction is a highly nonlinear phenomenon, thus the accuracy of the methods based on the potential flow theory decreases. The third component of the added resistance is related to the viscous damping of heave and pitch and is negligible compared with the hydrodynamic damping (radiation waves). Therefore, added resistance is considered as an inviscid phenomenon, which allows the application of methods and solvers based on the potential flow theory, as well as an extrapolation of the results from model to full scale, neglecting the scale effects. Added resistance is pressure driven force and can be extrapolated using Froude similarity law. The frictional part of added resistance, subjected to scale effects, is more emphasized in short incoming waves, when it can account for up to 20% of the added resistance. However, it should be emphasized that, in the latter case, added resistance accounts for only up to few percent of calm water resistance [3].

International maritime transport caused an average of 2.6% of global CO2 emission during the period from 2007 to 2015 [4,5]. For example, in 2015, ships burned about 298 million tons of fuel, of which 72% was heavy fuel oil (HFO), which caused the emission of about 932 million tons of CO2 into the atmosphere. Assuming that fossil fuels will remain dominant in the future, CO2 emissions from maritime transport are projected to increase by 50% to 250% by 2050, depending on economic and energy developments [4]. The containerization, which records the most significant increase in cargo transport, is the most energy efficient way to transport cargo and the share of CO2 emitted by ships is relatively low in global CO2 emission. However, there is a tendency to improve the energy efficiency of ships and limit greenhouse gas (GHG) emissions as maritime transport continues to evolve and increase. For this reason, IMO has set a goal in 2011 to reduce GHG emissions of ships by at least 50% by 2050 compared with 2008 by introducing mandatory technical measures for new built ships and operational measures for existing ships, with the aim to increase the energy efficiency of ships and reduce CO2 emissions [1]. EEDI of new built ships is the most important technical measure, aimed at promoting the use of energy efficient equipment and engines, as well as optimization of the propulsion system and ship hull [6]. In the optimization process, it is necessary to take into account the ship response as well as loads acting on the hull in waves. EEDI requires a minimum energy efficiency level per capacity mile for different ship types and sizes. It is expected to encourage innovation and technical development of all components of the ship system that affect fuel consumption, already in the preliminary design phase.

Considering the added resistance in waves as an inviscid phenomenon, numerical methods based on the potential flow theory can be successfully applied for its calculation. Bunnik et al. [7] have analyzed and compared the numerical results of motions and added resistance of the container ship and ferry with the available experimental data. Thus, the authors applied a linear potential flow theory based on the Green function with and without forward speed taken into account, that is, with approximate and exact forward speed, a nonlinear method based on Rankine singularities and computational fluid dynamics (CFD) based on the viscous flow theory. The authors concluded that the application of CFD requires incomparably greater computational resources and does not contribute significantly to the accuracy of the results when nonlinear effects are not pronounced. They also concluded that the methods based on a very similar or identical mathematical model give different results depending on the method of implementation, discretization, or boundary conditions, and that

it is more demanding to obtained accurate results for heave than for pitch. Hong et al. [8] have divided the added resistance of S175 container ship into a radiation and diffraction part and applied the Green function based on pulsating and translating sources to solve the boundary condition problem. The applied method provided greater accuracy compared with the experimental results at higher Froude numbers than the Rankine panel method. Linear potential flow theory has been applied to calculate added resistance of intact and damaged S175 ships in [9]. For considered sea states, two calculation models, that is, flooding simulated as increased displacement mass and sloshing inside flooded tank, gave very similar results. Riesner and el Moctar [10] have developed a partially nonlinear time domain method for predicting added resistance at a constant forward speed in regular waves. For the calculation of radiation forces in the time domain, the developed method is based on hydrodynamic coefficients obtained in the frequency domain. The pressure is integrated on the instantaneous wetted surface, which enables the determination of nonlinear Froude–Krylov force and restoring forces. Moreover, an additional force is applied to take into account a variation in wetted surface owing to the radiation and diffraction of incident waves. The additional term, added to take into account the viscous component of the added resistance, is in a form of an exponential function depending on the block coefficient, ratio of wavelength to ship length, and constants determined by the least squares method based on the available results from the literature. Yang et al. [11] have improved the Faltinsen asymptotic formula for calculating the added resistance in short waves. The authors took into account the ship draft, the alteration in the flow velocity around the hull, and the shape of the hull above the waterline. Furthermore, the authors introduced a correction factor associated with a variation in the hull cross section below the waterline, and variation in the hull shape above the waterline was taken into account based on the waterline bluntness coefficient. As new built ships are getting larger, the ratio of wavelength to ship length decreases, meaning that it becomes very important to estimate the added resistance in short waves as accurately as possible. Whether based on potential or viscous flow theory, numerical tools for calculating added resistance in short waves require a large number of panels or finite volumes in order to capture alterations in the flow around the fore part of a hull. Liu and Papanikolaou [12] have developed a relatively simple method for the determination of added resistance in accordance with IMO-MEPC.232(65) EEDI recommendations [13] for establishing the so-called level 1 methods for estimating the minimum required power for navigation and maneuverability of ships in severe weather conditions. The developed empirical method requires only main ship and wave characteristics as input values. The proposed method was extended and adapted for a wider speed range and variation of draft and trim at different loading conditions [14]. The authors analyzed the influence of the advancing speed, the pitch radius of gyration, the draft and the trim on the amplitude and position of the peak value of the added resistance in head waves. They concluded that the peak position is closely related to the natural frequency of heave in long waves and that the amplitude depends on the radius of gyration. Seo et al. [15] have numerically calculated added resistance in short waves using methods based on the potential and viscous flow theory for different ship forms. Given that the added resistance in short waves is greatly influenced by the size of the panel or finite volume, especially for full hull forms, the authors concluded that a convergency study is necessary. Liu and Papanikolaou [16] have utilized the 3D panel method to calculate second-order sway force, that is, drift force for various wave headings at low forward speeds, and concluded that the mean sway force is most significant in relatively short waves. Guha and Falzarano [17] have developed a numerical method for the determination of added resistance based on the direct pressure integration along the instantaneous wetted surface, taking into account the angle of bow flare above the linearized free surface. The authors concluded that it is of great importance to take into account the hull immersion angle when calculating the relative wave amplitude along the waterline on the example of container ships S175 and KCS (Kriso Container Ship), while it does not significantly affect the results of the analyzed Ro-Ro ship. Yang and Kim [18] have come to the similar conclusion by analyzing the bow shape of KVLCC2 (Korean Very Large Crude Carrier) above the waterline. Park et al. [19] have investigated the effect of draft on added resistance in waves experimentally and numerically

on the example of KVLCC2. On the basis of the experimental results, the authors concluded that, for short wave lengths, the largest added resistance occurs in the ballast condition and that the peak value shifts towards higher frequencies when reducing the draft. The influence of bow flare on added resistance in waves has been investigated by Fang et al. in [20]. The authors performed numerical calculations, based on the potential flow theory in time domain, and concluded that, by increasing the bow flare angle, the ship motion amplitudes decrease, while added resistance in waves increases. Trim optimization has been performed by various authors for calm water conditions [21,22]. However, the literature lacks studies regarding the added resistance in waves for different trim conditions and trim optimization for sailing in waves. Jung and Kim [23] have applied the multi-objective optimization method to minimize the ship resistance and speed loss caused by wind and waves by varying main dimensions and prismatic coefficient of full KVLCC2 hull form. Kim et al. [24] have proposed a reliable methodology for the estimation of speed loss of the S175 container ship for various sea states and wind conditions. The validation study of ship motions and added resistance obtained using 2D and 3D linear potential flow methods and CFD showed reasonably good agreement with the experimental data in regular head and oblique waves. The authors concluded that the estimated sea margin is significantly increased when the initial reference speed is decreased, mainly owing to added resistance in waves, despite the fact that the absolute value of the required increase in power is lower.

Within this paper, the effect of hull form characteristics of KCS, that is, prismatic coefficient and longitudinal position of the centre of buoyancy, as well as pitch radius of gyration, trim, and sailing speed on added resistance in regular waves and for different sea states is evaluated. To the best of the authors' knowledge, this problem has not been investigated in such a comprehensive way. Considering that added resistance in waves could cause a significant increase compared with calm water resistance, the obtained results provide a valuable insight into the possibility to reduce the added resistance in waves for actual sea states that the ship will encounter during her service. During the design process, certain ship characteristic could be optimal with respect to the calm water resistance, but could disrupt seakeeping characteristics and cause an increase in motion amplitudes or added resistance in waves. The obtained results show the effect of variation of certain ship characteristics on added resistance in waves, on the example of a typical container ship. Added resistance is calculated utilizing the 3D panel method based on the Kelvin type Green function. The numerically obtained results are validated against the experimental data available in the literature and the convergence study is performed to assess the numerical uncertainty.

#### **2. Methodology**

#### *2.1. Case Study*

Calculations of added resistance in regular and irregular waves are performed for the benchmark container ship, KCS, Figure 1. The main particulars of KCS are shown in Table 1. The original hull form is modified in terms of variation of the prismatic coefficient and the position of the longitudinal centre of buoyancy (LCB) by shifting the hull cross sections in the longitudinal direction, while the midship coefficient and main dimensions remain constant. Cross sections of original and modified hull forms have the same shape and area, but different longitudinal positions. The hull modification method can be applied to ships with and without a parallel mid-body, with the aim of adding a parallel mid-body (if the ship does not have one) or changing its length, while keeping the original value of the prismatic coefficient constant or varying its values [25]. The limits for the variation of the prismatic coefficient and position of LCB are defined according to [25–27]. In order to change the prismatic coefficient of the KCS hull form, cross sections are shifted in the longitudinal direction, without affecting the position of LCB. On the other hand, while changing the position of LCB, an original value of the prismatic coefficient is kept constant.

**Figure 1.** 3D model of Kriso Container Ship (KCS).

**Table 1.** The main particulars of Kriso Container Ship (KCS).


Cross sections of an example of a modified hull form of KCS with prismatic coefficient equal to 0.688 are shown in Figure 2. In order to analyze the effect of gyration radii on added resistance in waves, the pitch (and yaw) radius of gyration is set as 24%, 25%, and 26% of the length between perpendiculars. The roll radius of gyration is set as 35% of the ship beam. The variation of trim is performed for constant displacement volume with step equal to 0.2◦ with positive trim defined as bow up and negative trim as bow down. Added resistance in waves for the original hull form of KCS is calculated for three different speeds.

**Figure 2.** Cross sections of KCS with *CP* = 0.661 (black line) and *CP* = 0.680 (blue line).

#### *2.2. Numerical Assessment of Added Resistance in Waves*

Within this research, added resistance is numerically calculated utilizing the 3D panel method, by satisfying the Kelvin boundary condition on a free surface. The ship hull is discretized by quadrilateral flat panels with sources of constant and unknown strength located in the collocation points. The velocity field in the computational domain is described by a scalar function of the velocity potential Φ and the Green theorem is applied to obtain the velocity potential by distributing the singularities over the domain boundaries. Boundary integral equations (BIEs) are solved to determine the strength of distributed sources, in a way that the Laplace equation and all the necessary boundary conditions are satisfied. By developing the velocity potential above mean free surface and linearizing

the equations under the assumption of small amplitudes and wave steepness, that is, excluding the square terms, the boundary value problem (BVP) of the first order is defined as follows [28]:

$$
\nabla \mathbf{V} = \nabla^2 \Phi = 0,\tag{1}
$$

$$
\delta\_{\mathcal{S}} \Phi\_z + \Phi\_{tt} + \mu \Phi\_t = 0 \text{ for } z = 0,\tag{2}
$$

$$
\boldsymbol{\Phi}\_{\boldsymbol{n}} = \mathbf{X}\_{\boldsymbol{t}} \mathbf{n} \quad \text{on wetted surface}\_{\boldsymbol{t}} \tag{3}
$$

$$
\Phi\_z = 0 \quad \text{for } z = -h,\tag{4}
$$

where **X***<sup>t</sup>* is the velocity vector; **n** is the normal vector; and indices *t* and *z* denote time and space derivatives, respectively. μ represents a positive and small parameter that introduces the energy dissipation into the momentum equation by a fictitious force dependent on the flow velocity without affecting the inviscid and irrotational properties of the fluid. The result is an additional term in the free surface boundary equation in order to prevent an infinite response at resonant frequencies.

The radiation boundary condition, which requires the velocity potential to disappear at an infinite distance away from the ship hull, is automatically satisfied in the fairly perfect fluid. All velocity potential components satisfy the radiation condition except for the incoming wave potential. Wave elevation ζ evaluated at mean free surface and pressure *p* on hull are given as follows:

$$
\mathbb{X}\mathbb{g} = -\Phi\_t - \mu\Phi,\tag{5}
$$

$$\frac{p}{\rho} = -\mathbf{g}\mathbf{X}\_3 - \Phi\_t - \mu\Phi\_\prime \tag{6}$$

where *g* is the gravity constant and ρ is the fluid density.

The Green function, G, as a fundamental solution of the Laplace equation, is formulated in such way that it satisfies the linearized boundary condition on the free surface, boundary condition on the bottom, and the radiation condition as follows:

$$\nabla^2 \mathcal{G}(P, Q, t) = 4\pi \delta(P - Q) \, , \tag{7}$$

$$
\mathcal{g}\mathcal{G}\_z + \mathcal{G}\_{tt} + \mu'\mathcal{G}\_t = 0 \quad \text{for } z = 0,\tag{8}
$$

$$\mathcal{G}\_z = 0 \quad \text{for } z = -h,\tag{9}$$

The Green function represents the field of velocity potential in the computational domain at *P*(*x*, *y*, *z*) created by the source of unit strength located at *Q*(*x* , *y* , *z* ). δ(*P* − *Q*) = δ(*x* − *x* )δ(*y* − *y* )δ(*z* − *z* ) is the Dirac function and μ = μ. The free surface Green function is expressed by the Fourier–Hankel integral and approximated by Chebyschev polynomials. More details can be found in [28]. After applying the Green theorem to the domain boundaries (mean free surface, hull surface, sea bed, and cylindrical surface at infinity) and considering a complementary domain inside the ship limited by hull and interior free surface, the integral equation is defined as follows:

$$4\pi\Phi(P) = \iint\limits\_{H} [(\Phi\_{\text{il}} - \Phi\_{\text{il}}')\mathcal{G} - (\Phi - \Phi')\mathcal{G}\_{\text{il}}] \,\text{d}s,\tag{10}$$

where Φ represents the velocity potential in the interior domain. With the velocity potential and the Green function expressed as Φ(*P*, *t*) = Re, φ(*P*)e−iω*<sup>t</sup>* and <sup>G</sup>(*P*, *<sup>Q</sup>*, *<sup>t</sup>*) <sup>=</sup> Re, *G*(*P*, *Q*)e−iω*<sup>t</sup>* - , respectively; and velocity potential expressed by radiation φ*j*, diffraction φ7, and incoming wave potential φ<sup>0</sup> as <sup>φ</sup> <sup>=</sup> <sup>−</sup>i<sup>ω</sup> .<sup>6</sup> *j*=1 ζ*aj*φ*<sup>j</sup>* + ζ*a*0(φ<sup>0</sup> + φ7), with φ = φ and with source strength defined as σ = φ*<sup>n</sup>* − φ *n*. Moreover, by satisfying the boundary condition on the hull (derivation of the velocity potential in the normal direction), the strength of distributed sources can be determined as follows:

$$2\pi\sigma\_{\hat{\jmath}} + \iint\_{H} \sigma\_{\hat{\jmath}} \mathbf{G}\_{\text{n}} \, \text{ds} = \begin{cases} n\_{\hat{\jmath}} & j = 1 \dots 6 \\ \frac{-\partial \phi\_0}{\partial \mathbf{n}} & j = 7 \end{cases} \tag{11}$$

With the known velocity potential, it is possible to determine the pressure on each panel from the Bernoulli equation. When the point at which the velocity potential is determined is located at a singular point on the hull (*P* = *Q*), the gradient of the Green function becomes singular. In that case, the small area around the singular point is excluded from the integration. For this reason, a term 2πσ*<sup>j</sup>* is added to Equation (11) that contributes to the amount of normal derivative on the hull in the vicinity of the singular point on the hull. Integral Equation (11) has a unique solution except for the frequencies when the determinant disappears. As the Green function satisfies the free surface boundary condition in the entire computational domain, the solution in the exterior domain of the body is obtained simultaneously as a fictitious solution in the complementary domain. That may cause numerical error at so-called irregular frequencies corresponding to the eigenvalues of the homogenous Dirichlet problem. With φ = φ at these frequencies, the interior BVP has a nontrivial solution. One of the possible solutions to this problem is to define a different boundary condition on the interior free surface, and thus achieve a unique solution for all frequencies. The extended BIEM imposes Neumann's "rigid lid" condition on the interior free surface. In that way, an extended boundary condition is introduced without changing or adding any more parameters to the formulation. For this reason, it is necessary to adequately discretize the interior free surface as well, depending on the frequency range in which the occurrence of the first irregular frequency is expected. The integral equations on the hull and interior free surface read as follows:

$$2\pi\sigma(P) + \iint\limits\_{H\cup\mathcal{F}'} \sigma(Q)G\_{\mathcal{U}}(P,Q)\mathrm{d}\mathbf{s} = \begin{cases} \begin{array}{c} n\_j \end{array} & j = 1\ldots 6\\ -\partial\phi\_0/\partial\mathbf{n} \text{ j} = 7 \end{array} \tag{12}$$

$$4\pi\sigma(P) - \iint\limits\_{H\cup\mathcal{F}'} \sigma(Q) G\_{\mathcal{U}}(P, Q) \mathrm{d}s = 0,\tag{13}$$

More details on the extended BIEM can be found in [29], where it was shown that the second-order loads are much more sensitive to the irregular frequencies compared with the first-order quantities and that the effects of irregular frequencies should be removed. The second-order wave loads can be determined by integrating the second-order pressure over the mean wetted surface, taking into account the change in the first-order quantities owing to ship motions. In that way, second-order loads can be expressed through one part depending on the variables of the first order and the other part depending on the second-order velocity potential. The ship response to the incoming waves causes a shift from the equilibrium position and a change in the pressure distribution over the wetted surface. Nonlinear forces of a higher order are caused by the change in wetted surface due to the incoming waves, but also as a result of the ship response in waves (especially angular motion), as well as the nonlinear pressure term in the Bernoulli equation. If the force up to the second order is considered, it contains both a constant and an oscillatory part dependent on the square of the wave amplitude. The result of the direct pressure integration along the mean wetted surface and waterline is the quadratic transfer function (QTF) of high-frequency and low-frequency wave loads. The constant drift force, that is, the added resistance, is a time-averaged value of the second-order force, which amounts to about 5% of the first-order wave force and requires only the calculation of the first-order velocity potential. It is represented by the diagonal members of QTF calculated as follows [30,31]:

$$\mathbf{F}\_{1} = -\frac{1}{2}\rho g \oint\_{\text{WL}} \zeta\_{r}^{2} \mathbf{n} \,\mathrm{dl} + \frac{1}{2}\rho \iint\_{H} \nabla\Phi \nabla\Phi \,\mathbf{n} \,\mathrm{ds} + \rho \iint\_{H} \mathbf{X} \nabla\Phi\_{l} \mathbf{n} \,\mathrm{ds} + mR\ddot{\mathbf{X}}\_{\text{G}}.\tag{14}$$

where ζ*<sup>r</sup>* = ζ − *X*<sup>3</sup> is the relative wave elevation defined as coupled wave elevation and ship vertical motion and **X** = **X***<sup>G</sup>* + *R***x** is the motion vector of the first order in the ship coordinate system, where the position vector of the point on the hull is given as **x** and linearized rotation transformation matrix as *R*. Equation (14) includes the ship response, the first-order velocity potential, the gradient of the velocity potential, and the wave elevation on mean wetted surface and along the waterline. The first term on the right side of the equation refers to the relative wave elevation, that is, to the integration of the first-order pressure along the variable part of the wetted surface. If the hull emergence angle γ has a non-zero value in the waterline area, that is, the hull surface is not vertical, which is the case for most modern hull forms in the bow and stern area, the normal vector in the integral equation along the waterline needs to be modified according to **n** = **n**/ cos <sup>γ</sup> [32]. The second and third term in Equation (14) refer to the integration of the first-order pressure along the mean wetted surface, taking into account pressure correction due to ship displacement. The last term refers to the rotation of the first-order force vector and its horizontal component owing to pitch. Namely, first-order hydrodynamic and hydrostatic forces acting along the vertical axis of ship coordinate system will give a longitudinal second-order force component equal to *<sup>x</sup>*5*F*<sup>3</sup> <sup>=</sup> *<sup>x</sup>*5*<sup>m</sup>* .. *X*3*G*.

A commercial software package HydroSTAR [33] based on linear potential flow theory is used to perform hydrodynamic calculations. HydroSTAR provides a solution of the first-order wave diffraction and radiation problem, and second-order wave loads with and without forward speed. The immersed part of the ship hull and interior free surface are discretized by quadrilateral panels. Every panel model has approximately 4000 panels on the hull and 6000 panels on the interior free surface. The size of the panels on the interior free surface is equal to 50% of the panel size on the hull in order to shift irregular frequencies further towards a higher frequency range.

The reference coordinate system is located on the mean free surface at the stern with the positive direction of *z* axis upwards. Within HydroSTAR, mesh is generated using the so-called automatic mesh generator (AMG) based on an adaptive cosine rule, which reduces the panel size from the bottom of the ship towards the waterline, allowing for finer discretization in the area of the waterline. The pressure on the panels closest to the waterline is used to determine the wave elevation using the boundary condition on the free surface. The coordinate system for solving the diffraction and radiation problems is located at the ship centre of buoyancy, and the solution is transformed to the centre of gravity when solving the motion equation.

#### 2.2.1. Correction of the Results in Short Waves

Generally, added resistance in waves calculated by linear potential flow theory significantly underestimates the experimental results for higher incoming wave frequencies. In short waves, where the diffraction force is dominant, linear potential flow theory fails to properly include the diffraction component owing to highly non-linear effects. The results in short waves are thus corrected according to [34] based on the experimental data. With the correction term included, added resistance is calculated as follows:

$$R\_{AW} = R\_{AWm} + R\_{AWr\_{\prime}} \tag{15}$$

where the correction in short waves depends on the bluntness coefficient *Bf* , draft and incoming wave frequency α*d*, and forward speed (1 + α*U*) as follows:

$$R\_{AWr} = \frac{1}{2} \rho g \zeta\_a^2 B B\_f \alpha\_d (1 + \alpha\_{lI}) \, \tag{16}$$

$$B\_f = \frac{1}{B} \left\{ \int\_I \sin^2(\alpha + \beta\_w) \sin \beta\_w \mathrm{d}l + \int\_{\mathrm{II}} \sin^2(\alpha - \beta\_w) \sin \beta\_w \mathrm{d}l \right\},\tag{17}$$

$$\alpha\_d = \frac{\pi^2 I\_1^2(k\_\epsilon T)}{\pi^2 I\_1^2(k\_\epsilon T) + K\_1^2(k\_\epsilon T)},\tag{18}$$

$$k\_t = k \left( 1 + \frac{\alpha V}{\mathcal{g}} \cos \alpha \right)^2,\tag{19}$$

$$1 + a\_{\rm II} = 1 + \mathbb{C}\_{\rm II} \text{F} \mathbf{n}\_{\prime} \tag{20}$$

$$\mathcal{C}\_{L}(\alpha) = \max\left[10, -310\mathcal{B}\_f(\alpha) + 68\right],\tag{21}$$

where α is the wave heading, β*<sup>w</sup>* is the slope of line element along waterline, *I* and *II* are the domains of the integration, *I*<sup>1</sup> is the first-order Bessel function, *K*<sup>1</sup> is the second-order Bessel function, *k* is the wave number, *ke* is the encounter wave number, ω is the wave frequency, and *CU* is the forward speed coefficient.

#### 2.2.2. Added Resistance in Irregular Waves

The mean value of added resistance in irregular waves can be determined based on the spectral moment *m*0*<sup>R</sup>* or integral of the response spectrum curve. The response spectrum is determined as the product of the added resistance transfer function *RAW* ζ2 *a* (ω*e*) in regular waves and the encounter wave energy spectrum *S*ζ(ω*e*) as follows:

$$\overline{R\_{AW}} = 2 \int\_0^\infty S\_\zeta(\omega\_\ell) \frac{R\_{AW}}{\zeta\_\pi^2} (\omega\_\ell) \mathrm{d}\omega\_\ell = 2m\_{0\mathcal{R}\_\ell} \tag{22}$$

where ζ*<sup>a</sup>* is the wave amplitude.

In order to cover the wave energy distributed in the range of encounter frequencies, the wave energy spectrum is transformed into the wave encounter spectrum based on the encounter frequency <sup>ω</sup>*<sup>e</sup>* <sup>=</sup> <sup>ω</sup> <sup>−</sup> <sup>ω</sup><sup>2</sup> *<sup>g</sup> V* cos α as follows:

$$S\_{\zeta}(\omega\_{t}) = \frac{S\_{\zeta}(\omega)}{\sqrt{1 - \frac{4\alpha\_{t}V\cos\alpha}{\mathcal{J}}}},\tag{23}$$

It should be noted that Equation (23) is valid only for head and oblique waves. Namely, the following wave encounter spectrum can be double valued because the same encounter frequency corresponds to different incoming wave frequencies. In addition, as the ship speed increases, the wave encounter spectrum extends into the complex solutions range. Within this research, for fully developed sea, a modified two-parameter Pierson–Moskowitz or Bretschneider wave energy spectrum recommended by ITTC (International Towing Tank Conference) is used to calculate added resistance in head waves at different sea states [35]. The Bretschneider wave energy spectrum is defined as follows:

$$S\_{\zeta B}(\omega) = \frac{173H\_{S1/3}^2/\overline{T}^4}{\omega^5} \text{e}^{-\frac{602/\overline{T}^4}{\omega^6}},\tag{24}$$

With the characteristic wave period expressed through zero crossing period as *T* = 1, 086 *Tz*. *HS*1/3 is the significant wave height.

For limited fetch, the modified JONSWAP (Joint North Sea Wave Project) wave energy spectrum is used:

$$S\_{\zeta l}(\omega) = \frac{320H\_{S1/3}^2}{T\_p^4 \alpha^5} \mathbf{e}^{\frac{1050}{T\_p \omega^4}} \boldsymbol{\gamma}^{\left| - \left( \frac{\omega/\omega\_p - 1}{\omega \sqrt{2}} \right)^2 \right|}, \quad \boldsymbol{\sigma} = \begin{cases} 0, 07 \text{ za } \omega < \omega\_p \\\ 0, 09 \text{ za } \omega > \omega\_p \end{cases} \tag{25}$$

where ω*<sup>p</sup>* = <sup>2</sup><sup>π</sup> *Tp* , γ = 3, 3, and peak wave period is defined as 1, 073*Tz* = 0, 834*Tp*.

#### *2.3. Validation and Verification Study*

Numerically obtained results of the added resistance in waves are validated against the available experimental data for KCS [36–38] for *Fn* = 0.26, which corresponds to a speed of 24 knots. Relative deviation of the numerical results is calculated as follows:

$$RD = \frac{\mathbb{C}\_{AW, \text{HSTAR}} - \mathbb{C}\_{AW, \text{EXP}}}{\mathbb{C}\_{AW, \text{EXP}}} 100\%,\tag{26}$$

where added resistance coefficient is defined as follows:

$$\mathbb{C}\_{AW} = \frac{R\_{AW}}{\rho g \zeta\_a^2 B^2 / L}. \tag{27}$$

The mesh convergence study, as part of the verification procedure [39], is performed to evaluate the effect of the discretization on the obtained numerical results of added resistance, heave, and pitch amplitudes, as well as to quantify the numerical uncertainty. Although it seems that, using a larger number of panels, with sources of constant strength located in collocation points, a more accurate solution can be achieved, this may not be a case. However, a larger number of panels allows more detailed description of the hull geometry, which is of particular importance in the bow and stern area, where significant changes in pressure occur, Figure 3. Generated panel models have 1048, 4366, and 17,823 panels in total. In order to preserve the geometrical similarity between the panel models, each panel is divided into four panels when creating a model with finer discretization.

**Figure 3.** Panel models of KCS: coarse mesh (upper) and fine mesh (lower) [33].

The numerical uncertainty of the obtained results is performed according to the verification procedure recommended by ITTC [40] to estimate the uncertainty of numerical results obtained by viscous flow theory, but it can also be applied to the results of potential flow theory [41]. The numerical uncertainty is evaluated for added resistance coefficient in regular waves and irregular waves by applying the Bretschneider wave energy spectrum in the frequency range from 0.3 to 0.8 rad/s for sea state defined with *Hs* = 3.5 m and *TZ* = 10.5 s. According to [40], the type of convergence is determined based on the change between solutions obtained using medium and fine mesh (φ<sup>2</sup> − φ1) and coarse and medium mesh (φ<sup>3</sup> − φ2) as follows:

$$R = \frac{(\phi\_2 - \phi\_1)}{(\phi\_3 - \phi\_2)},\tag{28}$$

Monotonic convergence is achieved when 0 < *R* < 1, oscillatory convergence when −1 < *R* < 0, and divergence when |*R*| > 1. Numerical uncertainty can be calculated for monotonic convergence case using generalized Richardson extrapolation (RE). The order of accuracy is calculated as follows:

$$p = \frac{\ln\left(\left(\phi\_3 - \phi\_2\right)/\left(\phi\_2 - \phi\_1\right)\right)}{\ln(r)},\tag{29}$$

where *r* = *<sup>h</sup>*<sup>2</sup> *<sup>h</sup>*<sup>1</sup> *h*<sup>3</sup> *<sup>h</sup>*<sup>2</sup> is the uniform refinement ratio corresponding to the ratio of panel size determined based on the total number of panels as *hi* <sup>=</sup> 1/ <sup>√</sup> *Ni* [42,43]. With safety factor *FS* equal to 1.25, the normalized numerical uncertainty is calculated as follows:

$$
\overline{\mathcal{U}} = \frac{F\_S \left| \delta\_{RE} \right|}{\phi\_1 - \delta\_{RE}} \cdot 100,\tag{30}
$$

where RE error is defined as follows:

$$
\delta\_{RE} = \frac{\phi\_2 - \phi\_1}{\binom{h\_2}{h\_1}^p - 1}.
\tag{31}
$$

#### **3. Results and Discussion**

#### *3.1. Validation of the Numerical Results*

The obtained numerical results show satisfactory accuracy compared with the experimental data, except in the short wave region, Figure 4. After the correction of the results for high wave frequencies has been applied, significantly better agreement with the experimental data is obtained. The correction of the results, which accounts for the diffraction part of added resistance, was applied for the results up to λ/*L* = 0.85, because it may cause an overestimation of the numerically obtained added resistance values in moderate and long incoming waves. Some discrepancy in the experimental results can be seen as well, especially in the short wave region. By comparing the ship transfer functions of heave and pitch with the experimental data available in the literature, it can be noticed that the pitch transfer function shows much better agreement with the experimental results, with the relative deviation below 5%. The numerically obtained values of heave for moderate and long waves significantly exceed the experimental values. Relative deviation of heave at λ/*L* = 1474 is equal to 29%, which confirms that ensuring the accuracy of heave is more demanding than pitch, Figure 5.

**Figure 4.** Comparison of the numerically obtained added resistance coefficient with available experimental data [36–38].

**Figure 5.** Comparison of the numerically obtained heave (**left**) and pitch (**right**) transfer functions with available experimental data [36–38].

#### *3.2. Verification of the Numerical Results*

Table 2 shows the hydrostatic characteristics of KCS obtained using panel models with different mesh density. The differences between the hydrostatic characteristics for different panel models are not significant, even though the displacement volume of the finest discretization is closest to the actual value. It should be noted that the position of the vertical centre of buoyancy is given with the respect to the coordinate system located on the mean free surface. The effect of panel model discretization on the obtained results of heave and pitch motion amplitudes can be seen in Figure 6. Increasing the number of panels does not influence the obtained results, unlike in the case of added resistance. The largest change in the results is obtained for the highest frequency using fine mesh. Without the correction being applied, this value has the lowest relative deviation compared with the experimental results. The numerical uncertainty of the obtained added resistance coefficient in regular waves and coefficient of the mean added resistance in irregular waves can be seen in Table 3. It can be seen that the calculated numerical uncertainty, for the results with monotonic convergence, is below 2.1% for regular and equal to 1.64% for irregular waves.

**Table 2.** Hydrostatic characteristics of KCS obtained using different panel models.


**Figure 6.** Heave (**left**) and pitch (**right**) transfer functions obtained using different panel models.


**Table 3.** Convergence study for panel size used to discretize KCS.

#### *3.3. Results of the Sensitivity Analysis*

The effect of the prismatic coefficient of KCS on added resistance in regular and irregular waves (Table 4) is presented by the results in Table 5. The results are given as relative deviations with respect to the results obtained for the original hull form of KCS. It can be seen that, in long and moderate regular waves, except for the wave frequency corresponding to the peak value of added resistance, the increase in prismatic coefficient and block coefficient at the same time causes the decrease in added resistance. Namely, for large wavelengths, without the presence of diffraction component, an increase in the prismatic coefficient leads to smaller amplitudes of the ship motions. The relative motions become large when the length of the ship is approximately equal to the wavelength, causing the largest resistance. In very long waves, the relative motions tend to zero, and the force tends to a Froude–Krylov force that would ideally act on the ship in the absence of a diffraction component. Consequently, in very long waves, the added resistance tends to zero. The peak value of the added resistance decreases with the decrease in the prismatic coefficient. In order to analyze the effect of variation of the prismatic coefficient on added resistance in short waves, a correction for diffraction was applied. It can be noticed that the diffraction part of added resistance in short waves increases as the prismatic coefficient increases, which means that the fore part of ship is fuller. For irregular waves, the mean value of added resistance was calculated for six sea states by means of spectral analysis based on two theoretical wave energy spectra, Table 4. Spectral analysis based on the corrected results showed that added resistance in irregular waves decreases for hull modifications with a lower value of the prismatic coefficient, and vice versa, and that a change in added resistance is more pronounced for lower sea states when wave energy is concentrated in the high frequency range. In that way, the diffraction component of added resistance determines the trend of mean value in irregular waves. At higher sea states, with wave energy distributed along the entire frequency range, the obtained results are mainly within numerical uncertainty and the results in regular waves counteract each other.


**Table 4.** Significant wave height and wave period for considered sea states.


**Table 5.** Influence of the prismatic coefficient on added resistance in regular and irregular waves. JONSWAP, Joint North Sea Wave Project.

In Table 6, the obtained results of the effect of the longitudinal centre of buoyancy on added resistance in regular and irregular waves are shown. In order to change the position of *LCB*, for example, towards the bow, cross sections of the fore part were moved towards bow and cross sections of the aft part towards midship. In that way, the fore prismatic coefficient increases, while the aft one decreases, keeping the original value of the total prismatic coefficient constant. Two *LCB* positions are analyzed and the obtained numerical results are given as relative deviations with regard to results of the original KCS hull form. On the basis of the obtained numerical results, it is possible to conclude that the position of *LCB* does not significantly affect added resistance in long regular waves. It can be seen that the peak value of added resistance decreases by almost 6% with *LCB* shifted towards stern, and increases by over 3% with *LCB* shifted towards bow. In general, shifting the position of *LCB* towards the stern allows a fine shape of the fore part of the ship. On the other hand, shifting the position of *LCB* towards the bow increases the distance between the longitudinal centre of buoyancy and longitudinal centre of floatation, which reduces the ship relative motions, causing the decrease in added resistance in relatively short waves. This can be observed for λ/*L* = 0.5, where added resistance increases when *LCB* shifts towards stern. A noticeable effect of *LCB* on the mean value of added resistance in irregular waves can be observed for higher sea states, when the position of *LCB* is shifted towards the stern. The decrease in added resistance in that case equals to 3.5–5% for both theoretical wave energy spectra.

To perform seakeeping calculations, it is necessary to know the mass characteristics along with the hull form characteristics. As shown in [44], the influence of the vertical position of the centre of gravity on added resistance in head waves is almost negligible when calculated using linear potential flow theory and by keeping the values of gyration radii constant. Changing the position of the centre of gravity does not affect the terms of hydrodynamic mass matrix, given that the coordinate system is placed at the centre of gravity. On the other hand, an insignificant change in terms of hydrodynamic added mass and damping matrix, especially in the case of heave and pitch, is noticed. In the restoring forces matrix, the large value of longitudinal metacentric height changes only slightly by shifting the vertical position of the centre of gravity.


**Table 6.** Influence of the position of the longitudinal centre of buoyancy on added resistance in regular and irregular waves.

Pitch gyration radius in head waves has a more significant effect on added resistance in waves. As mentioned before, its value is varied as a percentage of ship length between perpendiculars. The initial value of pitch (and yaw) gyration radius for KCS container ship corresponds to 25% of the length between perpendiculars. Additionally, values corresponding to 24% and 26% of the length between perpendiculars are analyzed, Figure 7. It can be seen that, by increasing and decree asing thpitch gyration radius, added resistance in long and moderate waves increases and decreases, respectively, by approximately 15%, while the change in peak value amounts to about 6%. According to the obtained numerical results in short waves, there is no effect of change on the pitch gyration radius, which is expected considering that the amplitudes of absolute ship motions in that frequency range are small. By changing the pitch radius of gyration, it can be seen from Figure 8 that added resistance in irregular waves decreases or increases by up to 10% for higher sea states. The change in added resistance is more pronounced for the Bretchneider than for the JONSWAP sea spectrum, except for SS6, when the wave energy described by the Bretschneider spectrum shifts towards lower wave frequencies. The narrow banded JONSWAP spectrum has a higher peak value in the range of high wave frequencies, causing a larger change in mean added resistance for SS1, even though the mean value of added resistance in that case equals approximately 20 kN.

Figure 9 shows the obtained numerical values of added resistance in regular waves for three speeds: 20, 22, and 24 knots, which correspond to Froude numbers of 0.216, 0.238, and 0.260, respectively. Compared with the design speed equal to 24 knots, at the frequency corresponding to the peak value of added resistance, a decrease of speed by 2 knots causes a decrease in added resistance by 15%, and the decrease by 4 knots causes the decrease by an additional 14%. This decrease is even more pronounced at lower wave frequencies. For example, a decrease in added resistance by decreasing the speed by 2 and 4 knots is about 15% and 28% for λ/*L* = 1.33 and 22% and 36% for λ/*L* = 1.5, respectively. It can also be seen that, as the speed increases, the peak value of the added resistance slightly shifts towards lower frequencies. The corrected results in short waves show a significant decrease in added resistance for lower speeds: 25% and 35% for 22 knots and 20 knots, respectively.

**Figure 7.** Influence of pitch gyration radius on added resistance coefficient in regular waves.

**Figure 8.** Influence of pitch gyration radius on added resistance coefficient in irregular waves: *ryy*/*LPP* = 0.24 (**left**), *ryy*/*LPP* = 0.26 (**right**).

**Figure 9.** Influence of speed on added resistance coefficient in regular waves.

The decrease in the mean values of added resistance in irregular waves with decrease in speed is shown in Figure 10. It can be seen that, for lower sea states, the decrease in added resistance is more pronounced in the case of the JONSWAP spectrum. Even though a decrease in percentage of added resistance is more pronounced for lower sea states, the absolute differences between the added resistances are much larger at higher sea states. Decreasing the speed for two knots can lead to an approximately 10% reduction in added resistance based on the results obtained using the Bretschneider wave energy spectrum. The obtained numerical results of the effect of ship speed, as well as pitch gyration radius on added resistance in regular waves, show similar trends as in [14].

**Figure 10.** Influence of speed on added resistance coefficient in irregular waves: 20 knots (**left**) and 22 knots (**right**).

Figure 11 shows the obtained added resistance coefficients for different trim conditions at the investigated λ/*L*. Trim can significantly affect the fuel consumption during the operation and is often optimized in the design stage. In that way, finding the optimal trim is enabled, taking into account that the ship's weight distribution affects the trim. As already mentioned, trim is varied from 1◦ to −1◦ with a step equal to 0.2◦. The obtained numerical results are not reliable for a ship with trim equal to 1◦. The values of added resistance for that particular case show unphysical character. Namely, the discretization of the bow, with mean free surface very close to the bulb, in the most prominent trim condition with bow up, leads to almost horizontally positioned flat panels and singularities very close to the free surface, causing the numerical error. With that particular trim condition excluded, it appears that the optimum trim condition in short waves corresponds to 0.8◦, that is, bow up. For λ/*L* = 1.15, which corresponds to the peak value of added resistance, there is no significant decrease or increase in added resistance for a certain trim condition, except a significant increase equal to 12% for trim of 0.8◦. For trim equal to 0.6◦, the decrease in added resistance equal to 7.3% can be noticed for λ/*L* = 1.5. For λ/*L* = 1.33, there is a noticeable decrease in added resistance equal to 6.4% for the same trim as well. Generally speaking, it appears that trim conditions with bow up to certain measure are favorable for added resistance in regular waves despite the contribution of the projected area of flat bottom in the longitudinal direction, which contributes to the added resistance in waves.

**Figure 11.** Influence of trim on added resistance coefficient in regular waves.

¶

The results of sensitivity analysis regarding trim and its influence on mean value of added resistance in irregular waves for six sea states can be seen in Table 7. For the Bretschneider sea spectrum, except for SS1, there is no significant influence of negative trim on added resistance, which increases by 2–3% for trim conditions equal to −0.8◦ and −1◦. For higher sea states, a significant increase in added resistance can be observed for trim of 0.8◦ by stern. On the other hand, for trim equal to 0.6◦, the decrease in added resistance is equal to about 2–3%. In the case of the JONSWAP sea spectrum, the decrease of added resistance is also the largest for trim equal to 0.6◦. In irregular waves, slight trim by the stern seems to be most favorable for added resistance, even though it does not

contribute significantly to the decrease in added resistance. The largest change in added resistance calculated by means of both wave energy spectra can be observed for SS1 owing to the diffraction component. It should be noted that some of the obtained numerical results are within the calculated numerical uncertainty.


**Table 7.** Influence of trim on added resistance in regular and irregular waves.

The effect of ship characteristics on added resistance in regular head waves depends on the frequency of the incoming wave. By increasing the prismatic coefficient to the highest value considered within this study, added resistance in short waves and the peak value of added resistance increase, while in moderate and long waves, added resistance decreases. In very long waves, the change in the prismatic coefficient as well as the change in the position of the longitudinal centre of buoyancy have a negligible effect. The greatest effect of variation of the *LCB* position is observed at the peak value of the added resistance. As the speed decreases, the added resistance decreases for all analyzed λ/*L* values. Finally, as the pitch gyration radius increases, the value of the added resistance increases as well, more significantly at large and moderate wavelengths. Depending on the frequency range covered by spectral energy, for certain sea states, the influence of the pitch gyration radius will be negligible, while for the other sea states with wave energy distributed in the area of lower frequencies, this influence will be significant. A similar trend can be observed in the case of other varied characteristics. It can be concluded that the results are highly dependent on the incoming wave frequencies and characteristics of sea states, as well on the applied theoretical wave energy spectrum. The distribution of wave energy in the frequency range is dependent on the wave period and the intensity of change in added resistance by varying different ship characteristics on the wave height.

#### **4. Conclusions**

Hydrodynamic calculations of added resistance of KCS in regular waves were performed in the frequency range utilizing the 3D panel method based on the potential flow theory. The obtained numerical results were corrected for the diffraction component of added resistance in short waves. The numerical results were validated against the experimental data available in the literature and satisfactory agreement was achieved. Numerical uncertainty was evaluated for those results in regular

waves, whereas monotonic convergence was achieved and for the mean value of added resistance in irregular waves for certain sea states. At high wave frequencies, there is a potential problem of the appearance of irregular frequencies when calculating second-order forces and the results are unreliable in this area. In that frequency range, the results depend mostly on the applied correction based on the waterline shape, speed, draft, and wave characteristics. Within this research, the effect of ship characteristics on added resistance in regular waves and for certain sea states was investigated. By increasing the prismatic coefficient, while keeping the main dimensions and midship coefficient constant, the added resistance in long and moderate regular waves decreases, except for the peak value, owing to the smaller amplitudes of ship motions. In short waves, the change in added resistance due to the change in the prismatic coefficient is more pronounced. Namely, with the increase in the prismatic coefficient, added resistance, mostly caused by the diffraction of the waves, increases. Spectral analysis based on the corrected results showed that, despite the decrease in the motion amplitudes, the mean value of added resistance in irregular waves for lower sea states increases with the increase in the prismatic coefficient. On the basis of the results obtained by shifting the *LCB* position, the largest effect on added resistance in regular waves was observed for the peak value. Shifting *LCB* towards the bow increases the distance between the longitudinal centre of buoyancy and longitudinal centre of floatation, reducing the amplitudes of relative motions. The mean value of the added resistance in higher sea states decreases by shifting *LCB* towards the stern. As already known, by decreasing the ship speed, a significant reduction in added resistance can be achieved. A significant influence on added resistance was observed for the pitch gyration radius as well, especially for moderate wavelengths. Additionally, sensitivity analysis of trim was performed considering that optimum trim can lead to higher ship efficiency, not only in calm water, but in waves as well. It was observed that trim by stern up to 0.6◦ reduces added resistance for actual sea states by 2–3%.

As the effect of different ship characteristics depends on the distribution of wave energy in the frequency range, the obtained results provide a valuable insight into the effect of the variation of ship characteristics on added resistance for a particular sea state. When designing or optimizing a ship and its propulsion system for calm water conditions, it is important to evaluate the effect of the variation of ship characteristics on added resistance in waves as well, as it is one the main causes of an increase in required power in service. Ship characteristics optimized for calm water performance are not necessarily optimal for real sailing conditions. Considering that a ship in service encounters waves from various directions, it would be beneficial to investigate the effect of variations in ship characteristics for different headings. This will form a part of future research.

**Author Contributions:** Conceptualization, I.M. and N.D.; methodology, I.M., N.D., A.F., and I.G.; software, I.M.; validation, I.M.; formal analysis, I.M.; investigation, I.M., N.D., A.F., and I.G.; resources, I.M.; writing—original draft preparation, I.M., N.D., A.F., and I.G.; writing—review and editing, I.M., N.D., A.F., and I.G.; visualization, I.M; supervision, N.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
