3.1 Routes generation

The presented weather-routing tool has been developed in the MATLAB environment employing functions and toolboxes, such as the mapping toolbox which provides functions for analysing geographic data and creating map displays. In Zaraphonitis et al. (Reference Zaraphonitis, Kytariolou, Dafermos, Gatchell, Östman and Papanikolaou2021), a brief analysis of the developed tool has been presented. The vector format is used for the loaded geographic data, which form sequentially ordered pairs of latitude and longitude coordinates. The vector map can be displayed using a suitable map projection. For navigation purposes, the Mercator projection is the most suitable one, since on this display the loxodromic paths are presented by straight lines on the created map.

By having defined the map environment, an algorithm has been developed generating random routes and displaying them on the map. First, the starting and the ending point of a route are provided, which correspond to the coordinates of the departure and the arrival ports, respectively. Thereafter, there are two options for the generation of the routes. In the first one and based on the n waypoints, a developed algorithm randomly locates them within a predefined envelope. In the second one, the location of all or some of the waypoints is specified. In both options, the group of all these points creates n + 1 legs, where each one can be separately analysed. For instance, it is possible to divide each leg of the whole route into several equidistant points with respective coordinates. These coordinates are unique and based on the available functions, information for each point, such as whether it lies in the sea or at land, the sea depth at that location and its distance from the nearest coast, can be obtained.

With the developed algorithm, a targeted number of different routes are generated, as shown in Figure 2, having as input the starting and the ending points and the random position of the n waypoints. In addition, constraints can be defined such as not allowing the operation within a specific distance from the coast or specifying a minimum depth along the route. Any route that does not satisfy the constraints is characterised as unfeasible and it is not considered further in the process.

Figure 2. Generation of 10 random routes

3.2 Weather and seakeeping data

3.2.1 Weather data

Weather data, and specifically wave, wind and current data, can be acquired at different temporal and spatial resolution. Temporal resolution refers to the time needed to update the information for the same location. During this time step, the weather conditions are assumed to remain constant, which means that having a low temporal resolution can result in poor estimations considering that weather conditions are time sensitive. On the other hand, a high temporal resolution can lead to a time-consuming process with not a considerable change in the quality of the results. Similarly, spatial resolution refers to the distance between every two locations of the grid where weather data are provided.

In case the wave, wind and current data are not available at the same spatial or temporal resolution, the tool sets as the main grid for the optimisation problem the wave's spatial resolution. The intersection points of the main grid as the ship crosses each ‘grid square’ are identified for the examined route, as shown in Figure 3. For these points, shown as red stars, the wave values are calculated through linear interpolation based on the coordinates of these points and their distance from the four picks forming each ‘grid square’ (black dots). Similarly, wind and current values at the intersection points are calculated through linear interpolation based on the nearest grid square formed by their spatial resolution. Finally, between every two intersection points weather conditions are assumed to remain constant, which is an adequate assumption when a dense grid is used.

Figure 3. Main grid formed from wave data resolution equals to 0 ⋅ 5° × 0 ⋅ 5° superposed to 1° × 1° wind resolution. The yellow line represents a random ship route intersecting the example grid section and the red stars represent the intersection points of the random route with the main grid where the desirable values are calculated through linear interpolation

3.3 Fuel oil consumption calculation

3.3.2 Wind resistance

The wind acts as an external force on the vessel, resulting in the wind resistance component. There are several calculation methods, which are based either on semi-empirical formulas, experimental results of models in tunnel tests or CFD results. This study considers only the longitudinal component of the wind resistance. Nevertheless, it is worth mentioning that the transverse component could impose a drift motion that can also increase the resistance or affect ship's manoeuvrability, which, however, are not considered in this study. Figure 4 describes the calculation process.

Figure 4. Wind resistance flowchart

3.3.4 Propeller and main engine modelling

A ship's total resistance $({{R_{\textrm{tot}}}} )$ can be derived by the several components according to the following:

(1)\begin{equation}{R_{\textrm{tot}}} = {R_{\textrm{calm}}} + {R_{\textrm{wave}}} + {R_{\textrm{wind}}}\end{equation}

where ${R_{\textrm{calm}}}$, ${R_{\textrm{wave}}}$ and ${R_{\textrm{wind}}}\;$ refer to calm water, added wave and wind components of the resistance, respectively. Having determined the ship's total resistance, the brake power can be calculated (e.g. Wang et al., Reference Wang, Li, Huang, Ma, Jiang, Yuan, Mwero, Negenborn, Sun and Yan2020):

(2)\begin{equation}{P_\textrm{B}} = \frac{{{R_{\textrm{tot}}}{V_\textrm{S}}}}{{k{\eta _\textrm{s}}{\eta _0}{\eta _\textrm{H}}{\eta _\textrm{R}}}}\end{equation}

where ${R_{\textrm{tot}}}$ is the total resistance, ${V_S}$ is the ship's speed, k corresponds to the number of the propellers, ${\eta _s}$ corresponds to the shaft efficiency, ${\eta _0}$ refers to the propeller's open water efficiency which is equal to (k _TJ)/(k _Q2π), ${\eta _H}$ refers to the hull efficiency equals to (1–t)/(1–w), J is the advance coefficient, and t and w are the trust deduction and wake coefficient, respectively, while ${\eta _R}$ corresponds to the propeller's total rotational coefficient.

When a propeller's characteristics (diameter, pitch ratio, number of blades and the ratio of the expanded blade area) are known, and thrust and torque coefficient curves are provided, then the quantities ${k_\textrm{T}}$, k _Q, $J$, $n$ are possible to be determined from a propeller's open water diagram, which corresponds to calm water conditions. Firstly, the following quantity is calculated:

(3)\begin{equation}\frac{{{k_T}}}{{{J^2}}} = \frac{{T/({\rho {n^2}{D^4}} )}}{{{{({V/nD} )}^2}}} = \frac{T}{{\rho {V^2}{D^2}}} = CC\end{equation}

where $\rho$ is the water density, n is the propeller's revolutions, D is the propeller's diameter and T is the thrust. From the intersection of the curve ${k_\textrm{T}} = CC{J^2}$ with the curve (${k_\textrm{T}} - J$) of the open water diagram of the propeller, the values of $J$, ${k_\textrm{T}}$, k _Q and ${\eta _0}$ can be determined. Then,

(4)\begin{equation}n = \frac{{{\textrm{V}_{ad}}}}{{\textrm{JD}}}\end{equation}

where ${V_{\textrm{ad}}} = {V_s}({1 - w} )$ is the advance speed of the propeller when operating in the ship's wake.

The fuel oil consumption will be calculated based on the main engine's model by using the operating point (required brake power from Equation (2) and main's engine revolution from Equation (4) when no gear box is used) and the corresponding specific fuel oil consumption (SFOC). Using the manufacturer's manual,Footnote ¹ the engine's SFOC map is created following the instructions presented in pages 2 ⋅ 09–2 ⋅ 11. Specifically, the SFOC for an arbitrary load within the loading diagram is needed. Knowing the SFOC value for the nominal MCR L1 rating, reduction rates are provided in the manual for the propeller curve and the constant speed curve for a range of loads. Then, for the operating points that lie between these curves, interpolation is carried out, while for the others, extrapolation. Having calculated SFOC values for a grid of (P,n) pairs, the SFOC contours shown in Figure 5 along with the loading diagram can be derived.

Figure 5. Main engine's SFOC map

In our case, when the point of operation exceeds the main engine's diagram,Footnote ² the route is defined as unfeasible. The process is repeated and when, finally, the ship's speed is attainable, the SFOC value for the specific operating conditions is calculated by the engine's SFOC diagram. Note that when a shaft generator or a shaft power in system are used, respective corrections in the fuel oil consumption are needed. Furthermore, this analysis assumes no engine degradation, while propulsion characteristics are not changed under the effect of weather conditions (e.g. due to ship motions).

3.4 Speed and route management

As mentioned, wave-added resistance can lead to speed loss for the same engine load. Apart from the adverse weather conditions, the captain can also decide to reduce the speed to avoid unsafe effects, such as slamming or propeller immersion. On the other hand, the estimated time of arrival (ETA) is a major factor leading captains to speed up or slow down during the voyage to reach the destination port at the arranged time and avoid port traffic. Moreover, storm avoidance through speeding up or slowing down is also a well-known method for weather routing applications.

In the current framework, the speed management is carried out in the following ways:

• • The ship maintains the commanded speed for the entire route or for a route segment by examining whether it can be achieved under the encountered weather conditions.

• • For a specific ETA, the algorithm can optimise not only the selected route but also the speed between waypoints or any other selected points, aiming at avoiding storms or other dangerous phenomena and reaching the destination port on time.

• • The virtual arrival concept is also considered. Information about port delays can be sent on board during the voyage and an updated speed optimisation can be derived by taking advantage of the voluntary speed loss concept.

It is obvious from these scenarios that route and speed optimisation are interwoven issues that must be handled as one to succeed significant fuels savings.

4.1 Weather data

All weather data are obtained from Copernicus Marine Environment Monitoring Service (CMEMS). Wave data sets (including significant wave height and its direction as well as the wave peak period) are available at 0 ⋅ 2° × 0 ⋅ 2° spatial resolution and for 3-hour temporal resolution. In addition, wind data sets (including true wind speed and its direction) are available at 0 ⋅ 25° × 0 ⋅ 25° spatial resolution and for 6-hour temporal resolution, while wave current data sets (including current speed and its direction) at 0 ⋅ 25° × 0 ⋅ 25° spatial resolution and 12-hour temporal resolution.

4.2 Calm water resistance

Calm water resistance was determined by CFD calculations using FreSCo+ (Hafermann, Reference Hafermann2007). Note that the CFD calculations were significantly close to the ship's model tests, while being utilised to cover a speed range lower than the available one of the model tests (see also Zaraphonitis et al., Reference Zaraphonitis, Kytariolou, Dafermos, Gatchell, Östman and Papanikolaou2021).

4.3 Wind resistance

The wind speed experienced by a moving vessel at sea is called apparent wind speed and differs from the true wind speed that the weather forecasts provide. So, it is necessary to calculate the apparent wind speed ${V_\textrm{a}}$ and direction a, where according to Figure 6, V _s is the ship's speed, V_t is the true wind speed, and $\varTheta$ is the angle formed between the two vectors mentioned, Y = 90–$\varTheta$, and:

(5)\begin{gather}{V_a} = \sqrt {{{({{V_t}\cos Y} )}^2} + {{({{V_t}\sin Y + {V_s}} )}^2}} \end{gather}

(6)\begin{gather}a = {\tan ^{ - 1}}\left( {\frac{{{V_t}\cos Y}}{{{V_t}\sin Y + {V_s}}}} \right)\end{gather}

Figure 6. Apparent wind calculation

Apparent wind speed and its direction are two major parameters to estimate wind resistance. In the presented cases, the wind resistance was determined by using Blendermann's method and coefficients (Blendermann, Reference Blendermann1994), which uses a semi-empirical loading function based on wind tunnel tests. Wind resistance coefficients were also compared (Figure 7), with results obtained using Fujiwara's regression formula (ITTC, 2017) as well as the reference values from the sample containership appeared in Annex F.3 of ITTC (2017).

Figure 7. Comparison of wind resistance coefficients

4.4 Added wave resistance

To calculate the added wave resistance in an irregular sea, the spectral analysis technique (Lloyd, Reference Lloyd1998) was employed as presented next, while the JONSWAP wave spectrum was used. Given that the significant wave height ${H_\textrm{s}}$ and the modal wave period ${T_0}$ describe the sea state, the JONSWAP spectrum is defined as (e.g. Lloyd, Reference Lloyd1998):

(7)\begin{equation}{S_{J\zeta }}(\omega )= 0.658\; C\; {S_{B\zeta }}(\omega )\end{equation}

where $C = {\gamma ^J}$, $J = \textrm{exp}\left[ {\frac{{ - 1}}{{2{\sigma^2}}}{{\left( {\frac{{\omega {T_0}}}{{2\pi }} - 1} \right)}^2}} \right]$, while: $\sigma = 0.07\,\textrm{for}\,\omega < \frac{{2\pi }}{{{{\rm T}_0}}}$ and $\sigma = 0.09\,\textrm{for}\,\omega > \frac{{2\pi }}{{{{\rm T}_0}}}$, ${S_B}\zeta (\omega )= \frac{A}{{{\omega ^5}}}exp \left( {\frac{{ - B}}{{{\omega^4}}}} \right)$, while: $A = 487\frac{{H_s^2}}{{T_0^4}}$ and $B = \frac{{1949}}{{T_0^4}}$.

For the current application the value of the peak-enhancement parameter γ is equal to 3.3, which is the typical value when the JONSWAP spectrum is used, however it shall be noticed that for fully developed and open seas a more refined estimation of this parameter could be considered (e.g. Ochi, Reference Ochi1998).

Added wave resistance depends on the wave frequency, and thus for a vessel advancing through sea waves, the encountered wave frequency is needed:

(8)\begin{equation}{\omega _e} = \omega - \frac{{{\omega ^2}{V_s}cos\mu }}{g}\end{equation}

where $\omega$ is the wave frequency, $\mu$ is the angle defined by the ship's heading and wave direction (which is defined as the relative wave heading), and $\textrm{g}$ = 9 ⋅ 81 m/sec².

A frequency interval $\delta \omega$ around the frequency $\omega$ of the wave energy spectrum is transformed as:

(9)\begin{equation}\delta {\omega _e} = \left( {1 - \frac{{2\omega {V_s}}}{g}cos\mu } \right)\delta \omega \; \; \end{equation}

The area below the spectrum within the interval $\mathrm{\delta \omega }$ is relevant to the energy carried by that frequency range. This energy is not changed by the spectrum transformation and, as a result, the area bounded by the frequency range $\delta \omega$ in the wave energy spectrum will be equal to the area bounded by the frequency range $\delta {\omega _\textrm{e}}$ in the encounter wave spectrum. If $\delta \omega$ and $\delta {\omega _\textrm{e}}$ are infinitesimal, the encounter wave spectrum ordinates are:

(10)\begin{equation}{S_\zeta }({{\omega_e}} )= {S_\zeta }(\omega )\frac{{d\omega }}{{d{\omega _e}}} = {S_\zeta }(\omega )\frac{g}{{g - 2\omega {V_s}\cos \mu }}\end{equation}

The total areas under the two spectra are equal as the total wave energy as well as the significant wave height are not affected by this conversion. Using spectral analysis, the amplitude of a sine wave representing a small range of frequencies $\delta {\omega _\textrm{e}}$ is

(11)\begin{equation}{\zeta _0} = \sqrt {2{S_\zeta }({{\omega_e}} )\delta {\omega _e}} \end{equation}

and the single sine wave results in an added resistance equal to

(12)\begin{equation}2{C_{aw}}{S_\zeta }({{\omega_e}} )\delta {\omega _e} = 2{C_{aw}}{S_\zeta }(\omega )\mathrm{\delta \omega }\end{equation}

where ${C_{\textrm{aw}}} = {R_{\textrm{aw}}}/{\mathrm{\zeta }_0}^2$ is the added resistance response function. Thus, when considering all the wave components, the total added resistance is (Lloyd, Reference Lloyd1998)

(13)\begin{equation}{\bar{R}_{aw}} = 2\int_0^\infty {{C_{aw}}{S_\zeta }({{\omega_{}}} )d\omega } \end{equation}

The added resistance component that corresponds to Equation (12) is pre-calculated by NewDrift+ software, for a range of different wave frequencies, headings and speeds. NEWDRIFT+ is NTUA's in-house software for seakeeping analysis based on a 3D linear panel method and can calculate all ship's responses and added resistance, assuming that the ship is under the effect of regular waves (Liu et al., Reference Liu, Papanikolaou and Zaraphonitis2011).

For each one of the discrete frequencies of the generated wave spectrum, NewDrift+ data are interpolated first regarding the wave frequency, then the heading and finally regarding the ship's speed to derive the added resistance (Zaraphonitis et al., Reference Zaraphonitis, Kytariolou, Dafermos, Gatchell, Östman and Papanikolaou2021) Finally, this methodology is used to calculate the total added resistance in a sea state (Figure 8).

Figure 8. Added wave resistance

4.5 Seakeeping data

The results presented herein are based on the use of the NEWDRIFT+ software for the evaluation of the ship's responses in waves. To represent the wetted surface, the hull is discretised by using trilateral or quadrilateral panels, as shown in Figure 9.

Figure 9. Containership's panelisation

Following Lloyd (Reference Lloyd1998) for the irregular sea state, similarly to the methodology for added wave resistance, the motion energy spectra must be calculated to estimate statistical quantities of the ship motion responses. The resulting motion amplitudes and phases depend on the vessel's speed, the waves’ heading and the encounter frequency. Within linear seakeeping theory and considering the ship as a rigid body, ship motion amplitudes are proportional to the wave amplitude and can be expressed in the form of nondimensional transfer functions. Transfer functions are defined as ship motion amplitudes divided by the wave amplitude (ζ₀) for linear motions and by the wave slope amplitude (k ζ₀) for angular motions, where k is the wave number.

Motion energy spectra are obtained by multiplying the encounter wave spectrum with the squared motion transfer functions of the corresponding encounter frequency. The motion energy spectrum ordinates for each encounter frequency ${S_\mathrm{\zeta }}({{\omega_\textrm{e}}} )$ is calculated by

(14)\begin{equation}{S_{{x_i}}}({{\omega_e}} )= {S_\zeta }({{\omega_e}} ){\left( {\frac{{{x_{i0}}}}{{{\zeta_0}}}} \right)^2}\end{equation}

where i = 1,2,3 correspond to surge, sway and heave motion, respectively. For roll, pitch and yaw motions (i = 4,5,6), the wave slope spectrum must be calculated for the nth sine wave component:

(15)\begin{equation}{S_a}({{\omega_n}} )= \frac{{\omega _n^4}}{{{g^2}}}{S_\zeta }({{\omega_n}} )\end{equation}

Then, the motion energy for roll, pitch and yaw is

(16)\begin{equation}{S_{{x_i}}}({{\omega_e}} )= {S_a}({{\omega_e}} ){\left( {\frac{{{x_{i0}}}}{{k{\zeta_0}}}} \right)^2}\end{equation}

The respective variance and route mean square (rms) values for each motion are calculated according to

(17)\begin{gather}{m_0} = \int_0^\infty {{S_{{x_i}}}({{\omega_e}} )d{\omega _e}} \end{gather}

(18)\begin{gather}{\sigma _0} = \sqrt {{m_0}} \end{gather}

The significant value of the motion energy spectrum is equal to four times the rms motion, since the wave spectra as well as the motion energy spectra are considered narrow-banded spectra. For some critical locations onboard, such as the bridge and the bow, the information concerning the motion displacement, velocity and acceleration could be crucial. The motions on any point on the ship can be calculated based on the motions of the centre of gravity as calculated previously. The angular motions (roll, pitch and yaw) are constant on every point on the vessel, whereas local linear motions (surge, sway and heave) depend on the location on the ship. For example, Figure 10 shows the lateral and vertical accelerations for two critical locations onboard for several wave periods and for a specific value of the significant wave height and relative wave heading $\mu$.

Figure 10. Lateral and vertical accelerations at bridge and at the bow for ship's speed 16 kn, significant wave height 4 m and relative wave heading equal to μ = 150 deg. (Head seas correspond to 180 deg.)

4.6 Optimisation cases

The optimisation is performed through the optimisation toolbox available in MATLAB using a genetic algorithm (see for example Matlab, 2022). Specifically, 100 generations with a population of 100 are chosen for the optimisation processes presented next. The search space is predefined to save computational time, but the optimisation procedure is not based on a fixed-grid technique, as happens in other optimisation approaches (e.g. path finding algorithms, conventional dynamic programming, etc.). Optimisation is also performed for optimal heading (Case A) and both for optimal heading and for nonconstant main engine's and propeller's revolutions (Case C), to approach the route planning problem under realistic conditions. Moreover, the speed optimisation in Case C is combined with a supposed required ETA.

In Case A, the speed is set constant to 16 kn and no safety criteria are considered. Figure 11(a) shows the genetic algorithm evolution during 100 generation of 100 population each.Footnote ³. It is obvious from the figure that the algorithm converges to the optimal result during almost the 10th generation. On the other hand, in Case B, the speed is again kept constant to 16 kn, but the next safety criteria are added:

• • The threshold value for the significant heave is set to 3 ⋅ 5 m.

• • The respective value for the significant roll angle is set to 5°.

These values were chosen for the demonstration of the optimisation study, however operability criteria based on the ability for carrying out several activities onboard either for the crew or passengers can be implemented (e.g. Perera et al., Reference Perera, Rodrigues, Pascoal, Soares, Rizzuto and Guedes Soares2012).

Figure 11. Algorithm's evolution for all cases, where the red circle indicates the optimum ones

The evolution of the algorithm is shown in Figure 11(b). Note also that the number of population and the number of generations are the same as in Case A. As expected, the number of unfeasible routes is higher, and the algorithm converges to the optimal result later in comparison with Case A. Finally, in Case C five new parameters that represent the speed for each route's leg are added to the optimisation process. A maximum travel time of 152 ⋅ 5 h is also added as a constraint. The evolution of the algorithm is shown in Figure 11(c).

A common approach is to compare the results of the optimal routes with the respective ones of the orthodromic and the loxodromic paths, while considering for both of them a constant speed of 16 kn. In the following figures (Figures 12–16) characteristics of these routes are shown in comparison with the optimal ones from each case. The points listed in x-axes represent the intersection points with the main grid, as described in Section 3.2.1.

Figure 12. Significant wave height encountered for the orthodrome, the loxodrome, the optimal path with no safety criteria (Case A) and the optimal path with safety criteria (Case B)

Figure 13. Relative wave heading for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 14. Calm water resistance, accounting the effect of currents, for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 15. Added resistance due to wind for orthodrome, loxodrome and optimal routes (Case A & Case B)

Figure 16. Wave added resistance for orthodrome, loxodrome and optimal routes (Case A & Case B)

As derived by the optimisation procedure, to save fuel, it is preferable to select a route that either passes from regions with low values of the significant wave height (Figure 12) or from regions where the ship encounters following seas (Figure 13). The optimal route in Case B in order to satisfy the set constraints completely avoids head seas, but the avoidance of the bad weather leads to an extended journey with higher fuel oil consumption. Similarly, Figure 14 presents the calm water resistance for the same four paths. Again, the optimal route for Case B follows the favourable path concerning the encountered sea currents, whereas the orthodrome and the loxodrome encounters unfavourable sea currents that increase the calm water resistance and, consequently, the total resistance. Figures 15 and 16 present the added wind and wave resistance, respectively. Obviously, even though the contribution of the wind resistance component to the total resistance is much lower compared with those of calm water and wave, the variations derived by the selections of the optimal route when considering safety criteria reveal that it is a significant element to be considered. Another issue to be discussed is that although the sailing sea state is characterised by heavy waves, the wind resistance is low, since the analysis considers the combined effect of wind and swell waves. As expected, both the optimal routes avoid the passage from regions where the wave-added resistance is high to reduce the required fuel. In addition, Figure 17 shows the resultant optimised speeds at each leg for the optimal path found by the genetic algorithm for Case C. The algorithm suggests reducing the speed in regions where the significant wave height is high or when encountering head seas. Finally, Figure 18 presents the weather evolution each day of travel.

Figure 17. Optimized speed at each leg

Figure 18. Significant wave height evolution within the area of the optimal route (Case A) for each day of the voyage

4.7 Summary of optimisation results

Table 2 summarises the results for all the cases examined, where the maximum of the significant vertical acceleration values, as derived considering each route leg, at the bridge and at the bow have also been included. The optimal route with no safety criteria leads to 4 ⋅ 9% fuel savings compared with the orthodromic path, which is the shortest path and reaches the port of destination almost 2 h earlier than the optimal one. Moreover, when the safety constraints are considered, the optimal path found by the genetic algorithm requires 10% more fuel and 22 h more time than the orthodrome at the same constant speed (16 kn). When speed is optimised (Case C), while the total time duration is set to 152 ⋅ 5 h as a constraint, the optimal path ends in a slight decrease of fuel oil consumption compared with the optimal with constant speed and no safety criteria. All routes are shown in Figure 19.

Figure 19. The derived optimized routes, where the red line = loxodrome, red curve = orthodrome, green polyline = optimal with safety criteria, purple polyline = optimal with optimized speed, black polyline = optimal with no safety criteria

Table 2. Results from all optimisation cases