CFD-Powered Ship Trim Optimization: Integrating ANN for User-Friendly Software Tool Development

1. Introduction

With the aim to reduce emission of greenhouse gasses, in 2013 the International Maritime Organization (IMO) introduced the Energy Efficiency Design Index (EEDI) for new ships and later in 2023 the Energy Efficiency for Existing Ships Index (EEXI). These new measures initiated changes in the shipbuilding industry. By implementing these regulations, it is anticipated to reduce global CO₂ emissions by 40% by year 2030 and 70% (or at least 50%) by year 2050, compared to the levels in 2008. All ships that fall under MARPOL Annex VI and have over 400 GT must meet the set criteria which is based on nominal ship data. IMO determined the parameters and procedure for calculating the required EEDI and EEXI, as well as the formula for calculating the attained EEDI/EEXI that must be lower than the required values. The required value for EEDI is determined by the formula:

Required EEDI = a(DWT)^−c

(1)

where parameters a and c are based on the regression curve fit [1] of data from ships that were built between 1999 and 2009, and also depend on the type of ship. The required EEDI decreases over the years according to the formula:

Attained EEDI ≤ Required EEDI = (1 – X/100)·reference line value

(2)

where X is the specified reduction factor given in [2] through four phases. Values for EEXI are based on the required EEDI, taken into account with the adequate reduction factor from [3]. The requirements for EEXI are becoming slightly stricter with time, with the plan of EEXI and EEDI becoming equal in 2025. Beside the two mentioned energy efficiency parameters, IMO introduced operational efficiency indicators. In 2009 [4] the Energy Efficiency Operational Indicator (EEOI) was invented and is defined as the ratio of emitted CO₂ per unit of transport work. EEOI never became mandatory but it is a representative value of the ship’s operational energy efficiency level for a consistent period of time. Another indicator came into force in 2023, the Carbon Intensity Index (CII), which is calculated as the ratio of the total mass of emitted CO₂ to the total transport work undertaken in a specific calendar year, considering travelled distance, time spent underway and the amount of consumed fuel over a period of one year. As well as for EEDI and EEXI, attained [5] and required [6] values are also calculated for CII, where the attained CII has to be lower than the required, which decreases over time [7]. Based on the attained CII, the ship belongs to one of five energy efficiency categories (A, B, C, D, E). The boundaries between those categories and the calculation method for each ship type are defined in [8] and they depend on the ship’s DWT. Since CII is a mandatory parameter, all ships of 5000 GT or more are from 2019 obliged to record the relevant information (travelled distance, time spent underway and the amount of consumed fuel over a period of one year) under the Data Collection System (DCS) [9]. In case that a ship achieved class E for at least one year, or class D for three years in a row, it must create a plan for reducing CO₂ emissions and satisfy the requirements for at least class C [10]. In 2013 IMO introduced the Ship Energy Efficiency Management Plan (SEEMP), which consists of three parts: ship management plan to improve energy efficiency, ship fuel oil consumption data collection plan and ship operational carbon intensity plan. Part I applies to all ships that fall under MARPOL Annex VI and have over 400 GT, while Part II and Part III apply for ships that have over 5000 GT and fall under MARPOL Annex VI.

Various measures can be considered in order to meet the energy efficiency requirements developed by the International Maritime Organization and attain reduction in CO₂ emissions. Those measures are divided in two categories: operational methods and technical methods. Operational methods include: improvement in voyage execution [11,12], reduction of auxiliary power consumption [11,13], weather routing [10,12], “just in time” voyage [13], optimum ballast [10,12], optimum cargo distribution [12], energy saving utilities [14], optimum use of rudder and heading control systems [10], optimized hull and propeller maintenance [10], speed optimization [12,15], slow steaming [16,17] and trim optimization [12,18].

These measures can be applied for both existing ships as well as new ships. Technical measures are design related and therefore more favourable for new ships. These measures refer to wind-assisted propulsion, fuel type change [10,16,19], waste heat recovery [12], upgrading and maintenance of propulsion system, hull retrofit (bulb and/or stern modification, installation of energy saving devices [10], etc.

Optimum trim means that the angle of trim for a specific operating condition, regarding displacement and speed, provides minimum resistance which directly implies optimal efficiency level [12]. The resistance of the ship changes depending on the trim, although the displacement and speed stay the same [20]. Beneficial aspect of trim optimization is that neither hull modification nor engine upgrade is needed. The ship is trimmed if the draught at bow differs from the draught at the aft section of the ship. While negative trim indicates that the draught at the bow is greater than the draught at stern, positive trim implies the opposite. Every vessel is optimized for a number of conditions (even keel at full load and design speed, ballast condition, etc.), but the actual operating conditions usually differ from the expected ones [21]. It was confirmed in [22] that trim can affect the total resistance of a ship for various service speeds and the optimum trim for every speed is different. Trim optimization method is intended for minimizing the resistance in calm water and therefore minimizing fuel consumption which can be accomplished with a specially developed program for the ship. The results of [18] indicate the possibility of reducing total fuel consumption during a whole voyage by 1.2% by utilizing calculated optimum trim for whole voyage. 4250 TEU container ship sea trial with use of a trim optimizing program reported a main engine power reduction of 910 kW and energy saving rate of 9.2% [23]. Trim of the ship can be influenced by redistribution of ballast, fuel and/or load between tanks. Parameters that change when the ship is trimmed compared to even keel condition are wetted surface area, length of waterline and submerged hull form at bow as well as at stern [24]. Potential disadvantages of this method could be reduction of visibility, reduced freeboard, emergence of the propeller [25], underkeel clearance, seakeeping and shipped water on deck [26], and those should as well be taken into account while finding the optimum trim for an operating condition.

Calculating ship resistance in calm water by conducting model tests and numerical simulations provides data about resistance for different drafts and trims. As a result, a set of curves with highlighted lowest resistance for a specific draft is obtained, as in [22,27]. While onboard tools have this information at disposal, model tests are still being common and basins are equipped with traditional procedures and up-to-date insights to provide their proper execution, pursuing to measure small power variations within foreseen range of 0 to 4% of the total installed power [28], or 2% to 4% fuel consumption reduction [29]. With the advancement of technology and computers, CFD software calculations provide trim tables with an accuracy that can compete with the results gained from traditional model tests, with even less investments. Other methods that can be used to determine the optimum trim are sea trials and machine learning method. The shortcomings of these methods are that sea trials are fuel and time-consuming for determining the optimum ballast, while for obtaining information about the optimum trim with the machine learning method a lot of data from ship’s past voyages is needed [30]. In this study, the aim of work is focused only on applying CFD method on one RO-RO car carrier in order to optimize the trim for an energy efficient voyage and exploitation.

Application of CFD method for trim optimization for different ship types can be found in various existing studies. For example, in the research by [20] the effect of trim optimization for a container ship like MOERI (KCS) was a total resistance reduction of 2%, similar as the results of trim optimization for a US Navy ship [13]. A CFD investigation of the propulsion performance of a low-speed VLCC tanker at various initial trim angles was conducted by [31] and showed 1.76-2.12% total resistance reduction. [32] and [33] report even greater profit, such as fuel savings of up to 5% [32] and 6% reduction in delivered horse power [33]. Since the trim conditions can vary significantly, so do the results from trim optimization [34] highlights that in Series 60 the total resistance varies between the worst and best trims up to 11%. Trim optimization of bulk carrier in [35] results with total resistance reduction possibility by up to 14%. Savings for RO-RO ships were studied by [12] and indicated a possibility of up to 10.4% reduction in delivered power, or 1.2t fuel per day. Beside the calculation of optimal trim, specialized software that for input parameters quickly provide information about optimal trim have been mentioned [33,36].

Up to now, artificial neural networks (ANN) have found application in predicting fuel consumption, as demonstrated in studies by [37] and [38], which focused on utilizing NOON reports. Furthermore, ANN has been employed for forecasting ship speed, as evidenced by the work of [39]. Moreover, research efforts have extended to utilizing ANN for the joint prediction of ship speed and fuel consumption, leveraging data from sails, particularly in the case of a barquentine, as explored by [40]. The research [41] proposes a real-time hybrid electric ship energy efficiency optimization model considering time-varying environmental factors, aiming to optimize the EEOI under wind and wave conditions while maintaining speed limits, resulting in an average reduction of fuel consumption by 13.4% and real-time EEOI by 15.2%. In [42] is presented a real-time prediction model of ship fuel consumption through BP neural network training-related data, and further used it for ship speed optimization.

In today’s complex technological landscape, the integration of various disciplines is becoming increasingly vital for tackling engineering challenges. By bridging the fields of naval engineering, CFD, ANN and practical application development, this study not only underscores the importance of multidisciplinary association but also showcases its benefits in addressing real-world problems within the maritime industry. Through the synergy of these varied fields, novel solutions can be developed to optimize ship operations or design, enhance fluid flow analysis accuracy, and streamline decision-making processes. Moreover, the utilization of ANN enables the creation of sophisticated mathematical models that can effectively capture complex relationships and patterns within maritime systems, paving the way for more precise simulations and predictive analytics. Furthermore, the development of user-friendly applications facilitates the seamless implementation of these advanced methodologies, empowering stakeholders to leverage cutting-edge insights for improved vessel performance and operational efficiency. Essentially, this methodology not only promotes scientific comprehension but also encourages innovation and pragmatic breakthroughs in the field of maritime engineering, thus helping in a long-term growth of maritime technologies.

2. Methods

The conducted work consists of several phases: 1) 3D modeling; 2) First level verification of the reliability of the digital twin (3D model); 3) Open Water Test (OWT); 4) Second level verification of the reliability of the 3D model; 5) Trim optimization; 6) Determination of a mathematical model for assessing optimal trim, ship speed, engine (brake) power, daily fuel consumption, and propeller speed; 7) Programming an application based on previously conducted analyses. In this study various software have been used: StarCCM+ for CFD analysis, aNETka for mathematical model determining and MATLAB for application development.

In Table 1 are given ship’s (RO-RO car carrier) principal particulars.

Propeller geometry characteristics are presented in Table 2.

2.1. 3D Modeling

3D modeling is the first necessary step before conducting CFD simulations. Creating a 3D representation of the RO-RO vessel involved using existing documentation: hull construction drawings, rudder construction drawing and propeller drawing. This digital twin includes detailed modeling of the hull, rudder, and propeller, with the goal of faithfully capturing the ship’s attributes and performance.

2.2. First Level Verification

First level verification involves comparing hydrostatic data such as displacement (Δ), block coefficient (C_B), Longitudinal center of buoyancy (LCB), and wetted surface (WS) of the 3D model with data given in the Trim & Stability booklet for various drafts. No parameter should deviate more than 1% from (as per authors experience 1% deviation enables enough accuracy) its corresponding value given in the Stability Booklet.

2.3. Open Water Test

The OWT is a method for determining propeller characteristics such as thrust coefficient, torque coefficient, and propeller efficiency. The procedure for determining these coefficients through model tests is defined in [43]. In this case, model testing was not conducted; instead, propeller parameters were determined through CFD simulation, following the procedure described in [44]. According to this document, CFD results of the OWT are considered valid if they do not differ more than 3% from results obtained by model testing within the relevant propeller operating range. It is also emphasized that if model test results of the propeller are not known, an additional OWT with the appropriate B-series propeller model should be conducted. An appropriate B-series propeller is considered one with the same geometric characteristics (diameter, pitch, ratio of characteristic areas) as the original propeller installed on the considered ship. The 3D model of the B-series propeller was obtained using the B-series Propeller Generator available online [45]. OWT results of the B-series propeller should not deviate more than 3% from results obtained with mathematical model for the same geometric parameters. If this condition is met, the same CFD calculation methodology is applied to the original propeller, and results obtained this way are considered valid. A mathematical model can be found in [46].

The simulations are conducted using a steady state numerical flow model, solved in the rotating frame of reference. The numerical method used is called Multiple Reference Frames (MRF), where in one part of the domain the equations are solved in rotating frame of reference, while in the rest of the domain stationary frame of reference is used.

Figure 1 shows the domain of the open water test CFD setup. The outer cylinder shows the far-field boundary of the domain, while the inner, smaller cylinder represents the boundary between the rotational frame of reference and stationary frame of reference. Diameter and length of small cylinder are 4D, while diameter and length of large cylinder are 10D. The propeller is positioned downstream from the shaft. All simulations are conducted in full scale.

The same grid refinement options were used for both propellers grids (B series propeller and Original propeller) are presented in Figure 2.

Turbulence modelling is performed using the k-ω SST model. Inflow turbulence intensity was set to 1% relative to the inflow speed. The inflow speed conditions, and rotation rate of the propeller are selected to closely resemble the conditions expected at self-propulsion conditions.

2.4. Second Level Verification

Second level verification involves conducting CFD simulation for a draft for which there are available results from model tests or sea trials. The difference between results obtained by CFD analysis and model tests/sea trials should not be greater than 5%.

The equations are discretized using collocated Volume of Fluid (VOF) multiphase method implemented within software. It is used for large–scale two–phase (in this case sea water (ρ_w=1.025 t/m³) and air (ρ_A =1.212 kg/m³) at 18⁰C) flows encountered in naval hydrodynamics. A two–phase, incompressible, turbulent and viscous flow model is employed, governed by the continuity and Navier-Stokes equations:

$$\nabla \cdot V=0,$$

(3)

$$\frac{\partial V}{\partial t}\nabla \cdot \left(V\otimes V\right)=-\frac{1}{\rho }\nabla \cdot \left(pI\right)+\frac{1}{\rho }\nabla \cdot T_v+\frac{1}{\rho }{f}_b,$$

(4)

where V stands for the velocity vector field, p is the pressure, T_v is the viscous stress tensor and f_b is the resultant of the body forces.

The VOF multiphase Eulerian approach model implementation in Simcenter STAR-CCM+ belongs to the family of interface-capturing methods that predict the distribution and the movement of the interface of immiscible phases. This modeling approach assumes that the mesh resolution is sufficient to resolve the position and the shape of the interface between the phases.

An important quality of a system of immiscible phases (air and water) is that the fluids always remain separated by a sharp interface. The High-Resolution Interface Capturing (HRIC) scheme is used to mimic the convective transport of immiscible fluid components, resulting in a scheme that is suited for tracking sharp interfaces.

2.4.1. Actuator Disc Propeller Model

The body force propeller method is used to model the effects of a propeller such as thrust and thereby creating propulsion without actually resolving the geometry of the propeller. The method employs a uniform volume force distribution over the cylindrical virtual disk. The volume force varies in radial direction, so the total force can be calculated as volume integral for the whole modeled cylinder. The radial distribution of the force components follows the Goldstein optimum [47] and is given by:

$$f_{bx}=A_x{r}^{*}\sqrt{1-r^{*}},$$

(5)

where r* is the normalized disc radius defined as:

$$r^{*}=\frac{r^{\prime }-r_h^{\prime }}{1-r_h^{\prime }},$$

(6)

where r’=r/R_p and r_h’=R_h/R_p. Volume force f_bx is body force component in axial direction, r is the radial coordinate, R_h is the hub radius and R_p is the propeller tip radius. Coefficient A_x is defined as:

$$A_x=\frac{105}{8\pi }\frac{T}{t_d\left(R_p-R_h\right)\left(3R_h+4R_p\right)},$$

(7)

where T stands for propeller thrust and t_d is the virtual disk thickness. More about this can be found at [48]. The simulation is performed for a certain operating point specified by thrust (T). The advance ratio is calculated by solving the following equation numerically (J—advance coefficient, D—propeller diameter):

$$f\left(J\right)=K_T-{K_T}^{\prime },$$

(8)

where K_T is evaluated from the propeller performance curve and K_T’ is evaluated as:

$${K_T}^{\prime }=\frac{J^2{T}_{Operatingpoint}}{{\rho }_w{V}_A^2{D}^2}$$

(9)

With K_T available, the thrust is calculated for the propeller:

$$T=\frac{K_T{\rho }_w{V}_A^2{D}^2}{J^2}$$

(10)

With T (T_{(Operating point)}) available, the axial body force component can be calculated. Inflow plane is the plane inside actuator disk where volume-averaged velocity and density are computed. Input thrust in this part of the simulation is equal to the total resistance that is evaluated as sum of the pressure resistance and viscous resistance of the hull with rudder included.

The wake fraction (w) is extracted at propeller axis plane in steady resistance simulation.

Thrust deduction factor is calculated as per:

$$t=1-\frac{R_T}{T},$$

(11)

where R_T is total resistance of the ship evaluated in steady resistance simulation.

2.4.2. Steady Resistance Simulations

Total resistance, including viscous and pressure resistance, is obtained by re-running the self-propulsion simulation with disabled actuator disc. This procedure implies the same mesh as it was used earlier but only without an actuator disc.

2.4.3. Boundary Conditions

Fluid domain dimensions are defined as follows (from origin) this:

• Inlet: 2.5 Lpp;
• Outlet: 3 Lpp;
•  Bottom: 1.5 Lpp
•  Top: 1 Lpp;
• Port Side: 2 Lpp;
• Starboard Side: 2 Lpp.

2.4.4. Rigid Body Motion

Two degrees of freedom are allowed for the motion for the vessel, namely pitch and heave, allowing for dynamic trim and sinkage of the vessel.

2.4.5. Mesh

The triangle surface meshing method has been used to build the mesh. Prism layer mesher generates prismatic cell layers next to the boundary layer. The number of these layers is calculated as per:

$$m=\frac{\ln \left(1-\left(1-sf\right)\frac{\delta }{y_1}\right)}{\ln \left(sf\right)},$$

(12)

where sf=1.3 is stretch factor, δ is boundary layer thickness and y₁ is the first layer thickness. These cells help to captur0e the viscous effects in boundary layer correctly in the region which thickness is calculated as per [49]:

$$\delta =\frac{0.16L_{wl}}{\sqrt[7]{R{e}_{L_{wl}}}},$$

(13)

where L_wl is length at waterline of the ship and Re is Reynolds number (L_wl in subscription means that Reynolds number is calculated with waterline length as a reference value, otherwise L_pp is used). Trimmed cell mesher creates a volume mesh by cutting a template mesh with the surface geometry. First layer thickness is calculated for y⁺ = 150 by following guidelines [44] where it is emphasized that y⁺ should be in a range from 30 to 300. This value (150) is used as target factor to be ensured that the actual y+ values will fall within the required range.

$$y_1=\frac{150\vartheta }{u_{\tau }},$$

(14)

where υ is kinematic viscosity and u_τ is shear velocity. Shear velocity or frictional velocity is fictious quantity and it characterizes the shear at the boundary layer and it is given by:

$$u_{\tau }=\sqrt{\frac{{\tau }_w}{{\rho }_w}},$$

(15)

where τ_w is wall shear stress. This presents a contradiction: determining the thickness of the initial layer requires knowing the stress on the vessel’s hull, which in turn needs to be assessed using CFD. Hence, the recommendation from ITTC [43,50] is utilized in this scenario to compute the wall shear stress for the hull.

$${\tau }_w=0.5\left(\left(1+k\right)C_{FS}+C_A\right)V_s^2{\rho }_w,$$

(16)

where k is a form factor, C_FS is frictional resistance coefficient and C_A is the correlation allowance.

$$C_{FS}=\frac{0.075}{{\left({\log }_{10}\mathrm{Re}-2\right)}^2},$$

(17)

$$C_A=\left(5.68-0.6{\log }_{10}\mathrm{Re}\right){10}^{-3}$$

(18)

Form factor is calculated as mean value of the equations obtained by Grigson, Wright, and Conn and Ferguson given in [51]:

$$k=\frac{0.028+3.3\left(\frac{WS}{L_{pp}^2}\sqrt{C_B\frac{B}{L_{pp}}}\right)+18.7{\left(C_B\frac{B}{L_{pp}}\right)}^2+2.48C_B^{0.1526}{\left(\frac{B}{T_s}\right)}^{0.0533}{\left(\frac{B}{L_{pp}}\right)}^{0.3856}-1}{3}.$$

(19)

A grid convergence study is performed to assess the convergence of the results by systematically refining a specific input parameter in StarCCM+. In this study, three simulations are conducted, where the chosen input parameter is varied while keeping all other parameters constant. The refinement value used in this study aligns with the guidelines specified in [52] and is set to be uniform ($\sqrt{2})$. By systematically refining the input parameter and comparing the results from the different simulations, the grid convergence study aims to determine the level of convergence and establish the appropriate grid resolution for accurate and reliable results.

Convergence ratio is defined as:

$$R=\frac{{\epsilon }_{21}}{{\epsilon }_{32}},$$

(20)

where ${\epsilon }_{21}$is the difference between solution obtained using medium and fine mesh and ${\epsilon }_{32}$ is the difference between solution obtained using coarse and medium mesh. R is used for estimation of convergence conditions: monotonic convergence is achieved when 0<R<1, oscillatory convergence is achieved when -1<R<0 and divergence is achieved when |R|>1. Numerical uncertainty and error for monotonic convergence condition is estimated using generalized Richardson extrapolation (RE). The order of accuracy is calculated as per:

$$p=\frac{\ln \left({\epsilon }_{32}-{\epsilon }_{21}\right)}{\ln \sqrt{2}}$$

(21)

$${\delta }_{RE}=\frac{{\epsilon }_{21}}{{\sqrt{2}}^p-1}$$

(22)

The generalized RE solution (${\widehat{S}}_{RE}$) can be derived by:

$${\widehat{S}}_{RE}={\widehat{S}}_1-\left|{\delta }_{RE}\right|.$$

(23)

To estimate the grid uncertainty U_G, a safety factor approach is employed, defining U_G as:

$$U_G=FS\cdot \left|{\delta }_{RE}\right|,$$

(24)

where FS=1.25 is a safety factor. The normalized uncertainties are calculated as follows:

$$\overline{U}=\frac{U_G}{{\widehat{S}}_{RE}}\cdot 100\%.$$

(25)

2.4.6. Post Processing

In the numerical simulations conducted for this study, the direct consideration of roughness effects and air resistance resulting from the presence of a superstructure is not included. However, these effects are accounted for in the post-processing stage by following the recommended procedures and guidelines outlined in [53]. The ITTC guidelines provide specific methodologies for incorporating the effects of roughness and air resistance into the analysis, allowing for an assessment of the overall performance and characteristics of the ship. Roughness allowance (∆C_F) is calculated as per [54]:

$$\Delta C_F=0.044\left[{\left(\frac{k_s}{L_{WL}}\right)}^{\frac{1}{3}}-10{\mathrm{Re}}^{-\frac{1}{3}}\right]+0.000125$$

(26)

and air resistance coefficient ($C_{AAS}$) [50] is:

$$C_{AAS}=0.8\frac{{\rho }_A\cdot A_{VS}}{{\rho }_w\cdot WS},$$

(27)

where k_s indicates the roughness of hull surface equals to 150 ∙10^-6 m [54], and A_VS is the projected area of the ship superstructure to the transverse plane. New total resistance is calculated as:

$$F_T^{\prime }=\frac{1}{2}{V}_s^2\cdot WS\cdot {\rho }_w\cdot \left(C_T+\Delta C_F+C_{AAS}\right),$$

(28)

where C_T refers to total resistance coefficient obtained as a result from CFD simulation.

New thrust is calculated as per:

$$T^{\prime }=\frac{R_T^{\prime }}{1-t^{\prime }}$$

(29)

Together with propeller performance data (K_T, K_Q, η₀(J)), new operating point is obtained by numerical solving of equation:

$$\frac{T^{\prime }}{{\rho }_w\cdot V_A^2\cdot D^2}{J}^2-K_T\left(J\right)=0$$

(30)

Once the advance ratio coefficient is obtained, the propeller speed (n) and torque (Q₀) can be evaluated:

$$n=\frac{V_A}{JD}=\frac{\left(1-w^{\prime }\right)V}{JD},$$

(31)

$$Q_0=K_Q{\rho }_w{n}^2{D}^5,$$

(32)

where K_Q is evaluated from the propeller performance curve and Q₀ is torque in open water condition. Hence actuator disk approach has been used in this case, relative rotative efficiency is adopted to be 1, that means that torque is calculated as:

$$Q=\frac{Q_0}{{\eta }_R}\approx Q_0$$

(33)

Finally, brake power can be evaluated as:

$$P_b=\frac{2n\pi Q}{{\eta }_S},$$

(34)

where η_S is shaft efficiency, adopted to be 0.99. Hull efficiency is considered to be the same as in steady resistance simulations i.e., w’=w and t’=t, where w’ and t’ are wake fraction coefficient and thrust deduction factor in post processing.

Calibration factors C_n and C_Pb are calculated as the average ration between model test (MT) results and CFD results for propeller speed (n) and brake power (P_b).

$$C_n=\frac{n_{MT}}{n_{CFD}}$$

(35)

$$C_{Pb}=\frac{P_{b_{MT}}}{P_{b_{CFD}}}$$

(36)

2.5. Trim Optimization

Trim optimization was conducted for drafts of 7.5m, 8m, and 8.7m at speeds of 12.5kn, 15kn, and 18kn at 7 different conditions (trims). The range of trims goes from -1.5m to 1.5m with a 0.5m step. Negative sign means trim by bow and positive sign, trim by stern. The outcome of trim optimization in CFD analysis typically involves determining the optimal trim settings for a vessel at considered drafts and speeds. This optimization aims to minimize resistance and ensure optimal performance under different operating conditions. By adjusting the vessel’s trim, operators can achieve better hydrodynamic characteristics, leading to improved overall performance and operational efficiency. With available specific fuel oil consumption (SFOC) curve as a function of P_b, a possible reduction of fuel consumption was also calculated.

2.6. Mathematical Model for Assessing the Outcomes from Trim Optimization

First outcomes from CFD simulations are brake power and propeller speed for one ship speed and draft i.e., trim i.e., displacement. The results are presented in the form of relative differences between the brake power of the trimmed ship and the ship on an even keel. A negative sign indicates that there is a saving, i.e., a reduction in the required brake power for that case. A graphical example of the results is shown in Figure 3a, where the relative brake power savings are depicted as an area with ship speed on the x-axis, trim on the y-axis, and relative power savings on the z-axis. The optimal (favorable) sailing zone is depicted in green on the same diagram, meaning that at current trim and speed, less power is required compared to when the ship is sailing on an even keel. The yellow zone represents the transitional zone where the change in required brake power is negligible, while the red zone is the zone to be avoided because more power is needed for the ship to sail at the same speed compared to sailing on an even keel. Each point on the surface in Figure 3 a corresponds to a certain displacement, where it depends solely on the trim. For easier calculation, a surface Δ(trim, V_s) is formed from each line Δ(trim), and surfaces for all three considered drafts are shown in Figure 3b. In order to have smooth surfaces, all values between estimated from CFD analysis are linearly interpolated. Therefore, initial matrix 3x7 (3 speeds, 7 trims) is reshaped to 151x111 with a speed step 0.05kn and trim step 0.02m.

Additionally, values between drafts of 7.5m to 8m and from 8m to 8.7m are linearly interpolated with a step of 0.01m. Thus, from the initial matrix of 3x7x3 (3 speeds, 7 trims, 3 drafts), matrices (Σ_m) of size 151x111x120 are obtained where m is Δp_b or Δ.

The daily fuel oil consumption (DFOC) is estimated based on the determined brake power from CFD simulation and the available specific fuel oil consumption (SFOC) curve for the engine installed on the ship. The SFOC curve is approximated by a sixth-degree polynomial:

$$SFOC=-5.54\cdot {10}^{-23}{P_b}^6+3.75\cdot {10}^{-18}{P_b}^5-8.28\cdot {10}^{-14}{P_b}^4+8.00\cdot {10}^{-10}{P_b}^3-3.11\cdot {10}^{-6}{P_b}^2-9.63\cdot {10}^{-4}{P}_b+209$$

(37)

The mentioned polynomial provides SFOC values with an average error of 0.02% compared to the exact available values from reference document.

With evaluated DFOC for each brake power, the same diagrams as they are shown in Figure 4 can be created (Figure 5).

In operations, after loading cargo onto the ship, only the displacement is known. The goal of the conducted trim optimization is to find the optimal position, i.e., trim, and recommend the sailing speed based on that single parameter. The optimal trim and recommended speed are those at which the greatest savings in brake power, and consequently fuel consumption, are achieved. The mathematical expression of the above statement is as follows:

$$\delta P_{b_{min}}=\underset{ijk}{\min }\left(\delta P_b^{\left(ijk\right)}\left|\left|{\Delta }^{\left(ijk\right)}-{\Delta }_{LC}\right|\le \epsilon \cdot {\Delta }_{LC}\right.\right),$$

(38)

where ${\delta P}_{b_{min}}$is the biggest brake power reduction where the corresponding displacement ${∆}^{\left(ijk\right)}$ is within a tolerance ($ϵ$=0.0001) of the target displacement (${∆}_{LC}$). Index i represents speed, j represents trim and k represents brake power reduction i.e., displacement in matrix Σ_m. When ${\delta P}_{b_{min}}$is found, it is easy to extract indices i, j and k and therefore to derive optimum trim and speed for corresponded displacement i.e., mean draft. In practice, there may be a need for speed as an input parameter, which depends on operational requirements, and this case will be addressed in subsequent work.

2.6.1. Mathematical Model for Assessing Brake Power, DFOC and Propeller Speed

A sophisticated mathematical model utilizing Artificial Neural Networks (ANNs) has been developed based on CFD results. By specifically addressing the challenges posed by non-simulated loading conditions, this model offers a solution for estimating parameters such as brake power, propeller speed, and DFOC.

At its core, the model employs a feedforward artificial neural network (FFANN), a fundamental architecture in machine learning and deep learning. The FFANN’s structure comprises interconnected layers of neurons, including an input layer, one or more hidden layers, and an output layer (Figure 6). Each neuron is equipped with adjustable parameters called biases, which introduce flexibility by allowing for offsets or shifts in output. During training, these biases are fine-tuned to better fit the data and ensure accurate estimation of the desired parameters.

One of the key features of the FFANN is its one-directional flow of information, from the input layer through the hidden layers to the output layer. This design simplifies the learning process and enables efficient modeling of complex data relationships. Additionally, activation functions embedded within neurons introduce non-linearity, enabling the network to capture patterns present in real-world data.

Overall, this innovative approach represents a significant step forward in the advancement of engineering solutions. It not only addresses the inherent challenges of non-simulated loading conditions but also sets a precedent for the application of artificial intelligence in optimizing the performance of complex systems. As further research and development are conducted in this field, the potential for ANNs to revolutionize engineering practices and transportation systems continues to expand.

In the context of a feedforward artificial neural network (FFANN), the connections between neurons are governed by weights, which dictate the strength of the connections. These weights are adjusted during the training process using an algorithm called backpropagation. The name “backpropagation” stems from its mechanism of iteratively propagating errors backward through the network to adjust the weights and biases of neurons.

The training process ©nvolves two phases: the forward pass and the backward pass. During the forward pass, input data is fed into the network, and predictions are made based on the current weights and biases. These predictions are compared to the actual outputs, and the resulting error is calculated.

In the backward pass, this error is propagated backward through the network, layer by layer. The gradient of the error with respect to each weight and bias is computed using techniques like the chain rule of calculus. These gradients indicate the direction and magnitude of adjustment required to minimize the error. The weights and biases are then updated accordingly.

At each neuron, the weighted sum of inputs is computed by multiplying each input value by its corresponding weight and summing the results. This sum is then combined with the neuron’s bias:

$$S_{x,i}=f\left({\displaystyle \sum _i{S}_{x-1,i}\cdot w_{x,i}+b_{x,i}}\right),$$

(39)

where $f$—activation function, $S_{x,i}$—output signal, of $i$^th neuron in $x$^th layer, $S_{x-1,i}—$output of $i$^th neuron in $x-1$th layer, $w_{x,i}$—weights connected to neuron and $b_{x,i}—$corresponding bias.

The resulting value is passed through an activation function, which introduces non-linearity into the network. Common activation functions include sigmoid, tanh and others, each serving different purposes in facilitating learning and modeling complex data relationships. In this case, sigmoid function is used:

$$f=\frac{1}{1+e^x}.$$

(40)

The general form of the function is:

$$Y_u=\frac{f_p\left(p+{\displaystyle {\sum }_{w=1}^r\left(P_w\cdot f_{p-1}\left(\dots \cdot f_2\left(b_i+{\displaystyle {\sum }_{j=1}^m\left(B_{ij}\cdot f_1\left(a_j+{\displaystyle {\sum }_{k=1}^n\left(A_{jk}\cdot \left(P_k\cdot X_k+R_k\right)\right)}\right)\right)}\right)\dots \right)\right)}\right)-G_u}{L_u},$$

(41)

where:

• X_k—Input parameters (number of neurons in input layer);
• Y_u—Output parameters (number of neurons in output layer);
• f₁,f₂,…,f_(p-1),f_p—activation function;
• k,j,©,…,v,w—number of neurons in each layer;
• A,B,…,P—weights between layers;
• a,b,…,p—weights between neurons in each layer and bias.

The software that has been used for training the ANN in this study is called aNETka [55] developed in LabVIEW. The output file from this software is not direct mathematical model, just weights coefficients and four additional parameters called Down offset Input (DI_k), Down offset Target (DT_u), Up offset Input (UI_k) and Up offset Target (UT_u).

$$D{I}_k=\min \left(X_k\right)$$

(42)

$$D{T}_u=\min \left(Y_u\right)$$

(43)

$$U{I}_k=\max \left(X_k-\min \left(X_k\right)\right)$$

(44)

$$U{T}_u=\max \left(Y_u-\min \left(Y_u\right)\right)$$

(45)

The rest of parameters are calculated as follows:

$$P_k=\frac{0.9}{U{I}_k}$$

(46)

$$R_k=-\frac{0.9\cdot D{I}_k}{U{I}_k}+0.05$$

(47)

$$L_u=\frac{0.9}{U{T}_u}$$

(48)

$$G_u=-\frac{0.9\cdot D{T}_u}{U{T}_u}+0.05$$

(49)

In the trim optimization study, data obtained from CFD analysis, including brake power, DFOC, and propeller speed, were utilized as inputs for training an Artificial Neural Network (ANN). This ANN establishes correlations between draft, speed, trim, displacement, and the aforementioned performance parameters.

The applicability range of the model is defined within specific boundaries for each parameter: draft (T_s) ranges from 7.5 to 8.7 m, speed (V_s) ranges from 12.5 to 18 kn, trim ranges from -1.5 to +1.5 m and displacement (Δ) ranges from 18079 to 22884m³. The input dataset for training consists of a 63x7 matrix, where number 63 represents the total number of CFD simulations conducted and number 7 denotes to each parameter. Four parameters are designated as the input layer (draft, displacement, trim, speed), while the remaining three represent the output layer (propeller speed, brake power, DFOC).

Prior to training, all data were normalized within the range of 0.05 to 0.95 to ensure the stability of activation functions, which can be sensitive to extremely small or large input values. This normalization step enhances the convergence and effectiveness of the training process.

For testing purposes, 8% of the total input parameters (5 out of 63) were reserved, while the remaining data was allocated for training. The training process was programmed to iteratively adjust the weights and biases of the network, stopping when the Root Mean Square (RMS) percentage error reached the target threshold of 1%:

$$RMS=\sqrt{\left(\frac{{\displaystyle \sum \left[{\left(\frac{TV-OV}{TV}\right)}^2\right]}}{N}\right)},$$

(50)

where TV—target values, OV—ANN output values, $N$—number of values

When creating a neural network model, it’s vital to carefully consider the number of neurons in each layer and the presence of hidden layers. These decisions must be well-informed and tailored to match the complexity of the problem at hand and the specific attributes of the dataset.

2.7. Trim optimization software application

Formulas derived for a particular purpose may not always be user-friendly, leading to the creation of an application to facilitate their use. In this case, an application developed using MATLAB’s App Designer package offers a user-friendly interface for utilizing these formulas. The application allows users to input only one parameter i.e., displacement and obtain corresponding output values efficiently: optimum trim, speed, propeller speed, brake power and DFOC. This enables users to quickly obtain results without the need for manual calculations.

3. Results

The results in this section will be presented by subsections in the same order, i.e., following the phases outlined in the section 2 (Methods).

3.1. 3D Modeling

In Figure 7, a 3D model and characteristics sections of waterlines (a), buttocks (b) and frames (c) are depicted in various projections obtained based on construction drawings.

The propeller underwent a modeling process, aligning with the specifications outlined in the reference drawing. A comparison was conducted between the 3D model of the propeller and the reference drawing to ensure precise adherence to the intended design. The findings of this comparison are visually presented in Figure 8, offering a clear representation of the propeller’s details and shape. The 3D model of the propeller (Figure 9a) enables engineers and designers to perform CFD simulations and evaluate performance characteristics. The rudder was modeled to faithfully represent its original design but without any gaps between rudder shaft and rudder shell (Figure 9b).

3.2. First Level Verification

The reliability of the developed 3D model was assessed by comparing its calculated hydrostatic parameters against the corresponding ones given in the Trim & Stability booklet. The corresponding percentage deviations between C_B, V, LCB and WS are listed in Table 3.

It can be easily observed that deviations between calculated and reference values of hydrostatic data are less than 1% especially for the range of interest (drafts 7.5m, 8.0m and 8.7m).

3.3. Open Water Test

Graphical representation of CFD results (K_T, K_Q, η₀) from conducted OWT with B-series propeller with data obtained with mathematical model is presented in Figure 10a, while relative differences between each evaluated parameter is presented in Table 4.

As the results for each of the advance coefficients are within the prescribed 3%, the CFD calculation methodology is considered valid, therefore, with the same setup, original propeller OWT was conducted and results are presented in Figure 10b.

3.4. Second Level Verification

Obtained results from CFD OWT simulation with original propeller are used as an input propeller characteristics data in self-propulsion simulations. First group of simulations were conducted at design draft (7.8m) and speeds of 17kn, 19kn and 21kn because for these draft and speeds model test data are available.

In the Figure 11 meshed domain is shown. For better evaluation of free surface, additional refinement zones are defined in wave field zones and near the hull (Figure 12a,b).

Results from completed self-propulsion simulations and steady resistance simulations by following procedure described in section 2.4 are presented in Table 5 where R_p is pressure resistance and R_v is viscous resistance.

Results from CFD analysis and model test data are graphically presented in Figure 13. In Figure 13a brake power as a function of ship speed curve is depicted while in Figure 13b brake power as a function of propeller speed is presented.

Calibration factors are C_n=1.05 and C_Pb=1.04 which means that 5% deviation requirement in results compared to model test is met.

In the Figure 14 are presented y+ values for the underwater part of the hull and average y+ value is pointed out (average y+ = 70.61). Particular case is extracted from steady resistance simulation at speed of 18kn at 7.5m draft. For all other speeds y+ values do not deviate from aforementioned.

Three different meshes were considered for the grid uncertainty study with mesh sizes of 1.852 million cells (coarse), 3.628 million cells (medium), and 7.454 million cells (fine). The study was conducted at a speed of 18 knots at 7.5m draft (even keel). The results of the study, including the grid uncertainties for thrust are presented in Table 6.

3.5. Trim Optimization

Trim optimization CFD study was performed at 7.5m, 8m and 8.7m draft for a range of speeds 12.5kn, 15kn and 18kn at 7 different loading (trim) conditions and coarse mesh. Therefore, results for obtained propeller speed, brake power and estimated DFOC are presented in Table 7, Table 8 and Table 9, respectively.

3.6. Mathematical Model for Assessing the Outcomes from Trim Optimization

Brake power reduction and DFOC reduction are linearly interpolated for additional trims and speeds (values between considered in CFD analysis) and results are graphically presented for three considered drafts in Figure 15 and Figure 16. Optimum trim and speed are also depicted with red points in Figure 16. Black line in Figure 15 and Figure 16 represents neutral line, i.e., the state of even keel against which brake power savings and DFOC savings are calculated.

The biggest reduction in DFOC can be achieved at 1.5m trim by bow, 7.5m draft and speed of 15.17kn (10.5%), while for the 8m and 8.7m draft, the biggest reduction in brake power is at 12.5kn (8%, 7%, respectively). Sailing at 12.5kn instead of 15kn or 18kn, and maintaining an even keel, can result in greater fuel savings for higher drafts. Although greater fuel savings might be achievable with a greater trim by bow, but any condition more than -1.5m trim by bow cannot be attained on the current ship.

In Figure 17, wave fields are depicted for the case of V_s=15kn, T_s=7.5m at trim=0m (Figure 17a) and at trim=-1.5m (Figure 17b). It is evident that the total wave amplitude when trimmed towards the bow by 1.5m is reduced by 0.25m, indicating a decrease in pressure resistance—displacement ratio for 15% (see Table 10).

In the Figure 18 are presented free surfaces in center line and along the hull at two different conditions: trim = 0m (white hull) and trim = -1.5m (yellow hull). Black line corresponds to free surface level at zero trim condition while red line corresponds to free surface level at trim by bow for 1.5m. In the Figure 18a are presented bow waves, while in the Figure 18b are presented stern waves. Peak-to-peak amplitude of the bow wave is reduced by 0.64m while the stern wave is reduced for nearly 1m when the ship is trimmed towards the bow by 1.5m.

For others drafts and trims, the impact of bow trim remains similar, nevertheless a less noticeable reduction in pressure resistance. It seems that this kind of the bulb has better effect when it is more submerged into the water at lower drafts. This contradicts what is highlighted in numerous literature sources, including [56] and [57] but it is not the first time that the same conclusion with fully submerged bulb has been set [58,59].

3.6.1. Mathematical Model for Assessing Brake Power, DFOC and Propeller Speed

Unfortunately, a single mathematical model did not attain the desired RMS within a reasonable training timeframe and iteration count. As a result, two separate mathematical models were devised: one focusing on brake power and DFOC and the other dedicated to propeller speed prediction. The formula for brake power follows a general form:

$$P_b=\frac{f_3\left(c_{Pb}+{\displaystyle {\sum }_{i=1}^7\left(C_{iPb}\cdot f_2\left(b_i+{\displaystyle {\sum }_{j=1}^{10}\left(B_{ij}\cdot f_1\left(a_j+{\displaystyle {\sum }_{k=1}^4\left(A_{jk}\cdot \left(P_k\cdot X_k+R_k\right)\right)}\right)\right)}\right)\right)}\right)-G_{Pb}}{L_{Pb}}.$$

(51)

The general form of the function for DFOC is:

$$DFOC=\frac{f_3\left(c_{DFOC}+{\displaystyle {\sum }_{i=1}^7\left(C_{iDFOC}\cdot f_2\left(b_i+{\displaystyle {\sum }_{j=1}^{10}\left(B_{ij}\cdot f_1\left(a_j+{\displaystyle {\sum }_{k=1}^4\left(A_{jk}\cdot \left(P_k\cdot X_k+R_k\right)\right)}\right)\right)}\right)\right)}\right)-G_{DFOC}}{L_{DFOC}}.$$

(52)

The general form of the function for propeller speed is:

$$n=\frac{f_3\left(c_n+{\displaystyle {\sum }_{l=1}^5\left(C_{ln}\cdot f_2\left(b_l+{\displaystyle {\sum }_{m=1}^8\left(B_{lm}\cdot f_1\left(a_m+{\displaystyle {\sum }_{s=1}^4\left(A_{ms}\cdot \left(P_s\cdot X_s+R_s\right)\right)}\right)\right)}\right)\right)}\right)-G_n}{L_n}.$$

(53)

Configurations of neurons for both mathematical models are presented in the Figure 19a,b.

Obtained coefficients for Equations (51) and (52) are given in Table 11, Table 12, Table 13, Table 14 and Table 15, while coefficients for Equation (53) are given in Table 16, Table 17, Table 18 and Table 19.

The standard deviations of the relative differences in brake power, DFOC, and propeller speed calculated values using Equations (51)–(53), and the results obtained from the CFD analysis for the same three parameters are 0.6%, 0.6%, and 0.2%, respectively.

3.7. Trim Optimization Software Application

Development of software focusing on trim optimization is another step forward that demonstrates how engineering and technical practice can be more effectively connected. The equations obtained (51), (52), and (53) are very complex and not suitable for manual calculation. Additionally, the matrix obtained through interpolation of the CFD analysis results contains over 2 million elements (a 151x111x120 matrix), making tabular or graphical representation of the results will be impractical for the average user. Therefore, the most elegant solution is to develop an application with a graphical interface. Today, there are many software tools that operate based on writing in various programming languages, but this task is addressed to other engineering disciplines. This paper demonstrates the development of a software tool within MATLAB App Designer. MATLAB App Designer is highly suitable for rapid development of simple applications as it has certain functions predefined within the Component Library. In this specific case, only three predefined functions were used: “Button,” “Edit Field,” and “Label.” The interface design is depicted in Figure 20:

In the “Input displacement” label user should insert current displacement in cubic meters in the range from 18078 to 22884m³. The estimation of the optimal trim, corresponding mean draft, estimated speed, brake power, DFOC, and propeller speed is performed by simply clicking the “Calculate” button. In case the user inputs a displacement below or above the specified values despite the stated applicability limits of the mathematical model and initiates the calculation, a message will appear on the screen: “Displacement is out of range” (see Figure 21).

If a displacement within the specified limits is entered, and the “Calculate” button is clicked, all the data is instantly displayed (see Figure 22).

The application provides results based on previously obtained mathematical models for estimating the optimal trim, speed, corresponding mean draft, brake power, DFOC, and propeller speed. There is also a note that sign “-“ in Optimum Trim label corresponds to trim by bow. Due to the specified project task, which required trim optimization at certain speeds, trims, and drafts rather than displacements, it is clear that for the same draft and different trims, the vessel has different displacements. Considering the upper limit of the performed calculations with a draft of 8.7m and taking into account various trims, accordingly, the upper limit of displacements for which the calculations were performed at this draft (8.7m) is variable (refer to Figure 4b). Considering the obtained results, there is a displacement range where the optimal trim is 0m, meaning that sailing on an even keel is recommended (see Figure 23). This zone starts from the displacement corresponding to a draft of 8.7m and even keel position.

The mathematical model does not take draft as an input parameter but determines the mean draft by seeking the optimal trim for a given displacement. It’s logical to apply the recommendation of optimal trim obtained from the mathematical model only through additional ballasting of the vessel if achieving the optimal trim cannot be attained with the existing ballast water in the tanks (if any). This is manifested by further increasing the displacement. Therefore, it’s necessary to rerun the application with the new displacement. Through a few iterations, the optimal trim can be achieved, thus obtaining the estimated ship speed, propeller speed, brake power and DFOC in a reliable manner. It should be noted that the application is not part of the Loading Computer nor does it track stability data for a specific loading condition for recommended optimal trim. However, it can be integrated into the Loading Computer system to enhance its functionality.

4. Discussion and Conclusions

In this work, not only trim optimization aiming to improve the energy efficiency of the ship, specifically the CII parameter by reducing fuel consumption, was presented. The methodology for approaching this issue using CFD software, the application of results obtained from CFD analysis, i.e., using ANN to develop a mathematical model that provides estimates of parameters (brake power, DFOC, and propeller speed) for all conditions for which CFD calculations were performed, as well as for all intermediate conditions, was also introduced. Since the use of ANN resulted in a complex mathematical model that is not straightforward to use in engineering practice, it went a step further, incorporating it into a simple software tool (application).

The entire trim optimization process is described through 7 phases that should be applied to other types of ships in future research to develop a universal methodology, not only for trim optimization but also for CFD calculations based on which various engineering analyses such as optimization of bow or stern, determination of Vref for EEXI, efficiency of energy-saving devices, etc., can be conducted. For now, the applied methodology, which includes 3D modeling of a digital twin whose hydrostatic parameters differ by up to 1% when comparing the same parameters with data from the Trim & Stability booklet, has proven to be sufficiently good. With additional assumptions about determining the boundary layer thickness and shear wall stress, i.e., shear velocity, satisfactory accuracy of the y+ parameters, which are very significant for obtaining reliable results, was achieved. The obtained values of y+ (approximately 70-72) differ from the initial target y+ of 150, which is to be expected because the equations used to estimate both the boundary layer thickness and wall shear stress are not the only ones that can be applied and are very general. For example, the formula for boundary layer thickness applies to a flat plate, but in calculating shear stress, a form factor is included, which is itself obtained based on empirical expressions. Therefore, these formulas cannot provide exact parameter values because these parameters are actually calculated by CFD analysis. However, they can be used to determine enough fine mesh (parameters) so that the results obtained by CFD analysis fall within prescribed limits, according to currently published requirements (up to 5% difference in power brake and shaft speed) for a specific ship case.

The application of ANN to develop a mathematical model has proven to be very reliable. Generally speaking, artificial intelligence (AI) is increasingly finding applications in various engineering industries, where new design approaches can be developed and existing ones significantly accelerated based on machine learning. It should be noted that the obtained mathematical model applies to navigation in calm waters, which is practically an almost impossible scenario. This approach, which encompasses multidisciplinarity (CFD analysis, ANN application, software tool programming), can lay the foundation for further development of the shipbuilding industry from an engineering perspective.

The further idea is to collect real (measured) data on fuel consumption, engine power, ship speed, and navigation conditions during voyages. Monitoring and recording of some of these parameters have already become mandatory through [10]. Based on the stored parameters and the formed database, and with the help of machine learning and AI, it will be possible in the future to predict in advance to which energy class the ship belongs and what can be changed during real-time operation to make the ship “green.” This can lead to a global reduction in exhaust gas emissions instantly, without ship owners being exposed to penalties imposed by the system because the CII is determined retroactively.