On the Calculation of Propulsive Characteristics of a Bulk-Carrier Moving in Head Seas

: The present work describes a simpliﬁed Computational Fluid Dynamics (CFD) approach in order to calculate the propulsive performance of a ship moving at steady forward speed in head seas. The proposed method combines experimental data concerning the added resistance at model scale with full scale Reynolds Averages Navier–Stokes (RANS) computations, using an in-house solver. In order to simulate the propeller performance, the actuator disk concept is employed. The propeller thrust is calculated in the time domain, assuming that the total resistance of the ship is the sum of the still water resistance and the added component derived by the towing tank data. The unsteady RANS equations are solved until self-propulsion is achieved at a given time step. Then, the computed values of both the ﬂow rate through the propeller and the thrust are stored and, after the end of the examined time period, they are processed for calculating the variation of Shaft Horsepower (SHP) and RPM of the ship’s engine. The method is applied for a bulk carrier which has been tested in model scale at the towing tank of the Laboratory for Ship and Marine Hydrodynamics (LSMH) of the National Technical University of Athens (NTUA).


Introduction
The problem of a self-propelled ship is one of the most challenging in Naval Architecture since the Engineer is facing an unsteady Three-dimensional (3D) flow, characterized by high Reynold's number, significant viscous effects, the presence of free surface and the complex interaction between the hull and the propeller. Taking into account incident waves makes the problem even harder. Traditionally, the problem of the self-propelled ship in a seaway was tackled experimentally in model scale by separating the problem of the self-propelled ship in calm water from the problem of added resistance in waves. Thus, the solution could only be approximated. Even when a self-propelled scale model is used, the extrapolation of the results to the full-scale ship is problematic due to the scale effect on viscous resistance [1].
Computational Fluid Dynamics (CFD) methods solving the Reynolds Averaged Navier-Stokes (RANS) equations have become of age and are widely used to predict the resistance and propulsion characteristics of ships [2]. Time-dependent RANS solvers have also been developed and can yield promising results with regards to the 6-degrees-of-freedom (6DOF) motions of a ship as well as the added resistance in waves [3]. Solving the time-dependent self-propulsion problem of a ship in waves has proven to be quite time consuming, thus hybrid methods have been proposed. Such methods calculate the effect on the propulsive characteristics of the ship motions, via semi-empirical methods for the prediction of the effect on the ship wake [4].
In recent years, several methods have been proposed in order to solve the problem of the self-propelled ship in incident waves [3,[5][6][7][8]. In most cases, only regular waves are considered due to complete problem of the self-propelled ship in irregular waves [9]. However, these methods require significant computing power in order to obtain accurate results. Therefore, some simplifications are required in order to achieve reliable results and to reduce the computational cost. In this respect, we present a hybrid method for predicting the propulsion characteristics of a ship in irregular head seas in this work, which combines experimental data with CFD calculations.
Having as a starting point experimental model-scale resistance and sea-keeping measurements, we numerically solve the full-scale propulsion problem, where the thrust requirements on the propeller are based on the extrapolation of the model-scale resistance measurements. In this extrapolation, the scale effects are taken into account. Also, the required thrust is obtained by solving the self-propulsion problem, i.e., by calculating the wake fraction and thrust deduction ratio. The actual flow problem is solved by employing the in-house developed RANS solver in conjunction with an actuator disk model for the propeller. The method is used to calculate the propulsion characteristics of a bulk carrier traveling at a steady forward speed in irregular head waves.

The Numerical Method
The employed method has been developed at the Laboratory for Ship and marine Hydrodynamics (LSMH) of the National Technical University of Athens (NTUA). It solves the RANS equations in a body fitted structured mesh as shown in Figure 1. The three regions (I), (II) and (III) refer to three different grid topologies. Region (I) consists of a C-O grid following a sequence of transverse planes which intersect on the Zc-axis. Region (II) represents an H-O grid which is formed by transverse planes to the longitudinal symmetry x-axis, while region (III) forms a typical H-H grid which is generated when a dry or wet transom stern exists, [2]. The C-O topology is applicable for bulbus bow configurations. The H-O grid is the standard grid in most CFD applications around ship hulls. Finally, the H-H is suggested for transom sterns exhibiting discontinuity when the transom is wetted [10].
In any case, the 2D grids on the transverse planes are generated following the conformal mapping technique as described in [10]. In order to apply as far as possible high grid resolutions, the three blocks are created in a restricted domain underneath the predetermined free surface boundary (FS). The free surface is calculated according to the hybrid potential viscous flow method [2].  In any case, the 2D grids on the transverse planes are generated following the conformal mapping technique as described in [10]. In order to apply as far as possible high grid resolutions, the three blocks are created in a restricted domain underneath the predetermined free surface boundary (FS). The free surface is calculated according to the hybrid potential viscous flow method [2].
The RANS code employs a sequence of orthogonal curvilinear systems bounded on each transverse plane. In such a system the three velocity components are denoted as (u i , u j , u l ) and the u i -momentum equation reads: where (h i , h j , h l ) are the metrics, K ij = (1/h i h j )∂h i /∂x j the curvatures of the orthogonal system, ρ the fluid density, p the static plus the hydrostatic pressure, σ ij the stress tensor and the left-hand side of Equation (1) represents the convective terms: The stress tensor components are associated to the deformation tensor through The effective viscosity in Equation (3) is defined as where ν t stands for the turbulent (or eddy) viscosity which is calculated according to the two-equation k-ω-SST turbulence model. The equation for the kinetic energy of turbulence k reads while the equation for the specific dissipation ω reads According to the k-ω-SST model, the turbulent kinematic viscosity is calculated according to the formula In the above Equations (5)-(7), the constants σ k , σ ω , β* and the functions F 1 , F 2 are defined as per the original paper of Menter [11], while G represents the generation term of the turbulence kinetic energy. Essentially, this model separates the flow-field above the body surface in two regions, which are determined according to the local flow parameters and the normal distance from the surface. The k-ω-SST model has been developed in order to calculate effectively complex flow fields, like those appearing in stern flows past full ship forms.
The momentum (RANS) and the turbulence model equations are solved numerically by applying the finite volume method in a staggered grid arrangement. In the resulting non-linear algebraic equations, the convective terms are approximated according to the second order MUSCL scheme [12] associated with the minmod limiter function [13], while the diffusive and the source terms are approximated using central differences. The solution follows the SIMPLE algorithm, where the RANS equations are solved following a marching scheme, while the pressure (continuity) equation is treated as fully elliptic.
Dirichlet boundary conditions are applied on the external boundary N of Figure 1 for the flow variables. The velocity components and the pressure on this boundary are calculated by the potential flow solution performed around the hull and the real free surface [14]. On the free surface, FS in Figure 1, the kinematic condition is fulfilled (zero normal velocity) while for the pressure and the turbulence quantities Neumann conditions are applied. In order to calculate the flow variables above the solid boundary in a uniform way, the standard wall function method is employed. On the exit plane D, open-type conditions hold [15]. Finally, on the flow symmetry plane, along the ship symmetry axis, the normal velocity component is set equal to zero while the normal derivatives of all other variables vanish (Neumann conditions). The exchange of the flow-variables on the matching planes between different blocks is based on a second order interpolation.
A number of domain sweeps is required to obtain convergence of the momentum and continuity equations. Convergence is achieved when the sum of the absolute residuals of the finite volume equations over all nodes becomes lower than a pre-specified criterion. In the present study, this criterion is set equal to 5 × 10 −4 ρV S 2 A M for the momentum equations and equal to 5 × 10 −4 ρV S A M for the continuity equation, where V S represents the ship's speed and A M is equal to the transverse area covered by the grid at the mid-ship section. The self-propulsion problem is solved by adopting the actuator disk model. To introduce this model, the effective wake is approximated according to the open water numerical experiment as in [16,17]. The propeller thrust is calculated following an iterative procedure so that its final value equals the total resistance of the ship. This procedure is also followed when the free surface problem is solved, since the stern wave is influenced by the propeller operation.
Regarding modeling the propulsion in waves, the fundamental simplification adopted in the described method, is that the ship moves in steady forward speed, under rough head sea conditions without performing any other motion, Under this assumption, the propeller thrust should be equal to the hydrodynamic resistance of the ship at any time, i.e., T P (t) = R T,S (t) (8) where T P denotes thrust and R T,S (t) the resistance. In addition, it is assumed that R T,S (t) can be decomposed in two components as: In the above relation R CW,S (t) represents the hydrodynamic resistance of the body moving in still water and R AD,S (t) the added resistance due to the varying wave loads in sea route.
As an extension of the Froude hypothesis, R AD,S (t) is considered as a wave resistance component depending only on the Froude number and is calculated as where C W,AD (t) is the corresponding resistance coefficient, S is the wetted surface of the ship in still water and U S the ships' speed. Since C W,AD (t) depends on the Froude number, it can be calculated by performing towing tank experiments at model scale. In these experiments, we can measure independently the total resistance of the model R T,CW,M in calm water as well as the history of R T,M (t) at a given sea state, always at steady forward speed.
Then, the coefficient C W,AD (t) can be derived through Therefore, the added resistance in Equation (9) can be calculated at any time step, taking into account the scale difference in the time domain, i.e., where t s symbolizes the time in the ship scale, t m the time in the model scale and λ is the geometrical scale ratio. The reason that the component R CW,S (t) in Equation (9) is a function of time, is that any change in the propeller thrust induces a different thrust deduction ratio and has to be calculated by performing an iterative procedure for each time step. In this procedure, the actuator disk model is based on the adoption of body forces that are exerted on fluid particles passing through the plane of the propeller blades. In the present approximation, only the axial body forces are taken into account, being included as source terms in the axial momentum equation. The radial distribution of the body forces is based on the lifting line method, as described in [18]. The body forces are integrated in the actual control volumes as described in [15,17].
The lifting line code requires two fundamental parameters, i.e., the propeller thrust and the velocities at infinite. Both of them are defined by performing a self-propulsion iterative procedure at each time step. In each internal iterative step, the propeller thrust T P (t) is assumed to be equal to the body resistance of the previous step and changes the stern flow field accordingly. Then, a new total resistance R CW,S (t), Equation (9), is calculated and the procedure stops when for three successive time steps self-propulsion is achieved, i.e., Equation (8) holds. Since the flow rate through the propeller blades changes as the iterations are performed, the required lifting line input data are changed accordingly. In addition, at a certain internal step an open water numerical experiment is performed as described by [16] in order to calculate the effective wake fraction at infinite as a function of the flow rate passing through the propeller disk.
The abovementioned procedure is straight forward and provides the history of T P (t) almost independently of the exact distribution of body forces in each time step. Numerical experiments have shown that only an approximate prediction of the circulation distribution along the blades is required, which more or less remains similar under variable lifting line input conditions. However, when we need to calculate the required engine characteristics, we have to employ the lifting line output as concerned the required torque and rpm. To adapt these data, we assume that the propeller operates in a quasi-steady condition and therefore the steady lifting line approximation produces reasonable estimates for the required SHP and RPM.
In a further simplification of the problem, we can also consider the time averaged values for the resistance components of (9): Then the time-dependent propulsion problem is simplified into a steady-state one and the iterative procedure described in the previous paragraph is employed for a single time step. The overall convergence algorithm is presented in Table 1. At first, for time t = 0, the steady state problem is solved and, then, the solution proceeds to the next time steps to solve the unsteady flow at random waves. In each time step, the flow equations are solved under a constant propeller thrust for NS = 30 iterations and then the thrust is adjusted to be equal to the new ship resistance. This procedure is followed until the thrust remains unchanged for 3 successive K steps. Then, the time changes to t = t+ δt, where δt is specified for the particular case, and the solution stops when t becomes higher than a given tmax value.

Characteristics of the Bulk-Carrier
The bulk-carrier under investigation features a bulbous bow and a transom stern as shown in Figure 2. The main characteristics of the ship are shown in Table 2. In this Table, L BP represents the length between perpendiculars, L OA the length overall, L WL the waterline length at full load condition (FL), B the ship's beam, T the draft at full load condition (FL), ∆ the corresponding displacement and C B the block coefficient. The design speed of the ship is equal to 14.82 knots resulting in a Froude number based on L WL equal to 0.180 and a Reynolds number equal to 1.18 × 10 9 . iterations and then the thrust is adjusted to be equal to the new ship resistance. This procedure is followed until the thrust remains unchanged for 3 successive K steps. Then, the time changes to t = t+ δt, where δt is specified for the particular case, and the solution stops when t becomes higher than a given tmax value. Calculate the body resistance RT,S(t) and set propeller thrust TP = RT,S(t) 11 If TP remains constant for 3 successive K-steps GO TO 2 12 GO TO 4 13 END OF THE PROCEDURE

Characteristics of the Bulk-Carrier
The bulk-carrier under investigation features a bulbous bow and a transom stern as shown in Figure 2. The main characteristics of the ship are shown in Table 2. In this Table,

Model-Scale Experiments
In order to evaluate the resistance and sea-keeping characteristics of the ship, a λ = 1/35 scale model was constructed at LSMH. The main characteristics of the model are shown in Table 3. The model was tested in various displacements, trim and speeds [19,20] but in the present work we are concerned with the full load condition and the design speed. Following the Froude scaling method, where the Froude number Fn = 0.180 is kept equal between the ship and the model, the resulting model velocity is U M = 1.289 m/s as shown in Table 3. Table 3. Characteristics of the model. All experiments were conducted at the towing tank of LSMH-NTUA. The main characteristics of the facility are presented in Table 4. For resistance, propulsion and sea-keeping experiments, the models are towed along the tank via a carriage that houses all measuring equipment. The model is towed via a dynamometer, measuring the drag, also allowing for two degrees of freedom in pitch and heave. The tank is also equipped with a paddle-type wave generator able to create harmonic and irregular waves with a maximum amplitude of 0.20 m and a frequency range 0.3-2.0 Hz. The experimental uncertainty has been calculated for previous resistance measurements at the towing tank of LSMH-NTUA, that utilized the same experimental setup [21]. The uncertainty interval for the resistance measurements has been calculated to be less than 2% for the highest resistance values. The experiments were performed previously and presented in [19,20]. Initially the calm water resistance of the model was measured. Then, the resistance and the motions of the model were measured for several incident waves, both harmonic and irregular. The waves were generated using the calibrated paddle-type wave maker. The model velocity was fixed during each measurement. In the present work, we are concerned with the case of irregular waves, following the Bretschneider spectrum [22]:

Scale
where H S = 5.0 m is the significant wave height and T P = 12.0 s is the spectrum peak period, both at full scale. Additionally, ωis the angular frequency, g is the gravitational acceleration, π= 3.14 and α, β spectrum shape parameters. For stochastic procedures such as irregular waves and the corresponding motions of a ship, around 20 min of recordings are required to obtain an accurate representation of all wave components. In model scale thus, the experiments have to run for 200 s. Since the total length of the tank is 90 m and the carriage needs time to accelerate and decelerate, the time of the model at a speed of 1.289 m/s is 30 s. In order to obtain enough data, the measurements were conducted in a total of eight tests each under sequential segments of a specific wave time series corresponding to the abovementioned spectrum. In Figure 3, the λ = 1/35 scale model of the examined bulk-carrier is presented during resistance measurements in the Towing Tank of LSMH-NTUA. In Figures 4-11, the water elevation η is presented, as measured during experiments #1-#8. Finally, in Figure 12, the power density spectrum that results from the combined water elevation measurements is compared to the theoretical shape of the Bretschneider spectrum [22].
where HS = 5.0 m is the significant wave height and TP = 12.0 s is the spectrum peak period, both at full scale. Additionally, ω is the angular frequency, g is the gravitational acceleration, π = 3.14 and α,β spectrum shape parameters. For stochastic procedures such as irregular waves and the corresponding motions of a ship, around 20 min of recordings are required to obtain an accurate representation of all wave components. In model scale thus, the experiments have to run for 200 s. Since the total length of the tank is 90 m and the carriage needs time to accelerate and decelerate, the time of the model at a speed of 1.289 m/s is 30 s. In order to obtain enough data, the measurements were conducted in a total of eight tests each under sequential segments of a specific wave time series corresponding to the abovementioned spectrum. In Figure 3, the λ = 1/35 scale model of the examined bulk-carrier is presented during resistance measurements in the Towing Tank of LSMH-NTUA. In Figures 4-11, the water elevation η is presented, as measured during experiments #1-#8. Finally, in Figure 12, the power density spectrum that results from the combined water elevation measurements is compared to the theoretical shape of the Bretschneider spectrum [22].    (14) where HS = 5.0 m is the significant wave height and TP = 12.0 s is the spectrum peak period, both at full scale. Additionally, ω is the angular frequency, g is the gravitational acceleration, π = 3.14 and α,β spectrum shape parameters. For stochastic procedures such as irregular waves and the corresponding motions of a ship, around 20 min of recordings are required to obtain an accurate representation of all wave components. In model scale thus, the experiments have to run for 200 s. Since the total length of the tank is 90 m and the carriage needs time to accelerate and decelerate, the time of the model at a speed of 1.289 m/s is 30 s. In order to obtain enough data, the measurements were conducted in a total of eight tests each under sequential segments of a specific wave time series corresponding to the abovementioned spectrum. In Figure 3, the λ = 1/35 scale model of the examined bulk-carrier is presented during resistance measurements in the Towing Tank of LSMH-NTUA. In Figures 4-11, the water elevation η is presented, as measured during experiments #1-#8. Finally, in Figure 12, the power density spectrum that results from the combined water elevation measurements is compared to the theoretical shape of the Bretschneider spectrum [22].                                  In Table 5, the measured data for the calm water resistance R T,CW,M as well as the average value of the total resistance at waves R T,M (t) are presented. Additionally, in Table 5, the calculated values of the time-average added resistance in waves and the corresponding coefficient are presented, following Equations (9) and (11). In Figures 13-20, the total resistance R T,M (t) time histories are presented as measured during the experiments. The average resistance in waves, R T,M , and the calm water resistance, R T,CW,M , are also plotted.
In Table 5, the measured data for the calm water resistance RT,CW,M as well as the average value of the total resistance at waves T,M R (t) are presented. Additionally, in Table 5, the calculated values of the time-average added resistance in waves and the corresponding coefficient are presented, following Equations (9) and (11). In Figures 13-20    In Table 5, the measured data for the calm water resistance RT,CW,M as well as the average value of the total resistance at waves T,M R (t) are presented. Additionally, in Table 5, the calculated values of the time-average added resistance in waves and the corresponding coefficient are presented, following Equations (9) and (11). In Figures 13-20

Computations at Model Scale
Computations at model scale were, at first, performed in order to validate the numerical method with respect to the calm water resistance experiment in the towing tank. Since the hybrid method has been adopted for the calculation of the free surface, the potential flow solver was, at first, applied. The calculation domain on the free surface extended from X/L = −0.65 to X/L = 2.2, where L stands for the ship's length and X = 0 at the fore perpendicular. Besides, the external boundary was placed at a distance equal to L from the center line. A total number of 46,000 quadrilateral elements were used, of which 35,000 covered the free surface and 11,000 the wetted surface of the ship. During the calculations, the ship was free to sink. Once the potential solution has been accomplished, the viscous flow code has been applied in order to predict the free boundary around the transom stern and the wake region, in model-scale, thus obtaining a free surface corrected for the viscous effects at the stern. The same external boundaries were used in the computations, while the common boundary between the two solutions located at X/L = 0.85. The grid around the body consisted of 9.5M grid nodes, while 670 transverse sections and 180 grid points in the transverse direction were introduced. In Figure 21, the water elevation contour is presented for the ship at model scale, at a Froude number Fn = 0.180. In Figure 22, the free surface around the bow and the stern of the model is plotted.

Computations at Model Scale
Computations at model scale were, at first, performed in order to validate the numerical method with respect to the calm water resistance experiment in the towing tank. Since the hybrid method has been adopted for the calculation of the free surface, the potential flow solver was, at first, applied. The calculation domain on the free surface extended from X/L = −0.65 to X/L = 2.2, where L stands for the ship's length and X = 0 at the fore perpendicular. Besides, the external boundary was placed at a distance equal to L from the center line. A total number of 46,000 quadrilateral elements were used, of which 35,000 covered the free surface and 11,000 the wetted surface of the ship. During the calculations, the ship was free to sink. Once the potential solution has been accomplished, the viscous flow code has been applied in order to predict the free boundary around the transom stern and the wake region, in model-scale, thus obtaining a free surface corrected for the viscous effects at the stern. The same external boundaries were used in the computations, while the common boundary between the two solutions located at X/L = 0.85. The grid around the body consisted of 9.5M grid nodes, while 670 transverse sections and 180 grid points in the transverse direction were introduced. In Figure 21, the water elevation contour is presented for the ship at model scale, at a Froude number Fn = 0.180. In Figure 22, the free surface around the bow and the stern of the model is plotted.

Computations at Model Scale
Computations at model scale were, at first, performed in order to validate the numerical method with respect to the calm water resistance experiment in the towing tank. Since the hybrid method has been adopted for the calculation of the free surface, the potential flow solver was, at first, applied. The calculation domain on the free surface extended from X/L = −0.65 to X/L = 2.2, where L stands for the ship's length and X = 0 at the fore perpendicular. Besides, the external boundary was placed at a distance equal to L from the center line. A total number of 46,000 quadrilateral elements were used, of which 35,000 covered the free surface and 11,000 the wetted surface of the ship. During the calculations, the ship was free to sink. Once the potential solution has been accomplished, the viscous flow code has been applied in order to predict the free boundary around the transom stern and the wake region, in model-scale, thus obtaining a free surface corrected for the viscous effects at the stern. The same external boundaries were used in the computations, while the common boundary between the two solutions located at X/L = 0.85. The grid around the body consisted of 9.5M grid nodes, while 670 transverse sections and 180 grid points in the transverse direction were introduced. In Figure 21, the water elevation contour is presented for the ship at model scale, at a Froude number Fn = 0.180. In Figure 22, the free surface around the bow and the stern of the model is plotted.
The corrected free surface, considered now as a fixed boundary, is introduced in the potential flow code in order to calculate the velocity and the pressure components on the external boundary. Next, the viscous flow solver was employed to perform accurate resistance calculations underneath the free surface. The calculation domain extended up to X/L = 1.5, while the external boundary was placed at a distance of 0.33 L. Grid dependence tests have been conducted for three grids, in order to study their influence on the numerical solution. The grid nodes for each grid are presented in Table 6, where N 1 , N 2 and N 3 denote the number of nodes along the longitudinal, the circumferential and the transverse directions, respectively, whereas the second term stands for the transom grid. The coarse grid comprises about 5.344 M, the medium 9.208 M and the fine 15.799 M grid points and the ratio between the grid resolutions in successive refinements was about 1.2 in each direction. The corrected free surface, considered now as a fixed boundary, is introduced in the potential flow code in order to calculate the velocity and the pressure components on the external boundary. Next, the viscous flow solver was employed to perform accurate resistance calculations underneath the free surface. The calculation domain extended up to X/L = 1.5, while the external boundary was placed at a distance of 0.33 L. Grid dependence tests have been conducted for three grids, in order to study their influence on the numerical solution. The grid nodes for each grid are presented in Table  6, where N1, N2 and N3 denote the number of nodes along the longitudinal, the circumferential and the transverse directions, respectively, whereas the second term stands for the transom grid. The coarse grid comprises about 5.344 M, the medium 9.208 M and the fine 15.799 M grid points and the ratio between the grid resolutions in successive refinements was about 1.2 in each direction.
The calculated total resistance coefficients CF (friction), CP (pressure) and CT (total) with the three grids are compared in Table 7 and they are defined as The corrected free surface, considered now as a fixed boundary, is introduced in the potential flow code in order to calculate the velocity and the pressure components on the external boundary. Next, the viscous flow solver was employed to perform accurate resistance calculations underneath the free surface. The calculation domain extended up to X/L = 1.5, while the external boundary was placed at a distance of 0.33 L. Grid dependence tests have been conducted for three grids, in order to study their influence on the numerical solution. The grid nodes for each grid are presented in Table  6, where N1, N2 and N3 denote the number of nodes along the longitudinal, the circumferential and the transverse directions, respectively, whereas the second term stands for the transom grid. The coarse grid comprises about 5.344 M, the medium 9.208 M and the fine 15.799 M grid points and the ratio between the grid resolutions in successive refinements was about 1.2 in each direction.
The calculated total resistance coefficients CF (friction), CP (pressure) and CT (total) with the three grids are compared in Table 7 and they are defined as  The calculated total resistance coefficients C F (friction), C P (pressure) and C T (total) with the three grids are compared in Table 7 and they are defined as where S its wetted surface underneath the fixed calm water free surface and U S the ship's speed, while the total resistance R T is the sum of the friction component R F and the pressure component R P , i.e., R T = R F + R P . The difference in C T between the medium and the coarse grid is about 0.3% and about 0.67% between the fine and the medium, while convergence is not monotonic. Therefore, it was considered that the fine grid could be reliable for the comparison purposes. In this case, the calculated total resistance was equal to 24.426 Nt, whereas the measured at the towing tank was 24.340 Nt, resulting in a difference of 0.35% which can be considered satisfactory as regards the applicability of the employed method.

Calm Water Computations at Full Scale
Computations at full scale were, at first, conducted in calm water in order to calculate the resistance of the ship which is the reference for the following added wave resistance tests. The non-dimensional free surface was taken from the model scale calculations, assuming that the Reynolds effect on the stern waves does not play a significant role for the purposes of the present work. The Reynolds number at full scale was equal to 1.18 × 10 9 and the same grid densities of Table 6 were examined, changing suitably the distances of the adjacent to the wall nodes to lie in the desirable region regarding the wall function approach. Table 8 compares the three grid densities with respect to the integrated resistance coefficients. Unlike the model scale results, the differences in C T appear higher mainly owing to the coarser distribution of nodes in the normal direction, since many grid points are now concentrated near the body surface were the boundary layer is thinner due to the Reynolds effect. The convergence is monotonic and the difference in C T between the coarse and the medium grid is 0.86%, while between the medium and the fine grid is 0.55%. By applying the Richardson extrapolation as suggested by ASME [23], the uncertainty of C T at the fine grid is estimated lower than 0.03%. Next, self-propulsion calculations were carried out. The propeller body forces were calculated according to the lifting line approximation, as described previously, and 25 iterative steps were needed to obtain convergence of the thrust, so that it becomes equal to the ship resistance. The main propulsive characteristics are presented in Table 9. By comparing the values of the resistance coefficients to those of Table 8 with the fine grid, it is evident that the increase of C T is mainly due to the increase of C P which is about 38%. The thrust deduction factor 1-t is influenced accordingly. The nominal wake fraction 1-w n appears close to the effective wake fraction 1-w e which is derived according to the open water numerical experiment, [16]. Finally, the delivered horsepower (DHP) and the propeller revolutions (RPM) are calculated by the lifting line code where the input data are the propeller thrust and the mean effective wake. In the time-dependent computations when the added resistance is taken into account, it is assumed that the computational grid remains the same as the one examined in calm water, while only the propeller thrust changes according to the varying wave resistance. Among the deviations from the real situation that arise with this simplification, a significant problem is related to the operation of the propeller under a continuously varying wave surface. In order to explore the influence of the relevant approximation on the propulsive performance, a numerical experiment was performed by comparing two different situations: the ship moving in calm water and the double body case where the free surface is taken as a flat plane. Evidently, in the second situation the hydrodynamic resistance will be substantially lower since the wave resistance is not taken into account. Therefore, it should be interesting to calculate the propulsion characteristics of the double body, if this difference is taken as an additional load on the propeller thrust.
The grid density adopted for the double model was equal to that of the finer grid around the actual ship. Computations at full scale without propeller in operation have shown that the resistance of the double hull was about 39% lower than that of the actual ship, mainly due to the difference in the integrated pressure resistance component. Then, by adding the corresponding load on the propeller thrust, self-propulsion computations were performed and the percentage differences between the ship and the double body case, regarding the thrust, the wake fractions, the advance coefficient J, the propeller efficiency η P , DHP and RPM are depicted in Table 10. It is noticeable that all these propulsive characteristics vary below 2%, indicating that the propeller thrust is the dominant parameter which affects the flow development well underneath the free surface. In order to support this argument, the isowakes (i.e., the ratio of the axial velocity component to the ship's speed) on the propeller plane are plotted in Figure 23, corresponding to the computations without (a) and with propeller in operation (b). On each plot, the data on the left-hand side correspond to the double model while on the right-hand side, the data correspond to the case with free surface. Evidently, the differences between the two cases, double hull vs. free surface, are insignificant on the propeller disk, while deviations are observed close to the free surface. Therefore, these trends justify the small differences in Table 10, concerning the nominal and the effective wake fractions. Naturally, this comparison is only indicative since in reality the ship moves in random waves, which may affect more drastically the flow past the propeller.

Computations at Head Sea State
In order to study the effect of the added resistance when the ship is advancing in steady speed under random head-sea conditions, computations in the time domain have been performed following the assumptions presented in Equations (8)-(11) utilizing the resistance time series as measured in the towing tank. Thus, eight numerical experiments were conducted, each lasting 30 s. The time step in the calculations was Δt = 1 s and in each time step self-propulsion was achieved in lower than 10 iterative steps, while 15 domain sweeps were performed per iteration. Naturally, at the end of this procedure per time step, the propeller thrust equals the hydrodynamic resistance of the ship. In Figure 24, corresponding to the numerical experiment #1, the time histories of the thrust coefficient CT = (TP/(0.5ρSUS 2 )) and the added wave resistance CW = (CW,AD(t)) are plotted. Evidently, CT follows the variations of CW but their difference is not constant, since the thrust deduction factor varies in time. It should be noticed here that a cut-off for the negative thrust components is applied, since negative thrust values imply reverse propeller operation which is impossible in real life. Thus far, the results depend on the approach of the actuator disk and represent the "real" effects under the basic assumptions of the proposed method. Next, the quasi-steady approximation is followed in order to estimate the powering demands. In Figure 25, the flow rate through the propeller disk 1-wd and the effective wake fraction 1-we are

Computations at Head Sea State
In order to study the effect of the added resistance when the ship is advancing in steady speed under random head-sea conditions, computations in the time domain have been performed following the assumptions presented in Equations (8)-(11) utilizing the resistance time series as measured in the towing tank. Thus, eight numerical experiments were conducted, each lasting 30 s. The time step in the calculations was ∆t = 1 s and in each time step self-propulsion was achieved in lower than 10 iterative steps, while 15 domain sweeps were performed per iteration. Naturally, at the end of this procedure per time step, the propeller thrust equals the hydrodynamic resistance of the ship. In Figure 24, corresponding to the numerical experiment #1, the time histories of the thrust coefficient C T = (T P /(0.5ρSU S 2 )) and the added wave resistance C W = (C W,AD (t)) are plotted. Evidently, C T follows the variations of C W but their difference is not constant, since the thrust deduction factor varies in time. It should be noticed here that a cut-off for the negative thrust components is applied, since negative thrust values imply reverse propeller operation which is impossible in real life. Thus far, the results depend on the approach of the actuator disk and represent the "real" effects under the basic assumptions of the proposed method.
(a) (b) Figure 23. Axial velocity isowake curves on the propeller plane; on the left is the double body case and on the right the ship in calm water, for (a) the towed ship and (b) the self-propelled ship.

Computations at Head Sea State
In order to study the effect of the added resistance when the ship is advancing in steady speed under random head-sea conditions, computations in the time domain have been performed following the assumptions presented in Equations (8)-(11) utilizing the resistance time series as measured in the towing tank. Thus, eight numerical experiments were conducted, each lasting 30 s. The time step in the calculations was Δt = 1 s and in each time step self-propulsion was achieved in lower than 10 iterative steps, while 15 domain sweeps were performed per iteration. Naturally, at the end of this procedure per time step, the propeller thrust equals the hydrodynamic resistance of the ship. In Figure 24, corresponding to the numerical experiment #1, the time histories of the thrust coefficient CT = (TP/(0.5ρSUS 2 )) and the added wave resistance CW = (CW,AD(t)) are plotted. Evidently, CT follows the variations of CW but their difference is not constant, since the thrust deduction factor varies in time. It should be noticed here that a cut-off for the negative thrust components is applied, since negative thrust values imply reverse propeller operation which is impossible in real life. Thus far, the results depend on the approach of the actuator disk and represent the "real" effects under the basic assumptions of the proposed method. Next, the quasi-steady approximation is followed in order to estimate the powering demands. In Figure 25, the flow rate through the propeller disk 1-wd and the effective wake fraction 1-we are Next, the quasi-steady approximation is followed in order to estimate the powering demands. In Figure 25, the flow rate through the propeller disk 1-w d and the effective wake fraction 1-w e are plotted versus time. The first one is the direct result of the solution of the RANS equations, while the second one is found by an open-water numerical experiment, so that the velocity at infinite produces w d through the disk under the known T P [16]. Apparently, the two curves do not present the same trends, while 1-w e is always lower, as expected, than 1-w d . Besides, 1-w d exhibits unusually high values, which are due to the high values of the added wave resistance and their memory effects. In any case, 1-w e and T P are the crucial parameters being the input data in the lifting line code or in the open water propeller characteristics whenever they are available. In Figures 26 and 27, the isowake contours that are associated with the calculation of the fraction w d are plotted correspondingly for a local maximum and a local minimum of the added wave resistance component. Evidently, the thrust following the resistance history affects drastically the wake field and enhances strongly the axial velocity component when it takes high values. plotted versus time. The first one is the direct result of the solution of the RANS equations, while the second one is found by an open-water numerical experiment, so that the velocity at infinite produces wd through the disk under the known TP [16]. Apparently, the two curves do not present the same trends, while 1-we is always lower, as expected, than 1-wd. Besides, 1-wd exhibits unusually high values, which are due to the high values of the added wave resistance and their memory effects. In any case, 1-we and TP are the crucial parameters being the input data in the lifting line code or in the open water propeller characteristics whenever they are available. In Figures 26 and 27, the isowake contours that are associated with the calculation of the fraction wd are plotted correspondingly for a local maximum and a local minimum of the added wave resistance component. Evidently, the thrust following the resistance history affects drastically the wake field and enhances strongly the axial velocity component when it takes high values.   plotted versus time. The first one is the direct result of the solution of the RANS equations, while the second one is found by an open-water numerical experiment, so that the velocity at infinite produces wd through the disk under the known TP [16]. Apparently, the two curves do not present the same trends, while 1-we is always lower, as expected, than 1-wd. Besides, 1-wd exhibits unusually high values, which are due to the high values of the added wave resistance and their memory effects. In any case, 1-we and TP are the crucial parameters being the input data in the lifting line code or in the open water propeller characteristics whenever they are available. In Figures 26 and 27, the isowake contours that are associated with the calculation of the fraction wd are plotted correspondingly for a local maximum and a local minimum of the added wave resistance component. Evidently, the thrust following the resistance history affects drastically the wake field and enhances strongly the axial velocity component when it takes high values.    The results of the lifting line code concerning the propeller DHP and RPM are plotted on Figures  28-35 for numerical experiments #1-#8. As observed, a quite unrealistic behavior is observed for DHP, exhibiting prohibitive high values at the positive peaks of CT. The same conclusion holds for RPM which presents peaks of about 300% as regards the steady case. It should be taken into account that this behavior represents the propeller demands, whereas it is impossible for the engine to respond, because it has not the necessary margin and, besides, it cannot follow the extremely sudden peaks. Therefore, sustaining a constant speed seems to be an impossible situation, rendering the corresponding towing tank experiments questionable unless approximate estimations are made.  The results of the lifting line code concerning the propeller DHP and RPM are plotted on Figures 28-35 for numerical experiments #1-#8. As observed, a quite unrealistic behavior is observed for DHP, exhibiting prohibitive high values at the positive peaks of C T . The same conclusion holds for RPM which presents peaks of about 300% as regards the steady case. It should be taken into account that this behavior represents the propeller demands, whereas it is impossible for the engine to respond, because it has not the necessary margin and, besides, it cannot follow the extremely sudden peaks. Therefore, sustaining a constant speed seems to be an impossible situation, rendering the corresponding towing tank experiments questionable unless approximate estimations are made. The results of the lifting line code concerning the propeller DHP and RPM are plotted on Figures  28-35 for numerical experiments #1-#8. As observed, a quite unrealistic behavior is observed for DHP, exhibiting prohibitive high values at the positive peaks of CT. The same conclusion holds for RPM which presents peaks of about 300% as regards the steady case. It should be taken into account that this behavior represents the propeller demands, whereas it is impossible for the engine to respond, because it has not the necessary margin and, besides, it cannot follow the extremely sudden peaks. Therefore, sustaining a constant speed seems to be an impossible situation, rendering the corresponding towing tank experiments questionable unless approximate estimations are made.          In order to investigate the effect of the selected time step, the numerical experiment #1 was repeated using a time step of 0.1 s. In Figure 36, the time histories of the thrust coefficient CT are plotted for the Δt = 1 s and the Δt = 0.1 s tests. It is obvious that by using the lower value of Δt some peaks of the CW,AD can be captured that have a direct influence on the total resistance coefficient. This is clearly shown around t = 15 s and t = 25 s. It is noticeable that the existence of such peaks in CT do not affect the rest of the time history. As a direct consequence of this CT behavior, the calculated values for the required DHP and propeller RPM are influenced accordingly, as shown in Figure 37. However, the observed differences due to Δt do not essentially change the basic conclusions In order to investigate the effect of the selected time step, the numerical experiment #1 was repeated using a time step of 0.1 s. In Figure 36, the time histories of the thrust coefficient C T are plotted for the ∆t = 1 s and the ∆t = 0.1 s tests. It is obvious that by using the lower value of ∆t some peaks of the C W,AD can be captured that have a direct influence on the total resistance coefficient. This is clearly shown around t = 15 s and t = 25 s. It is noticeable that the existence of such peaks in C T do not affect the rest of the time history. As a direct consequence of this C T behavior, the calculated values for the required DHP and propeller RPM are influenced accordingly, as shown in Figure 37. However, the observed differences due to ∆t do not essentially change the basic conclusions concerning the prohibitive values of DHP and RPM. In order to investigate the effect of the selected time step, the numerical experiment #1 was repeated using a time step of 0.1 s. In Figure 36, the time histories of the thrust coefficient CT are plotted for the Δt = 1 s and the Δt = 0.1 s tests. It is obvious that by using the lower value of Δt some peaks of the CW,AD can be captured that have a direct influence on the total resistance coefficient. This is clearly shown around t = 15 s and t = 25 s. It is noticeable that the existence of such peaks in CT do not affect the rest of the time history. As a direct consequence of this CT behavior, the calculated values for the required DHP and propeller RPM are influenced accordingly, as shown in Figure 37. However, the observed differences due to Δt do not essentially change the basic conclusions concerning the prohibitive values of DHP and RPM.  Finally, it is interesting to study the self-propulsion of the ship at steady conditions, regarding the mean added wave resistance as an additional load on the propeller thrust. According to the experimental data, the mean added wave resistance component is calculated equal to 1.764 × 10 -3 and the corresponding propulsive characteristics are presented in Table 11. Comparing the results with those of Table 9 for the calm water predictions, we can find an average increase of about 77% in DHP and 17.5% in RPM. These differences are still substantially higher than the usual margins of the ship's engine, implying that there will be an actual loss of speed at the examined sea state.

Conclusions
A simplified numerical method has been employed for calculating the propulsion characteristics of a bulk-carrier in waves. Under the assumption that the ship moves at steady forward speed, the added resistance in time is estimated through the towing tank experiments and the propeller thrust is adjusted accordingly, while propeller action is modeled according to the actuator disk approach. Finally, it is interesting to study the self-propulsion of the ship at steady conditions, regarding the mean added wave resistance as an additional load on the propeller thrust. According to the experimental data, the mean added wave resistance component is calculated equal to 1.764 × 10 -3 and the corresponding propulsive characteristics are presented in Table 11. Comparing the results with those of Table 9 for the calm water predictions, we can find an average increase of about 77% in DHP and 17.5% in RPM. These differences are still substantially higher than the usual margins of the ship's engine, implying that there will be an actual loss of speed at the examined sea state.

Conclusions
A simplified numerical method has been employed for calculating the propulsion characteristics of a bulk-carrier in waves. Under the assumption that the ship moves at steady forward speed, the added resistance in time is estimated through the towing tank experiments and the propeller thrust is adjusted accordingly, while propeller action is modeled according to the actuator disk approach.
At first, computations were performed at model scale in still water in order to validate the method. Then, similar computations at full scale were conducted in order to calculate the ship's resistance which was the basis for determining the added resistance in waves and the required propeller thrust.
Then, unsteady RANS computations were applied at a random sea state characterized by a prescribed power spectrum. In order to study the variation of propulsive parameters (including DHP and RPM) in time, it was assumed that the ship moved steadily and, besides, the propeller performance could be approximated by a quasi-steady procedure. The estimation of these parameters in a certain time period indicated that the propeller demands vary in an unrealistic way as regards the engine response. This behavior is rather due to the abrupt changes of the wave resistance which impose extremely high loads on the propeller. Therefore, the basic assumption of the steady forward speed in waves is subject to criticism.
Another crucial issue, related to the aforementioned study, is that the towing tank experiments in random waves at constant speed appear questionable as regards the self-propulsion behavior of a real ship. However, they are very useful to estimate the mead added resistance and, consequently, to predict the mean propulsion demands in a seaway. Such a numerical experiment was also performed in order to compute approximately the required engine margins for the examined bulk carrier.
Author Contributions: S.P. and G.T. conceived and designed the numerical method; performed the simulations; analyzed the data and co-wrote this paper. All authors have read and agreed to the published version of the manuscript.