Design of Energy-Saving Duct for JBC to Reduce Ship Resistance by CFD Method

Abstract

The present work is to modify the stern duct of Japan bulk carrier (JBC) to reduce the ship total resistance for Froude number 0.142. The research method is to use the CFD (computational fluid dynamics) tool OpenFOAM to predict the total resistance and viscous flow field around the ship hull considering the free surface. First, the V&V (verification and validation) analysis was conducted to ensure our CFD method is reliable, and to estimate its error and uncertainty range for the original JBC with and without the original stern duct. Second, the angle of attack (AoA) of the foil section of the stern duct was adjusted, thus forming a series of different ducts axisymmetric along the propeller shaft axis. The CFD simulations were performed for each duct appended behind the ship hull. By comparing the total resistance with and without the duct, the ship hull resistance reduction effect caused by the duct would be obtained. Finally, the duct with 7$°$ AoA achieved the best resistance reduction of 2.49%. The original duct (with 20$°$ AoA) could only provide a 0.6~0.7% reduction. The nominal wake was also analyzed with the ducts, in order to understand the detailed flow phenomena which had an essential influence on propeller inflow.

1. Introduction

Based on the awareness and fact of worsening global warming and limited natural resources, stricter environmental regulations have been addressed to shipbuilding industries and maritime transportations. If ships can advance in sea more efficiently, this is a great concern for ship owners and designers; reduction of ship resistance or improvement of propulsion efficiency would be the solution. To achieve the above-mentioned goal, many energy-saving devices (ESDs) have been designed and installed ahead, after or even on the propeller for commercial ships. It can also be a combination of several ESDs [1,2]. One type of ESD is a stern duct which is also called wake equalization duct (WED) or pre-duct, i.e., the duct in front of the propeller. Normally, this is a nozzle-shaped, accelerating duct to increase the inflow velocity and make the inflow more uniform for the propeller behind, which would improve propulsion efficiency. According to the experimental data released in the workshop on CFD (computational fluid dynamics) in ship hydrodynamics in Tokyo 2015 (T2015) [3] for Japan bulk carrier (JBC) with a stern duct, not only was the propulsion efficiency improved, but the bare hull resistance could also be reduced. Without the propeller, the difference between total resistances with and without the stern duct was 0.6%, i.e., 0.6% resistance reduction.

JBC is a relatively new design of ship hull form compared with many other geometries such as Wigley hull, S60 container ship, KCS (KRISO container ship), KVLCC2 (KRISO very large crude-oil carrier 2) and DTC (Duisburg test case). With the higher block coefficient (CB = 0.858; 0.8098 for KVLCC2, 0.6505 for KCS), studying JBC is more challenging and necessary, since its resistance would be larger and its ship wake would be more non-uniform. Furthermore, it is worth understanding if an ESD works effectively on such full hull form. The ESD on JBC is a stern duct, which is axisymmetric along the propeller axis with NACA 4420 foil section.

Maasch et al. [4] developed a general-purpose wake analysis tool (WAT) to analyze ship’s nominal wake. The CFD velocity field for JBC with and without duct in calm water was an input for WAT. Their CFD tool was Star-CCM+. The V&V (verification and validation) analysis focused on the total resistance coefficient Ct and wake factor (1-w) of 0.7–1.0 propeller radius. Through the V&V analysis, the error E%D against experimental data for the total Ct was 2–3%, which was slightly higher than the experimental uncertainty, which was 1%. Furthermore, E%D = 2–4% for the ship sinkage and 4% for the ship trim. Five B.S.R.A. (British Ship Research Association) criteria were used to evaluate the quality of nominal wake. The wake non-uniformity is evaluated by Δw, which is the average difference of wake fraction w in different propeller radius. As many studies [3,5] for JBC had pointed out, the occurrence of a large flow separation area behind and below the duct was also confirmed. They further concluded that with the duct, the boundary layer thickness of the ship bottom around the stern region was smaller. The shear force coefficient was also relatively small. Within the 0.6% total Ct reduction, the hull part of Ct was reduced by 1% and the pressure component of the total Ct decreased by 3.4%. Because of the flow separation, only the flow velocity around the outer part of the duct was slightly uniform, but the nominal wake quality became worse. The B.S.R.A. criteria could not all be satisfied, implying a certain level of cavitation and vibration that the propeller would suffer. Lowering ship speed and draft from the design condition, i.e., 14.5 knots and 16.5 m draft for full scale, to 4.5 knots and 2 m, for the average Δw the 10% change was found without the duct, but the change was only 3% with the duct. This meant that by installing the duct, the influence of different ship operating conditions on the ship wake became less. However, the B.S.R.A. analysis showed that all propeller operating conditions became worse. It was concluded that operating JBC should not be far away from the design condition.

Furcas et al. [6] proposed a SBDO (simulation-based design optimization) method using CFD to design WED, i.e., the stern duct on the JBC. The channel of the duct reduced the wake velocity loss, which would create a positive and beneficial interaction between the ship and propeller, and would eventually provide extra thrust. The design method was based on the parameterized duct geometry and the simulation of Star-CCM+ with the target of maximal thrust generated by WED. The result showed the automatic design method could be applied in practical use. They found that designing the pre-duct against the ship stern wake profile would also improve the propulsion efficiency in general.

In fact, other than JBC, the stern duct has been investigated for different types of commercial hull forms for a long time. Çelik [7] first designed a WED using lifting line method for a chemical ship with CB = 0.858, and then utilized the CFD software Fluent to model the viscous flow field. The design was considered with the duct location, inner diameter, foil section and its angle. Eventually, the duct’s axisymmetric axis was located above the propeller shaft center. Thus, the duct outlet and more uniform flow behind the duct would encounter the propeller blades, while rotating up to the upper area. The maximal 9.7% increase for propulsion efficiency was gained. Guiard & Leonard [8] patented the Becker Mewis Duct^®, which achieved 6.3% average power reduction in the ship model towing tank test. The design tool was also CFD software Star-CCM+. They pointed out that the greatest challenge was to predict the ship wake in full scale, especially for the full hull form. In their work, Bertram’s correction formula [9] was used as a benchmark to predict the full-scale ship wake from the model-scale result done by CFD. Among many simulation setups, the full-scale CFD with SST (shear stress transport) k-ω and 0.5 mm roughness height provided the closest result. The difference of nominal wake fraction was 0.5%. The new design named Becker Twisted Fins^® (BTF) was also proposed. This was an axisymmetric duct in front of the propeller, but equipped with pre-swirl fins extending from the inside of the duct to the outside. Additionally, Winglet was in the fin tips. Those fins were distributed asymmetrically between port and starboard sides. Chang et al. [10] used Star-CCM+ to study a Mewis duct on a bulk cargo ship. The angle of attack (AoA) for the duct section was 15$°$, and only covered the 214$°$ upper part of sector. There are four inner guide fins inside the duct, and the intersection angles between neighbor fins were different. The duct and fin chord length were 0.28 and 0.18 propeller radius, respectively. The fin foil section was NACA0015 and the proper AoA for each fin was investigated. For the fin with 5–12$°$ AoA, the propeller efficiency would be higher than that without the duct. Although the propeller thrust was increased by the duct, the ship total resistance was also higher. Thus, the whole system needed to be further optimized.

Therefore, the present work will use CFD to analyze the hull resistance and flow field around the ship and duct, and to enhance the resistance reduction effect by designing the JBC’s stern duct. Furthermore, the nominal wake velocity profile on the propeller plane will be analyzed as a reference for propeller inflow velocity.

The objective of present work is to modify the stern duct of JBC to reduce the ship total resistance by using CFD software OpenFOAM. The total resistance and viscous flow field considering the free surface effect around the ship hull were analyzed. Two main stages of our study were included. The first stage was based on the original JBC with and without the original stern duct to conduct V&V analysis. The simulation error and uncertainty were obtained for a reliable and convincing CFD method. In the second stage, a series of axisymmetric ducts were constructed with AoA variation of the foil section, and then appended behind the hull for simulation. The total resistance reduction was calculated by the difference between the results with and without duct. The optimal duct with 7$°$ AoA was found for 2.49% resistance reduction. In contrast, for the original duct, the 20$°$ AoA could only achieve 0.6–0.7% reduction. To reveal the physical reason of the resistance reduction, the nominal wake was discussed among those ducts. Understanding the detailed flow field phenomena will also assist us in knowing how those ducts influence the propeller inflow.

2. CFD Methods

2.1. Geometry

Following the T2015 workshop instruction [11], the ship model of JBC hull with model scale 1:40 had the ship length L = L_PP = 7 m, beam B_WL = 1.125 m, depth 0.625 m and draft t = 0.4125 m. The JBC model was tested in a towing tank under Froude number Fr = 0.142, corresponding to 14.5 knots in full scale and model speed U = 1.177 m/s. The simulation in this study was conducted in model scale with Reynolds number Re = 7.44 $\times {10}^6$.

The ESD of JBC is an axisymmetric duct mounted on the top of the stern tube through a strut in front of the propeller. The technical drawing of the original JBC stern duct and strut is shown in Figure 1. The foil section is NACA4420, and the suction side facing inward indicates an accelerating duct. In the original design, the AoA is 20$°$. In the present work, AoA is defined based on the axis of symmetry for the duct, i.e., propeller shaft axis, as shown in Figure 2a.

In foil theory, lift coefficient is almost linearly proportional to AoA. Thus, increasing foil/duct lift by adjusting its AoA is our design method. Since the duct is in the non-uniform wake velocity field behind the ship stern, it is possible that at a certain AoA, its foil section will generate lift. The forward force component of lift will provide an extra thrust to counteract the ship resistance. As illustrated in Figure 2a, the foil trailing edge was set as a pivot point to change AoA in order to maintain the outlet area facing the propeller plane. The AoA ranged from $4°$ to $24°$ with 1$°$ or 2$°$ variation, and finally 17 different geometries of JBC hull were built with different ducts. The minimal $4°$ AoA was the smallest angle to be an accelerating duct and without contact with the stern surface, i.e., the duct inlet face area was still larger than outlet area. Through 17 CFD simulations (summary of the result in Section 3.3.1), an optimal AoA was found for the 7$°$ duct. Figure 2b shows the geometry of the 7$°$ duct. The 23$°$ AoA was the angle that the duct design no longer reduced the total resistance as AoA increased. Thus, the maximal AoA was selected as 24$°$ to confirm the situation.

2.2. Grids on the Geometry

2.2.1. Computational domain and JBC hull

Figure 3 shows the computational domain built for JBC. The domain size was based on the consideration of ship length. Only half of the computational domain, i.e., the flow field in the ship starboard side, was simulated, due to the symmetric condition. The origin (0,0,0) was set at the intersection of the ship AP (aft perpendicular) and baseline. The ship FP (forward perpendicular) was located along the +x direction, i.e., x = L. L was the ship length L_PP. The height of the undisturbed free surface would be the ship draft, i.e., z = t. The domain lengths for upstream, downstream, top, bottom and side were 1.285 L, 3.714 L, 0.536 L, 2.247 L and 2.714 L, respectively. This satisfied the recommendation of ITTC 7.5-03-02-03 guideline [12].

Figure 3 also indicated the boundary grid distribution approach. The constant spacing was implemented in the x and y direction. The geometric distribution with finer grid near free surface and ship body along z direction was applied. Initially, the simple grid block was built by the grid generator blockMesh, based on the above-mentioned boundary grid distribution. Next, the Cartesian grid was generated level by level toward the ship surface from the output of blockMesh. As the result of 6 volume refinements covering the whole ship one by one from the far field to the space around the hull body shown in Figure 2, more levels of finer grid were provided, surrounding the hull surface. In each level of refinement, one grid cell was split into two in the x and y directions, respectively. A process named Castellation in the grid generator SnappyHexMesh removed the cells inside the surface. Furthermore, the cells intersecting the surface would be trimmed and corrected to form a body-fitted grid. After that, the layered grids to resolve the boundary layer were grown upon the surface to 3 layers. In the end, the unstructured grid was constructed by SnappyHexMesh, and mainly consisted of hexahedra with some prisms, wedges, tetrahedral and polyhedral. Figure 4 illustrates that through the grid generation procedure, the finer and more levels of grid, the more could be seen distributed locally around the free surface and ship hull body.

The V&V (Verification and Validation) analysis was based on ITTC 7.5-03-01-01 guideline [13]. To perform V&V analysis, three different grid sizes were constructed: coarse, medium, and fine grid.

Table 1. Grid generation result.

		Grid Number on Boundaries *			Initial Grid Number (BlockMesh Result)	Final Grid Number (Snappyhex Result)	Layers (Coverage)	y+
		Nx	Ny	Nz	${\mathit{N}}_{\mathit{x}}\times {\mathit{N}}_{\mathit{y}}\times {\mathit{N}}_{\mathit{z}}$
No duct	Coarse	30	13	122	47,580	590,884	3 (94.5%)	87.3
	Medium	42	19	168	134,064	1,571,444	3 (99.7%)	56.6
	Fine	60	27	239	387,180	4,537,454	3 (99.2%)	37.9
With duct	Coarse	30	13	122	47,580	707,575	3 (93.4%)	84.5
	Medium	42	19	168	134,064	2,118,382	3 (94.9%)	52.4
	Fine	60	27	239	387,180	6,186,049	3 (98.4%)	37.4

* N_x: Grid number in x direction, Ny: Grid number in y direction, Nz: Grid number in z direction.

provides the details. The grid number on the boundaries increases following the refinement ratio $\sqrt{2}$ along x, y, and z direction, respectively. Therefore, the initial grid number generated by blockMesh was $\sqrt{8}$ times between two successive grid sizes. The final grid number after snapplyHexMesh generated the unstructured grid was also controlled with the ratio as close to $\sqrt{8}$ as possible. The total grid numbers for JBC hull (without duct) were 0.6 M (M: million), 1.6 M and 4.5 M from coarse to fine grid.

Figure 5 shows the grid distribution around the ship bow and stern region, and on the surface of those ship parts for three different grid sizes. Basically, the coverages of the layered grid covering the hull surface were all more than 94% (see Table 1). With the finer grid, the 3 layers of the grids can be seen to be more intact, covering the stern tube and stern part which had complicated geometry changes. Thus, their layered grid coverage could reach higher than 99%. Since the wall function was used to treat near wall flow, the thickness of the first layer, i.e., y⁺, was concerned. With the grid becoming denser, the average y⁺ reduced within the range of 38–88 (also listed in Table 1).

2.2.2. JBC Hull with Duct

The duct and strut were then added into the grid system constructed in Section 2.2 for three different grid sizes. The parameter setup for the domain boundary and grid generation remained the same. The grid generation result is presented in Figure 6 for coarse, medium, and fine grid. Table 1 shows the consistency of the grid number for the initial grid number with and without the duct. Since the higher level and finer grid density were set around the duct and strut surface, the final grid number would increase after the duct and strut were appended. The grids around the ship main body with and without the duct were nearly the same, so Figure 7 only displays the grid distribution around the ship stern and duct. The sharp and thin geometries were well-captured, such as the leading and trailing edge of the duct and strut, and the junction between them and with the stern tube (see Figure 6 for the medium grid as an example). The grid number increase due to the duct was around 0.12 M (=116,691), 0.55 M (=546,938) and 1.65 M (=1,648,595) cells for coarse, medium and fine grid, respectively. They were around 16%, 26% and 27% of the total grid number after adding the duct, which raised to 0.7 M, 2.1 M and 6.2 M. The grid details are also in Table 1. Finally, for the study of the hull with and without the original stern duct for three grid sizes, a total of 6 grid sets were built.

As explained in Section 2.1, with the AoA range from $4°$ to $24°$, a total of 16 + 1 grid systems with the JBC hull were built for duct design. The “+1” was the $20°$ duct which was the original design and its grid had existed in Section 2.2.1. Through 17 CFD simulations, an optimal AoA was found for the 7$°$ duct (details in Section 3.3.1). The grid distribution around the 7$°$ duct as it was being appended on the ship was illustrated in Figure 8.

2.3. Boundary Conditions

For boundary conditions, on the domain of the inlet, the uniform inflow velocity was specified as U. The non-slip condition was used on hull surface with the wall function for viscous turbulence ${\mu }_t$, turbulent kinematic energy k and turbulent dissipation rate ω. Wall roughness was not considered. On the inlet, the constant ${\mu }_t$, k, ω, and total pressure $p_0$ were estimated based on the usual assumption of the far field freestream for external flow. Volume fraction α = 0.5 was set at z = L (ship draft) to model the free surface. The symmetry condition was applied along the y = 0 plane (mid-plane) for flow symmetric condition. It also served as a zero-gradient condition in the surface normal direction on the domain side and bottom. The zero-gradient condition was utilized for static pressure p and α on the hull, p on the outlet, ${\mu }_t$ on the outlet and top. Some variations of the zero-gradient condition, including reverse flow treatment and water height conservation, were provided in OpenFOAM and applied on the other boundaries.

2.4. CFD Algorithms

The analysis and simulation tool in the present work is the open source CFD software OpenFOAM (open-source field and manipulation) version 6 for viscous flow simulation considering the free surface. An incompressible, two-phase flow model VOF (volume of fluid [14]) was applied to solve the flow field around the ship, including velocity, pressure field, and free surface elevation. The turbulence model is SST (shear stress transport) k–ω, i.e., SST-2003 [15], and the velocity and pressure coupling algorithm is PIMPLE [16]. It was exclusively developed by OpenFOAM to combine PISO (pressure implicit with splitting of operator) and SIMPLE (semi-implicit method for pressure linked equations).

The numerical schemes selected to discretize Navier–Stokes equations and the above-mentioned mathematical model of flow field modeling were briefed as below. The temporal discretization method was the 1st order Euler implicit with local time stepping. The gradient term was calculated by the 1st order central difference. The 2nd order upwind method was used mainly for the divergence term. The Laplacian term was computed by linear interpolation. The gradient of surface normal direction was obtained by a 2nd order explicit scheme with non-orthogonal corrections. The control point value of a volume or surface was also interpolated linearly.

3. Results

3.1. Verification and Validation

Table 2. Verification result for JBC hull without/with the original duct.

1000Ct	S1	S2	S3	RG	Verified or Not
Without duct	4.400	4.577	5.961	0.128	Yes
With duct	4.393	4.544	5.886	0.113	Yes

shows the verification result for the ship resistance of JBC hull without/with the original duct. In the Table, S₁, S₂ and S₃ present the total ship resistance coefficient 1000Ct for the fine, medium, and coarse grid, respectively. Ct is calculated as follows

$$Ct=\frac{Rt}{0.5\rho S{U}^2},$$

(1)

where Rt is the force value (N) of resistance and $\rho$ is the water density (998.2 kg/m³). The wetted surface area S = 12.2225625 m² for JBC bare hull. S = 12.2711875 m² is for JBC with duct including the original duct. It is also used later in Section 3.3.1 for different duct design, since the area difference among ducts is small compared to the hull.

To check the grid convergence, a ratio RG is referred as below:

$$RG=\frac{S_2-S_1}{S_3-S_2},$$

(2)

Once RG is less than 1, it means the resistance difference between the medium and fine grid is smaller than the difference between medium and coarse grid. The so-called monotonic convergence is achieved: as the grid number increases, the resistance difference decreases between two grid densities. Therefore, our CFD method is verified as shown in Table 2, and the grid independence is confirmed.

For validation, the grid uncertainty U_G is estimated by the suggestion of ITTC 7.5-03-01-01 guideline [13] and listed in

Table 3. Validation result for JBC hull without/with the original duct.

E%D	S1	S2	S3	D	UD	Uv%D	Validated or Not
Without duct	−2.59%	−6.72%	−38.99%	4.289	1%	5.26%	Yes
With duct	−3.05%	−6.60%	−38.07%	4.263	1%	4.55%	Yes
Rd (%)	0.16%	0.72%	1.27%	0.6%

. Additionally, U_G is compared against the error E%D calculated by the following equation

$$E\%D=\frac{D-S_i}{D}, i=1,2,3,$$

(3)

wherein D is the experimental resistance value provided by T2015 [11]. If |U_G| is greater than |E%D| of S₁, the validation is satisfied, i.e., the uncertainty level between CFD and experimental value is below CFD itself. Table 3 shows that our CFD method is also validated.

All CFD resistance results over-predict. The error will reduce to around 3% as grid number increases. The energy saving effect of the duct, i.e., ESD, can be evaluated by the resistance reduction Rd, which is calculated by Equation (4) below:

$$Rd=\frac{C{t}_{with duct}-C{t}_{without duct}}{C{t}_{without duct}}$$

(4)

The Rd result is also shown in Table 3. The medium and fine grid can capture the energy-saving effect better than the coarse grid can. Rd of S₁, S₂ and D is less than 1%, but for S₂ Rd is larger than 1%. In consideration of computational time consumption and flow field resolution, S₂ grid, i.e., medium grid with grid number 2.1 M with duct, is selected for the duct design. Furthermore, S₂ error is over-predicted at around 7%, which is acceptably small, and its Rd is 0.72%, which is much closer to experimental 0.6% than S₁ value (0.16%).

The coarse, medium and fine grid are appended with the optimal 7$°$ duct (the detailed result is in Section 3.3.1) to conduct the grid dependence test. However, the experimental data is not available. The U_D is assumed to be 1%, and the uncertainty is compared with the fine grid result S₁. i.e., Uv%S₁. As shown in

Table 4. Verification analysis for the JBC with the 7° duct.

	S1	S2	S3	RG	UD	Uv%S1
1000Ct	4.310	4.463	5.836	0.111	1%	4.55%
Grid number	6.02 M	2.20 M	0.49 M	Verified		(Validated)

, the grid dependence is confirmed by RG < 1, which means it is verified. The Uv = 4.55%S₁ with the 7$°$ duct coincides with Uv = 4.55%D with the original duct (Table 3). We assume that the S₁ for the 7$°$ duct can provide the similar error level as the S₁ of the original duct does (3.05%). Therefore, the validation can be considered as satisfied.

3.2. Nominal Wake Comparison for the Hull without/with the Original Duct

In the present work, the propeller is not considered, and only bare hull condition is investigated with and without the duct. The nominal wake analysis studies the velocity field on the propeller plane without a propeller. Although there is no propeller effect here, it still helps to understand how the complicated flow field phenomena of ship wake will influence propeller inflow. Once the propeller is put back in place, those phenomena are no longer observable as the flow field is physically blocked and is occupied by the rotating propeller itself. On the other hand, those phenomena will also be affected by propeller effect, which mainly is propeller suction, i.e., propeller induced velocity.

To have the comparison reference, some related values from the T2015 [14] release are listed, and are calculated into non-dimensional results as follows. The radius of the duct outlet is 156.72 mm/7000 mm/2 = 0.0112. The radius of the duct inlet is 111.56 mm/7000 mm/2 = 0.00797. Propeller radius R_P/L = 0.203/2/7 = 0.0145. The propeller center is below the water line at z/L = −0.0404214 and y/L = 0. Propeller plane in longitudinal direction is at x/L = 0.985714; thus, the nominal wake is analyzed on x/L = 0.985714 in our study. The 3-dimensional flow velocity includes u/U for axial direction, v/U for side direction, and w/U for vertical direction. For u/U distribution, i.e., nominal wake factor (1 – w_n) distribution, the contour lines are plotted in Figure 9 for the experimental measurement without and with duct. The flooded contours of u/U are drawn for CFD result without the duct in Figure 10 and with the duct in Figure 11. In both figures, the coarse, medium, and fine grid results which pass the V&V analysis in Section 3.1 (Table 2 and Table 3) are listed side by side, sharing the same contour legend. In those three figures, the vector fields of (v/U, w/U) are topped with u/U for observing vortex behavior.

In the comparison of Figure 9a and Figure 10, the flow field phenomena without the duct are well resolved by CFD results such as the bilge vortex beside the propeller center, low-speed area below the propeller center, and upward flow in the far field. All CFD and experiments show the agreement for the location of bilge vortex to be located between z/L = −0.05 and −0.03 vertically, i.e., almost right next to the propeller center, and around y/L = −0.01 horizontally. The bilge vortex is accompanied by the region of u/U = 0.3–0.5. The contour line of u/U = 0.4 forming a U-shaped pattern surrounds the propeller center on the propeller plane. From the coarse to fine grid, the U-shaped patterns become smoother with a clearer bilge vortex. However, in the experiment, near the core of the bilge vortex u/U reaches 0.2, which cannot be preserved by CFD. It may require a much finer grid to capture the u/U = 0.2 core. The secondary vortex is induced by the bilge vortex inside the low-speed area u/U < 0.3. As the grid densities increase, the pattern of the low-speed area with clearer secondary vortex will also get closer to the experimental pattern. Inside the low-speed area, the medium and fine grid result are both capable of capturing the u/U = 0.2 contour line, and the contour pattern of the fine grid result agrees with the experimental result better. In the far field, the upward flow with the vector component points toward y = 0, due to the ship stern geometry and bilge vortex rotation.

Obviously, the nominal wake would be influenced by the appearance of the original duct comparing Figure 9a to Figure 9b for the experiment, and Figure 10 to Figure 11 for the simulation. In Figure 11, the simulation vector fields are interpolated to 2D uniform spacing grids, which are the same grid used for coarse, medium, and fine grids. This is because the grid is relatively dense around the duct (see Figure 6 and Figure 7) compared with the other part of the flow field. The original vector field will block the view on the propeller plane even in the coarse grid. With the duct effect, the bilge vortex becomes smaller and moves upward to z/L = −0.04~−0.03, but its horizontal location is still around y/L = −0.01 horizontally. The low-speed area below the propeller center disappears and the flow velocity is accelerated evenly to u/U = 0.5–0.6. However, a severe flow separation occurs behind the lower part of the duct, which implies that the duct can be further improved. The CFD captures the above-mentioned flow field characteristics that the experiment measures. The main deviation of CFD is the area of reverse flow inside the flow separation, i.e., u/U < 0 or negative flow velocity, is larger than experimental result. It also explains that the CFD resistance is over-predicted in Table 2. The fine grid predicts a closer flow separation profile showing one isolated area for u/U < 0 and the continuous contour lines for u/U > 0. For the medium grid result, the contour lines of u/U < 0 and u/U = 0.1 are isolated. For the coarse grid result, contour lines of u/U < 0 and u/U = 0.1 are two isolated areas. Unlike the coarse grid, the medium and fine grid can resolve the profile of the experiment-like contour lines of u/U = 0.5–0.8 below the duct. Those contour lines protrude downward along the center plane with the secondary vortex, which is now outside the propeller radius.

Although the finer grid resolves the flow field and agrees with the experiment in more detail, it may reveal new or numerical phenomena. In Figure 10c, for the case without the duct, the u/U = 0.9 contour line of the fine grid result shows a bulge in the far field around y/L = −0.02 and z/L = −0.05 inserting into u/U = 0.8–0.9 area. The similar contour bulge can also be observed for the with-duct result for u/U = 0.8 and 0.9 in Figure 11c (fine grid). Without or with the duct, this does not exist in the experiment, coarse and medium grid (see Figure 9, Figure 10a,b and Figure 11a,b). In the experiment results, y/L = −0.02 and z/L = −0.05 are around the measurement boundary in the portside. Thus, the contour profile in the starboard side shall be considered.

To understand the mechanism of the energy saving and resistance reduction caused by the duct, Figure 12 analyzes the velocity difference between the CFD results with and without the duct on the propeller plane, i.e., nominal wake velocity difference. The medium grid results are selected here, since their Rd is close to the experimental value (Section 3.1 and Table 3). Furthermore, the medium grid is selected for the duct design after the V&V analysis. Table 3 also indicates S_2’s error difference (6.72% − 6.60% = 0.12%) is smaller than S_1’s (3.05% − 2.59% = 0.46%). In addition, the fine grid is not selected here because of the bulge of u/U contour line in the far field, which doesn’t exist in the experiment as discussed previously. On the other hand, the grids without and with the duct are different. Thus, the u/U was first interpolated from the without-duct to with-duct grid, and then subtracted from the u/U with the duct to obtain the velocity difference Δu/U:

Δu/U = u/U_{with duct} − u/U_{no duct(interpolated)}

(5)

The most interesting finding in Figure 12 is the Δu/U = 0–0.1 area. The axial flow velocity is accelerated by the duct not only inside the duct (excluding the area behind the propeller shaft), but also in the far field behind the main hull body. This means the majority of the velocity increases less than 0.1U, but covers a very large area especially behind the ship hull. As a dimension reference, see the left end of the horizontal axis of Figure 12: y/L = −0.08 is around the portside length of the ship beam (0.5B_WL/L = 1.125 m/7 m/2 ~ 0.080). This velocity increases up to 0.2–0.4U in the lower part inside the duct. In the lower position, a band region with 0–0.1U increase is observed in y/L = −0.01~0 and z/L = −0.055~−0.052. Thus, the low-speed area below the propeller center without the duct has been eliminated successfully inside the duct and behind the duct bottom. However, a small area of Δu/U = −0.2~0 is found around z/L = −0.06 near the center plane, implying that the low-speed area remains below the duct bottom. In other words, the low-speed area is not eliminated sufficiently by the duct. A major area of Δu/U < 0 as low as −0.4 corresponds to the bilge vortex location and separation flow behind the duct. This is proof that the nominal wake can be further improved by the duct design in the next section.

The small flow separation occurs after the truncation end of the stern tube. Since it is out of measurement range in experiment, it can only be seen in CFD results. With duct (Figure 12), a clear core is shaped by circular contour lines and low flow speed (u/U < 0.2) at z/L = −0.04 and y/L = 0, i.e., propeller shaft end or propeller center. Without the duct (Figure 11), the core is less clear, connecting with the top of the low-speed area and located slightly below the propeller center. As the grid density increases, the separation core is clearer and the u/U with the duct (Figure 11) is lower than without the duct (Figure 10). This is evidence that the flow separation deteriorates due to the duct. The accelerating duct causes higher velocity passing the truncation end, and the flow separates more severely with lower u/U inside. The Δu/U < 0 core around z/L = −0.04 and y/L = 0 in Figure 13 also illustrates this flow situation.

3.3. Duct Design

3.3.1. Energy Saving Result

In the present work, the JBC energy-saving device (ESD), i.e., stern duct, is designed to obtain the optimal energy saving effect estimated by resistance reduction Rd defined previously in Equation (4). Ct_{without duct} = 4.577 × 10⁻³, which is the S₂ value in Table 2. The design method is CFD simulation using the medium grid for the duct with different AoA of NACA4420 foil section. Table 4 lists the design results of the resistance coefficients, where Ct for the total, i.e., the JBC hull with the duct, Ct_d for the duct only and Ct_h for the hull only. Ct = Ct_h + Ct_d. The trend of Ct, Ct_d, Ct_h, and Rd to different AoA is plotted separately in Figure 13. The 20$°$ AoA is the original design.

The main design process is to decrease the AoA. The trend of Rd, Ct, and Ct_d shown in Figure 14a–c, respectively, is basically consistent. By decreasing AoA, the resistance reduction (Rd) increases correlating to the decrease of the total resistance (Ct). This is because the component of duct resistance Ct_d drops and eventually becomes negative values, i.e., extra thrust. The duct starts producing extra thrust when the AoA < 15$°$. The decreasing trend of Ct, Ct_d, and increasing trend of Rd are significant, but oscillate around 2% when AoA < 10$°$. Since energy saving is our main objective, the optimal AoA is found for 7$°$ to gain the lowest Ct and the highest Rd. In conclusion, the 7$°$ duct design would reduce up to 2.49% resistance, equivalent to a 2.49% energy saving.

There is no sign to possibly increase Rd beyond the minimal and maximal AoA. As AoA increases up to 23$°$ and 24$°$, the duct will no longer reduce the resistance, but instead increase the total resistance. Thus, maximal AoA is set to 24$°$. The minimal 4$°$ is explained near the end of Section 2.2.

The Ct_h disturbs much scatter in Figure 13d, unlike the other trends. Its trend is not clear, but as AoA decreases, it slightly increases. Based on

Table 5. Energy saving result of duct design.

AoA	1000Ct	1000Ctd	1000Cth	Rd (%)
4°	4.486	−0.019	4.505	1.99%
5°	4.475	−0.030	4.505	2.24%
6°	4.484	−0.017	4.501	2.04%
7°	4.463	−0.028	4.491	2.49%
8°	4.480	−0.027	4.507	2.12%
9°	4.476	−0.028	4.504	2.21%
10°	4.479	−0.025	4.504	2.14%
11°	4.479	−0.022	4.501	2.15%
12°	4.483	−0.017	4.500	2.06%
13°	4.494	−0.022	4.516	1.82%
15°	4.512	0.015	4.497	1.43%
16°	4.528	0.025	4.503	1.07%
18°	4.536	0.042	4.494	0.90%
20°	4.544	0.061	4.483	0.72%
21°	4.552	0.072	4.480	0.55%
23°	4.578	0.095	4.483	−0.02%
24°	4.589	0.107	4.482	−0.25%

, all Ct_h are lower than Ct_{without duct} by 1–2%. For AoA lower than 13$°$, Ct_h is almost equal to 4.5 × 10⁻³. For instance, Ct_h = 4.491 × 10⁻³ for the optimal 7$°$ duct is higher than 4.483 × 10⁻³ for original 20$°$ duct. The 7$°$ duct generates the most sufficient extra thrust to overcome the hull resistance increase. For AoA = 15$°$–23$°$, the Rd value is still positive providing the energy saving effect, but the duct itself causes resistance together with the ship hull. For example, the energy saving effect of 20$°$ duct mainly relies on the lower Ct_h instead of the positive Ct_d. This is supported by the analysis for Figure 12 in Section 3.2.

3.3.2. Nominal Wake Analysis

The nominal wake, i.e., axial velocity u/U distribution in ship wake, is analyzed in Figure 14 for several selected duct designs to understand the detailed flow field, energy saving mechanism, and the trend conclusion in Section 3.3.1. The flooded contour legend is the same with Figure 10 and Figure 11. The vector field is not shown here, since the flow field features corresponding to different u/U patterns had been discussed in Section 3.2.

As discussed for Table 3, when the AoA is greater than 15$°$, the duct itself only produces resistance. In Figure 14, the 18$°$, 20$°$ and 23$°$ ducts all have large and severe flow separation behind the lower part of the ducts. Their lowest u/U reaches negative value, i.e., u/U < 0, which is the sign of reverse flow. When the AoA is lower than 15$°$, the duct starts generating the extra thrust. Since the 15$°$ duct is the critical example, its flow separation shows a similar pattern to 18$°$, 20$°$, and 23$°$ ducts with u/U < 0, but it is clearly smaller. Finally, the 7$°$ and 10$°$ ducts no longer have u/U < 0 and much smaller separation. However, for the purpose of accelerating duct, the 7$°$ and 10$°$ ducts perform worse than the 15$°$, 18$°$, 20$°$ and 23$°$ ducts. Inside the radius of the ducts, the u/U can reach 0.5–0.6 for 15$°$, 18$°$, 20$°$ and 23$°$. Instead, inside the 7$°$ and 10$°$ ducts, u/U = 0.5–0.6 only occupy the upper area, and u/U = 0.4–0.5 surrounds the rest of the area. In other words, the flow acceleration of 7$°$ and 10$°$ ducts are lower and less uniform. The u/U profile of outer flow area looks no different among different ducts.

Figure 15 explains how the optimal duct reduces the resistance by overlapping the u/U contour profile of the 7$°$ and 20$°$ duct results in Figure 15a, and then subtracting them in Figure 15b. In Figure 15a, the red contour lines are for the 7$°$ duct and the blue ones for the 20$°$ duct. The velocity difference Δu/U in Figure 15b is calculated as:

$$\mathsf{\Delta }uU=u{U}_{7°\mathrm{duct}}-u{U}_{20°\mathrm{duct\; (interpolated).}}$$

(6)

Since the grid of the 7$°$ and 20$°$ ducts are different, the u/U of 20$°$ duct is interpolated from the 20$°$ duct to 7$°$ duct grid, and then the u/U is subtracted as Equation (6) on the 7$°$ duct grid to obtain Δu/U distribution in Figure 15b.

As observed in Figure 15a, the major flow field improvement due to the 7$°$ duct is the much smaller flow separation. In fact, behind the 7$°$ duct, u/U will not be lower than 0, while u/U < 0, i.e., reverse flow, occupies a considerable large area behind the 20$°$ duct. It is around z/L = −0.05~−0.04 and y/L = −0.01~−0.003 (also see Figure 11b). In the corresponding region for 7$°$ duct, u/U is between 0.1–0.3, and mainly u/U = 0.3 with a very small circle of u/U = 0–0.1. Since the flow velocity is much higher behind the 7$°$ duct in this region, Δu/U > 0.3 (in orange color) can be seen in Figure 14b, surrounded by a large area of Δu/U = 0–0.2 (in green color). The u/U = 0.3 contour line of the 7$°$ duct in Figure 15a extends to the bottom of the duct and downward along the centerline, but for the 20$°$ duct its u/U changes from 0.1 to 0.3. Thus, below the duct bottom, an area with Δu/U = 0–0.2 (in green color) appears in Figure 15b. This means the axial flow velocity below the 7$°$ duct is higher than the 20$°$ duct by 0–0.2U.

However, the 7$°$ duct performs lower axial flow velocity generally corresponding to the Δu/U = −0.2~0 area outside the duct in the outer flow field and inside the duct in Figure 15b. Since the hull with 7$°$ duct produces greater momentum loss behind the hull, it supports the discussion for Ct_h in Table 5 and its trend in Figure 13d in Section 3.3.1, where Ct_h of 7$°$ duct is larger than the 20$°$ duct value. In Figure 15a, in the outer field, the contour lines of 7$°$ duct are located slightly outward than those of 20$°$ duct. For instance, the u/U = 0.9 contour line of the 7$°$ duct is in the area of the u/U = 0.9–1.0. Inside the duct radius, the u/U = 0.5 line of the 7$°$ duct is only in the upper part, while the u/U = 0.5 line of the 20$°$ duct expands to the upper and lower part. As discussed previously, the 7$°$ duct is not a better accelerating duct than the 20$°$ duct. This is the reason the negative Δu/U is observed inside the duct here.

4. Conclusions

The present work has successfully provided an optimal ESD (energy-saving device), i.e., stern duct, for Japan bulk carrier (JBC) at Froude number 0.142 using OpenFOAM as the design tool. Our CFD method is verified and validated for JBC hulls with/without ducts. For verification, the grid independence is confirmed with monotonic convergence. For validation, the uncertainty is ensured larger than the error of the fine grid. In a compromise between time consumption, result resolution and resistance reduction value, the medium grid is selected to design the duct.

The stern duct is axisymmetric along the propeller axis with NACA4420 foil section. In the original duct, the angle of attack (AoA) of the foil section is 20$°$. As AoA decreases, the resistance reduction increases. Finally, the optimal 7$°$ duct is found. The original duct can only reduce 0.6% of ship resistance experimentally, and 0.72% numerically. The 7$°$ duct, however, can achieve the resistance reduction of up to 2.49%. In other words, the energy saving is 2.49% by using the 7$°$ duct. When AoA decreases lower than 15$°$, the duct produces negative resistance, i.e., extra thrust. The hull resistance trend is scattered, but slightly increasing. The 7$°$ duct generates the most sufficient extra thrust to overcome the hull resistance increase. The grid independence is also confirmed for the JBC hull with the optimal 7$°$ duct. By considering the similar uncertainty level, it can be assumed as validated, since the experimental data is not available. Overall, the uncertainty of our CFD method is around 5%. For the error against experiment, fine grid can reach approximately 3% and less than 7% for a medium grid.

For the nominal wake without the duct, three main flow field features will be captured as grid density increases: (1) bilge vortex; (2) low-speed area with the secondary vortex; and (3) upward flow in the far field. Appending the duct, the flow inside the duct will be accelerated with a more uniform velocity distribution. The smaller bilge vortex shifts upward slightly. The low-speed area is eliminated, and the secondary vortex is outside the propeller radius. However, the duct causes a severe flow separation downstream. In fact, the original duct itself produces resistance. The resistance reduction is mainly due to the hull resistance decrease related to the slightly higher axial velocity in the outer wake region. The axial velocity inside the bilge vortex without the duct and the separation flow behind the duct is slightly lower in the CFD results compared to the experiment.

The nominal wake is improved obviously by the optimal 7$°$ duct. Compared to the original duct, the flow separation is much smaller, with no sign of reverse flow. The axial flow velocity is generally much higher behind the 7$°$ duct. However, the 7$°$ duct does not accelerate the axial flow velocity as high and uniform as the duct with the AoA greater than and equal to 15$°$. Additionally, it induces slightly larger hull resistance, corresponding to the slightly slower axial flow velocity in the outer wake region.

In the future work, the duct design will consider the thickness and camber of the foil section, as well as non-axisymmetric geometry. Furthermore, the ship self-propulsion simulation will be conducted with those duct designs to investigate the energy saving effect on the propulsion efficiency. The duct design will also include the propeller effect.